NAMD: Molecule.C Source File

00001 
00007 // Molecule.C is compiled twice!
00008 // MOLECULE2_C undefined only compiles first half of file
00009 // MOLECULE2_C defined only compiles second half of file
00010 // This is shameful but it works.  Molecule needs refactoring badly.
00011 
00012 /*
00013    The class Molecule is used to hold all of the structural information
00014    for a simulation.  This information is read in from a .psf file and
00015    cross checked with the Parameters object passed in.  All of the structural
00016    information is then stored in arrays for use.
00017 */
00018 
00019 #include <stdio.h>
00020 #include <string.h>
00021 #include <stdlib.h>
00022 #include <ctype.h>
00023 
00024 #include "InfoStream.h"
00025 #include "Molecule.h"
00026 #include "strlib.h"
00027 #include "MStream.h"
00028 #include "Communicate.h"
00029 // #include "Node.h"
00030 #include "ObjectArena.h"
00031 #include "Parameters.h"
00032 #include "PDB.h"
00033 #include "SimParameters.h"
00034 #include "Hydrogen.h"
00035 #include "UniqueSetIter.h"
00036 #include "charm++.h"
00037 /* BEGIN gf */
00038 #include "ComputeGridForce.h"
00039 #include "GridForceGrid.h"
00040 
00041 #include "MGridforceParams.h"
00042 /* END gf */
00043 
00044 #define MIN_DEBUG_LEVEL 3
00045 //#define DEBUGM
00046 #include "Debug.h"
00047 
00048 #include "CompressPsf.h"
00049 #include <deque>
00050 #include <algorithm>
00051 
00052 #ifndef M_PI
00053 #define M_PI 3.14159265358979323846
00054 #endif
00055 
00056 using namespace std;
00057 
00058 class ResidueLookupElem
00059 {
00060 public:
00061   char mySegid[11];
00062   ResidueLookupElem *next;      // stored as a linked list
00063   int firstResid;               // valid resid is >= firstResid
00064   int lastResid;                // valid resid is <= lastResid
00065   ResizeArray<int> atomIndex;   // 0-based index for first atom in residue
00066 
00067   ResidueLookupElem(void) { next = 0; firstResid = -1; lastResid = -1; }
00068   ~ResidueLookupElem(void) { delete next; }
00069   int lookup(const char *segid, int resid, int *begin, int *end) const;
00070   ResidueLookupElem* append(const char *segid, int resid, int aid);
00071 };
00072 
00073 #ifndef MOLECULE2_C  // first object file
00074 
00075 #ifdef MEM_OPT_VERSION
00076 template int lookupCstPool<AtomSignature>(const vector<AtomSignature>&, const AtomSignature&);
00077 template int lookupCstPool<ExclusionSignature>(const vector<ExclusionSignature>&, const ExclusionSignature&);
00078 #endif
00079 
00080 int ResidueLookupElem::lookup(
00081         const char *segid, int resid, int *begin, int *end) const {
00082     const ResidueLookupElem *elem = this;
00083     int rval = -1;  // error
00084     while ( elem && strcasecmp(elem->mySegid,segid) ) elem = elem->next;
00085     if ( elem && (resid >= elem->firstResid) && (resid <= elem->lastResid) ) {
00086       *begin = elem->atomIndex[resid - elem->firstResid];
00087       *end = elem->atomIndex[resid - elem->firstResid + 1];
00088       rval = 0;  // no error
00089     }
00090     return rval;
00091 }
00092 
00093 ResidueLookupElem* ResidueLookupElem::append(
00094         const char *segid, int resid, int aid) {
00095     ResidueLookupElem *rval = this;
00096     if ( firstResid == -1 ) {  // nothing added yet
00097       strcpy(mySegid,segid);
00098       firstResid = resid;
00099       lastResid = resid;
00100       atomIndex.add(aid);
00101       atomIndex.add(aid+1);
00102     } else if ( ! strcasecmp(mySegid,segid) ) {  // same segid
00103       if ( resid == lastResid ) {  // same resid
00104         atomIndex[lastResid - firstResid + 1] = aid + 1;
00105       } else if ( resid < lastResid ) {  // error
00106         // We can work around this by creating a new segment.
00107         iout << iWARN << "Residue " << resid <<
00108           " out of order in segment " << segid <<
00109           ", lookup for additional residues in this segment disabled.\n" << endi;
00110         rval = next = new ResidueLookupElem;
00111         next->append(segid,resid,aid);
00112       } else {  // new resid
00113         for ( ; lastResid < resid; ++lastResid ) atomIndex.add(aid);
00114         atomIndex[lastResid - firstResid + 1] = aid + 1;
00115       }
00116     } else {  // new segid
00117       rval = next = new ResidueLookupElem;
00118       next->append(segid,resid,aid);
00119     }
00120     return rval;
00121 }
00122 
00123 
00124 //  Lookup atom id from segment, residue, and name
00125 int Molecule::get_atom_from_name(
00126         const char *segid, int resid, const char *aname) const {
00127 
00128   if (atomNames == NULL || resLookup == NULL)
00129   {
00130     NAMD_die("Tried to find atom from name on node other than node 0");
00131   }
00132 
00133   int i = 0;
00134   int end = 0;
00135   if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
00136   for ( ; i < end; ++i ) {
00137     #ifdef MEM_OPT_VERSION    
00138     Index idx = atomNames[i].atomnameIdx;
00139     if(!strcasecmp(aname, atomNamePool[idx])) return i;
00140     #else
00141     if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i;
00142     #endif
00143   }
00144   return -1;
00145 }
00146 
00147 //  Lookup number of atoms in residue from segment and residue
00148 int Molecule::get_residue_size(
00149         const char *segid, int resid) const {
00150 
00151   if (atomNames == NULL || resLookup == NULL)
00152   {
00153     NAMD_die("Tried to find atom from name on node other than node 0");
00154   }
00155   int i = 0;
00156   int end = 0;
00157   if ( resLookup->lookup(segid,resid,&i,&end) ) return 0;
00158   return ( end - i );
00159 }
00160 
00161 //  Lookup atom id from segment, residue, and index in residue
00162 int Molecule::get_atom_from_index_in_residue(
00163         const char *segid, int resid, int index) const {
00164 
00165   if (atomNames == NULL || resLookup == NULL)
00166   {
00167     NAMD_die("Tried to find atom from name on node other than node 0");
00168   }
00169   int i = 0;
00170   int end = 0;
00171   if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
00172   if ( index >= 0 && index < ( end - i ) ) return ( index + i );
00173   return -1;
00174 }
00175 
00176 /************************************************************************/
00177 /*                  */
00178 /*      FUNCTION initialize  */
00179 /*                  */
00180 /*  This is the initializer for the Molecule class.  It simply sets */
00181 /*  the counts for all the various parameters to 0 and sets the pointers*/
00182 /*  to the arrays that will store these parameters to NULL, since they  */
00183 /*  have not been allocated yet.          */
00184 /*                  */
00185 /************************************************************************/
00186 
00187 void Molecule::initialize(SimParameters *simParams, Parameters *param)
00188 {
00189   if ( sizeof(int32) != 4 ) { NAMD_bug("sizeof(int32) != 4"); }
00190   this->simParams = simParams;
00191   this->params = param;
00192 
00193   /*  Initialize array pointers to NULL  */
00194   atoms=NULL;
00195   atomNames=NULL;
00196   resLookup=NULL;
00197 
00198   // DRUDE
00199   is_drude_psf = 0;  // assume not Drude model
00200   drudeConsts=NULL;
00201   lphosts=NULL;
00202   anisos=NULL;
00203   lphostIndexes=NULL;
00204   // DRUDE
00205 
00206   //for compressing molecule info
00207   atomSegResids=NULL;
00208 
00209   if ( simParams->globalForcesOn ) {
00210     resLookup = new ResidueLookupElem;
00211   }
00212 
00213   #ifdef MEM_OPT_VERSION
00214   eachAtomSig = NULL;
00215   atomSigPoolSize = 0;
00216   atomSigPool = NULL;
00217   massPoolSize = 0;
00218   atomMassPool = NULL;
00219   eachAtomMass = NULL;
00220   chargePoolSize = 0;
00221   atomChargePool = NULL;
00222   eachAtomCharge = NULL;
00223   #else
00224   bonds=NULL;
00225   angles=NULL;
00226   dihedrals=NULL;
00227   impropers=NULL;
00228   crossterms=NULL;
00229   #endif
00230 
00231   donors=NULL;
00232   acceptors=NULL;
00233   
00234 
00235   #ifndef MEM_OPT_VERSION      
00236   tmpArena=NULL;
00237   exclusions=NULL;
00238   bondsWithAtom=NULL;
00239   bondsByAtom=NULL;
00240   anglesByAtom=NULL;
00241   dihedralsByAtom=NULL;
00242   impropersByAtom=NULL;
00243   crosstermsByAtom=NULL;
00244   #endif
00245 
00246   #ifdef MEM_OPT_VERSION
00247   exclSigPool = NULL;
00248   exclChkSigPool = NULL;
00249   exclSigPoolSize = 0;
00250   eachAtomExclSig = NULL;
00251   #else
00252   exclusionsByAtom=NULL;
00253   fullExclusionsByAtom=NULL;
00254   modExclusionsByAtom=NULL;
00255   all_exclusions=NULL;
00256   #endif
00257 
00258   langevinParams=NULL;
00259   fixedAtomFlags=NULL;
00260 
00261   #ifdef MEM_OPT_VERSION  
00262   clusterSigs=NULL;
00263   #else
00264   cluster=NULL;  
00265   #endif
00266   clusterSize=NULL;
00267 
00268   exPressureAtomFlags=NULL;
00269   rigidBondLengths=NULL;
00270   consIndexes=NULL;
00271   consParams=NULL;
00272   /* BEGIN gf */
00273   gridfrcIndexes=NULL;
00274   gridfrcParams=NULL;
00275   gridfrcGrid=NULL;
00276   numGridforces=NULL;
00277   /* END gf */
00278   stirIndexes=NULL;
00279   stirParams=NULL;
00280   movDragIndexes=NULL;
00281   movDragParams=NULL;
00282   rotDragIndexes=NULL;
00283   rotDragParams=NULL;
00284   consTorqueIndexes=NULL;
00285   consTorqueParams=NULL;
00286   consForceIndexes=NULL;
00287   consForce=NULL;
00288 //fepb
00289   fepAtomFlags=NULL;
00290 //fepe
00291 
00292   nameArena = new ObjectArena<char>;
00293   // nameArena->setAlignment(8);
00294   // arena->setAlignment(32);
00295   #ifndef MEM_OPT_VERSION
00296   arena = new ObjectArena<int32>;
00297   exclArena = new ObjectArena<char>;
00298   #endif
00299   // exclArena->setAlignment(32);
00300 
00301   /*  Initialize counts to 0 */
00302   numAtoms=0;
00303   numRealBonds=0;
00304   numBonds=0;
00305   numAngles=0;
00306   numDihedrals=0;
00307   numImpropers=0;
00308   numCrossterms=0;
00309   numDonors=0;
00310   numAcceptors=0;
00311   numExclusions=0;
00312 
00313   // DRUDE
00314   numLonepairs=0;
00315   numDrudeAtoms=0;
00316   numLphosts=0;
00317   numAnisos=0;
00318   // DRUDE
00319 
00320   numConstraints=0;
00321   numStirredAtoms=0;
00322   numMovDrag=0;
00323   numRotDrag=0;
00324   numConsTorque=0;
00325   numConsForce=0;
00326   numFixedAtoms=0;
00327   numFixedGroups=0;
00328   numExPressureAtoms=0;
00329   numRigidBonds=0;
00330   numFixedRigidBonds=0;
00331   numMultipleDihedrals=0;
00332   numMultipleImpropers=0;
00333   numCalcBonds=0;
00334   numCalcAngles=0;
00335   numCalcDihedrals=0;
00336   numCalcImpropers=0;
00337   numCalcCrossterms=0;
00338   numCalcExclusions=0;
00339 
00340 //fepb
00341   numFepInitial = 0;
00342   numFepFinal = 0;
00343 //fepe
00344 
00345   //fields related with pluginIO-based loading molecule structure
00346   occupancy = NULL;
00347   bfactor = NULL;
00348 }
00349 
00350 /*      END OF FUNCTION initialize */
00351 
00352 /************************************************************************/
00353 /*                  */
00354 /*      FUNCTION Molecule        */
00355 /*                  */
00356 /*  This is the constructor for the Molecule class. */
00357 /*                  */
00358 /************************************************************************/
00359 
00360 Molecule::Molecule(SimParameters *simParams, Parameters *param)
00361 {
00362   initialize(simParams,param);
00363 }
00364 
00365 /************************************************************************/
00366 /*                  */
00367 /*      FUNCTION Molecule        */
00368 /*                  */
00369 /*  This is the constructor for the Molecule class from CHARMM/XPLOR files. */
00370 /*                  */
00371 /************************************************************************/
00372 
00373 Molecule::Molecule(SimParameters *simParams, Parameters *param, char *filename, ConfigList *cfgList)
00374 {
00375   initialize(simParams,param);
00376 
00377   if(simParams->useCompressedPsf)
00378       read_compressed_psf_file(filename, param, cfgList);
00379   /*else if(simParams->genCompressedPsf){      
00380       compress_psf_file(this, filename, param, simParams, cfgList);
00381   }*/      
00382   else
00383       read_psf_file(filename, param);
00384 }
00385 
00386 /************************************************************************/
00387 /*                                                                      */
00388 /*      FUNCTION Molecule                                               */
00389 /*                                                                      */
00390 /*  This is the constructor for the Molecule class from plugin IO.      */
00391 /*                                                                      */
00392 /************************************************************************/
00393 Molecule::Molecule(SimParameters *simParams, Parameters *param, molfile_plugin_t *pIOHdl, void *pIOFileHdl, int natoms)
00394 {
00395 #ifdef MEM_OPT_VERSION
00396   NAMD_die("Sorry, plugin IO is not supported in the memory optimized version.");
00397 #endif
00398     initialize(simParams, param);
00399     numAtoms = natoms;
00400     int optflags = MOLFILE_BADOPTIONS;
00401     molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*sizeof(molfile_atom_t));
00402     memset(atomarray, 0, natoms*sizeof(molfile_atom_t));
00403 
00404     //1a. read basic atoms information
00405     int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
00406     if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
00407         free(atomarray);
00408         NAMD_die("ERROR: plugin failed reading structure data");
00409     }
00410     if(optflags == MOLFILE_BADOPTIONS) {
00411         free(atomarray);
00412         NAMD_die("ERROR: plugin didn't initialize optional data flags");
00413     }
00414     if(optflags & MOLFILE_OCCUPANCY) {
00415         setOccupancyData(atomarray);
00416     }
00417     if(optflags & MOLFILE_BFACTOR) {
00418         setBFactorData(atomarray);
00419     }
00420     //1b. load basic atoms information to the molecule object
00421     plgLoadAtomBasics(atomarray);    
00422     free(atomarray);
00423 
00424     //2a. read bonds
00425     //indices are one-based in read_bonds
00426     int *from, *to;
00427     float *bondorder;
00428     if(pIOHdl->read_bonds!=NULL) {
00429         if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder)){
00430             NAMD_die("ERROR: failed reading bond information.");
00431         }
00432     }    
00433     //2b. load bonds information to the molecule object
00434     if(numBonds!=0) {
00435         plgLoadBonds(from,to);
00436     }
00437 
00438     //3a. read other bonded structures
00439     int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
00440     int ctermcols, ctermrows;
00441     double *angleforces,  *dihedralforces, *improperforces, *ctermforces;
00442 
00443     plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
00444     if(pIOHdl->read_angles!=NULL) {
00445         if(pIOHdl->read_angles(pIOFileHdl,
00446                                &numAngles, &plgAngles, &angleforces,
00447                                &numDihedrals, &plgDihedrals, &dihedralforces,
00448                                &numImpropers, &plgImpropers, &improperforces,
00449                                &numCrossterms, &plgCterms, &ctermcols, &ctermrows,
00450                                &ctermforces)) {
00451             NAMD_die("ERROR: failed reading angle information.");
00452         }
00453     }
00454     //3b. load other bonded structures to the molecule object
00455     if(numAngles!=0) plgLoadAngles(plgAngles);
00456     if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
00457     if(numImpropers!=0) plgLoadImpropers(plgImpropers);
00458     if(numCrossterms!=0) plgLoadCrossterms(plgCterms);
00459 
00460   numRealBonds = numBonds;
00461   build_atom_status();
00462 
00463 }
00464 
00465 /*      END OF FUNCTION Molecule      */
00466 
00467 /************************************************************************/
00468 /*                  */
00469 /*        FUNCTION Molecule      */
00470 /*                  */
00471 /*  This is the destructor for the class Molecule.  It simply frees */
00472 /*  the memory allocated for each of the arrays used to store the       */
00473 /*  structure information.            */
00474 /*                  */
00475 /************************************************************************/
00476 
00477 Molecule::~Molecule()
00478 {
00479   /*  Check to see if each array was ever allocated.  If it was   */
00480   /*  then free it            */
00481   if (atoms != NULL)
00482     delete [] atoms;
00483 
00484   if (atomNames != NULL)
00485   {
00486     // subarrarys allocated from arena - automatically deleted
00487     delete [] atomNames;
00488   }
00489   delete nameArena;
00490 
00491   if (resLookup != NULL)
00492     delete resLookup;
00493 
00494   // DRUDE: free arrays read from PSF
00495   if (drudeConsts != NULL) delete [] drudeConsts;
00496   if (lphosts != NULL) delete [] lphosts;
00497   if (anisos != NULL) delete [] anisos;
00498   if (lphostIndexes != NULL) delete [] lphostIndexes;
00499   // DRUDE
00500 
00501   #ifdef MEM_OPT_VERSION
00502   if(eachAtomSig) delete [] eachAtomSig;
00503   if(atomSigPool) delete [] atomSigPool;
00504   #else
00505   if (bonds != NULL)
00506     delete [] bonds;
00507 
00508   if (angles != NULL)
00509     delete [] angles;
00510 
00511   if (dihedrals != NULL)
00512     delete [] dihedrals;
00513 
00514   if (impropers != NULL)
00515     delete [] impropers;
00516 
00517   if (crossterms != NULL)
00518     delete [] crossterms;
00519 
00520   if (exclusions != NULL)
00521     delete [] exclusions;
00522   #endif
00523 
00524   if (donors != NULL)
00525     delete [] donors;
00526 
00527   if (acceptors != NULL)
00528     delete [] acceptors;  
00529 
00530   #ifdef MEM_OPT_VERSION
00531   if(exclSigPool) delete [] exclSigPool;
00532   if(exclChkSigPool) delete [] exclChkSigPool;
00533   if(eachAtomExclSig) delete [] eachAtomExclSig;
00534   #else
00535   if (bondsByAtom != NULL)
00536        delete [] bondsByAtom;
00537   
00538   if (anglesByAtom != NULL)
00539        delete [] anglesByAtom;
00540   
00541   if (dihedralsByAtom != NULL)
00542        delete [] dihedralsByAtom;
00543   
00544   if (impropersByAtom != NULL)
00545        delete [] impropersByAtom;
00546   
00547   if (crosstermsByAtom != NULL)
00548        delete [] crosstermsByAtom;  
00549 
00550   if (exclusionsByAtom != NULL)
00551        delete [] exclusionsByAtom;
00552   
00553   if (fullExclusionsByAtom != NULL)
00554        delete [] fullExclusionsByAtom;
00555   
00556   if (modExclusionsByAtom != NULL)
00557        delete [] modExclusionsByAtom;
00558   
00559   if (all_exclusions != NULL)
00560        delete [] all_exclusions;
00561   #endif
00562 
00563 
00564   if (fixedAtomFlags != NULL)
00565        delete [] fixedAtomFlags;
00566 
00567   if (stirIndexes != NULL)
00568     delete [] stirIndexes;
00569 
00570 
00571   #ifdef MEM_OPT_VERSION
00572   if(clusterSigs != NULL){      
00573       delete [] clusterSigs;
00574   }  
00575   #else
00576   if (cluster != NULL)
00577        delete [] cluster;  
00578   #endif
00579   if (clusterSize != NULL)
00580        delete [] clusterSize;
00581 
00582   if (exPressureAtomFlags != NULL)
00583        delete [] exPressureAtomFlags;
00584 
00585   if (rigidBondLengths != NULL)
00586        delete [] rigidBondLengths;
00587 
00588 //fepb
00589   if (fepAtomFlags != NULL)
00590        delete [] fepAtomFlags;
00591 //fepe
00592 
00593 
00594   #ifndef MEM_OPT_VERSION
00595   delete arena;
00596   delete exclArena;
00597   #endif
00598 }
00599 /*      END OF FUNCTION Molecule      */
00600 
00601 /************************************************************************/
00602 /*                  */
00603 /*        FUNCTION read_psf_file      */
00604 /*                  */
00605 /*   INPUTS:                */
00606 /*  fname - Name of the .psf file to read        */
00607 /*  params - pointer to Parameters object to use to obtain          */
00608 /*     parameters for vdWs, bonds, etc.      */
00609 /*                  */
00610 /*  This function reads a .psf file in.  This is where just about   */
00611 /*   all of the structural information for this class comes from.  The  */
00612 /*   .psf file contains descriptions of the atom, bonds, angles,        */
00613 /*   dihedrals, impropers, and exclusions.  The parameter object is     */
00614 /*   used to look up parameters for each of these entities.    */
00615 /*                  */
00616 /************************************************************************/
00617 
00618 void Molecule::read_psf_file(char *fname, Parameters *params)
00619 
00620 {
00621 #ifdef MEM_OPT_VERSION
00622     return;
00623 #else
00624   char err_msg[512];  //  Error message for NAMD_die
00625   char buffer[512];  //  Buffer for file reading
00626   int i;      //  Loop counter
00627   int NumTitle;    //  Number of Title lines in .psf file
00628   FILE *psf_file;    //  pointer to .psf file
00629   int ret_code;    //  ret_code from NAMD_read_line calls
00630 
00631   /* Try and open the .psf file           */
00632   if ( (psf_file = Fopen(fname, "r")) == NULL)
00633   {
00634     sprintf(err_msg, "UNABLE TO OPEN .psf FILE %s", fname);
00635     NAMD_die(err_msg);
00636   }
00637 
00638   /*  Read till we have the first non-blank line of file    */
00639   ret_code = NAMD_read_line(psf_file, buffer);
00640 
00641   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00642   {
00643     ret_code = NAMD_read_line(psf_file, buffer);
00644   }
00645 
00646   /*  Check to see if we dropped out of the loop because of a     */
00647   /*  read error.  This shouldn't happen unless the file is empty */
00648   if (ret_code!=0)
00649   {
00650     sprintf(err_msg, "EMPTY .psf FILE %s", fname);
00651     NAMD_die(err_msg);
00652   }
00653 
00654   /*  The first non-blank line should contain the word "psf".    */
00655   /*  If we can't find it, die.               */
00656   if (!NAMD_find_word(buffer, "psf"))
00657   {
00658     sprintf(err_msg, "UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
00659        fname);
00660     NAMD_die(err_msg);
00661   }
00662 
00663   // DRUDE: set flag if we discover Drude PSF
00664   if (NAMD_find_word(buffer, "drude"))
00665   {
00666     is_drude_psf = 1;
00667   }
00668   // DRUDE
00669 
00670   /*  Read until we find the next non-blank line      */
00671   ret_code = NAMD_read_line(psf_file, buffer);
00672 
00673   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00674   {
00675     ret_code = NAMD_read_line(psf_file, buffer);
00676   }
00677 
00678   /*  Check to see if we dropped out of the loop because of a     */
00679   /*  read error.  This shouldn't happen unless there is nothing  */
00680   /*  but the PSF line in the file        */
00681   if (ret_code!=0)
00682   {
00683     sprintf(err_msg, "MISSING EVERYTHING BUT PSF FROM %s", fname);
00684     NAMD_die(err_msg);
00685   }
00686 
00687   /*  This line should have the word "NTITLE" in it specifying    */
00688   /*  how many title lines there are        */
00689   if (!NAMD_find_word(buffer, "NTITLE"))
00690   {
00691     sprintf(err_msg,"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
00692        fname);
00693     NAMD_die(err_msg);
00694   }
00695 
00696   sscanf(buffer, "%d", &NumTitle);
00697 
00698   /*  Now skip the next NTITLE non-blank lines and then read in the*/
00699   /*  line which should contain NATOM        */
00700   i=0;
00701 
00702   while ( ((ret_code=NAMD_read_line(psf_file, buffer)) == 0) && 
00703     (i<NumTitle) )
00704   {
00705     if (!NAMD_blank_string(buffer))
00706       i++;
00707   }
00708 
00709   /*  Make sure we didn't exit because of a read error    */
00710   if (ret_code!=0)
00711   {
00712     sprintf(err_msg, "FOUND EOF INSTEAD OF NATOM IN PSF FILE %s", 
00713        fname);
00714     NAMD_die(err_msg);
00715   }
00716 
00717   while (NAMD_blank_string(buffer))
00718   {
00719     NAMD_read_line(psf_file, buffer);
00720   }
00721 
00722   /*  Check to make sure we have the line we want      */
00723   if (!NAMD_find_word(buffer, "NATOM"))
00724   {
00725     sprintf(err_msg, "DIDN'T FIND \"NATOM\" IN PSF FILE %s",
00726        fname);
00727     NAMD_die(err_msg);
00728   }
00729 
00730   /*  Read in the number of atoms, and then the atoms themselves  */
00731   sscanf(buffer, "%d", &numAtoms);
00732 
00733   read_atoms(psf_file, params);
00734 
00735   /*  Read until we find the next non-blank line      */
00736   ret_code = NAMD_read_line(psf_file, buffer);
00737 
00738   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00739   {
00740     ret_code = NAMD_read_line(psf_file, buffer);
00741   }
00742 
00743   /*  Check to make sure we didn't hit the EOF      */
00744   if (ret_code != 0)
00745   {
00746     NAMD_die("EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
00747   }
00748 
00749   /*  Look for the string "NBOND"          */
00750   if (!NAMD_find_word(buffer, "NBOND"))
00751   {
00752     NAMD_die("DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
00753   }
00754 
00755   /*  Read in the number of bonds and then the bonds themselves  */
00756   sscanf(buffer, "%d", &numBonds);
00757 
00758   if (numBonds)
00759     read_bonds(psf_file, params);
00760 
00761   /*  Read until we find the next non-blank line      */
00762   ret_code = NAMD_read_line(psf_file, buffer);
00763 
00764   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00765   {
00766     ret_code = NAMD_read_line(psf_file, buffer);
00767   }
00768 
00769   /*  Check to make sure we didn't hit the EOF      */
00770   if (ret_code != 0)
00771   {
00772     NAMD_die("EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
00773   }
00774 
00775   /*  Look for the string "NTHETA"        */
00776   if (!NAMD_find_word(buffer, "NTHETA"))
00777   {
00778     NAMD_die("DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
00779   }
00780 
00781   /*  Read in the number of angles and then the angles themselves */
00782   sscanf(buffer, "%d", &numAngles);
00783 
00784   if (numAngles)
00785     read_angles(psf_file, params);
00786 
00787   /*  Read until we find the next non-blank line      */
00788   ret_code = NAMD_read_line(psf_file, buffer);
00789 
00790   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00791   {
00792     ret_code = NAMD_read_line(psf_file, buffer);
00793   }
00794 
00795   /*  Check to make sure we didn't hit the EOF      */
00796   if (ret_code != 0)
00797   {
00798     NAMD_die("EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
00799   }
00800 
00801   /*  Look for the string "NPHI"          */
00802   if (!NAMD_find_word(buffer, "NPHI"))
00803   {
00804     NAMD_die("DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
00805   }
00806 
00807   /*  Read in the number of dihedrals and then the dihedrals      */
00808   sscanf(buffer, "%d", &numDihedrals);
00809 
00810   if (numDihedrals)
00811     read_dihedrals(psf_file, params);
00812 
00813   /*  Read until we find the next non-blank line      */
00814   ret_code = NAMD_read_line(psf_file, buffer);
00815 
00816   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00817   {
00818     ret_code = NAMD_read_line(psf_file, buffer);
00819   }
00820 
00821   /*  Check to make sure we didn't hit the EOF      */
00822   if (ret_code != 0)
00823   {
00824     NAMD_die("EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
00825   }
00826 
00827   /*  Look for the string "NIMPHI"        */
00828   if (!NAMD_find_word(buffer, "NIMPHI"))
00829   {
00830     NAMD_die("DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
00831   }
00832 
00833   /*  Read in the number of Impropers and then the impropers  */
00834   sscanf(buffer, "%d", &numImpropers);
00835 
00836   if (numImpropers)
00837     read_impropers(psf_file, params);
00838 
00839   /*  Read until we find the next non-blank line      */
00840   ret_code = NAMD_read_line(psf_file, buffer);
00841 
00842   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00843   {
00844     ret_code = NAMD_read_line(psf_file, buffer);
00845   }
00846 
00847   /*  Check to make sure we didn't hit the EOF      */
00848   if (ret_code != 0)
00849   {
00850     NAMD_die("EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
00851   }
00852 
00853   /*  Look for the string "NDON"        */
00854   if (!NAMD_find_word(buffer, "NDON"))
00855   {
00856     NAMD_die("DID NOT FIND NDON AFTER ATOM LIST IN PSF");
00857   }
00858 
00859   /*  Read in the number of hydrogen bond donors and then the donors */
00860   sscanf(buffer, "%d", &numDonors);
00861 
00862   if (numDonors)
00863     read_donors(psf_file);
00864 
00865   /*  Read until we find the next non-blank line      */
00866   ret_code = NAMD_read_line(psf_file, buffer);
00867 
00868   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00869   {
00870     ret_code = NAMD_read_line(psf_file, buffer);
00871   }
00872 
00873   /*  Check to make sure we didn't hit the EOF      */
00874   if (ret_code != 0)
00875   {
00876     NAMD_die("EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
00877   }
00878 
00879   /*  Look for the string "NACC"        */
00880   if (!NAMD_find_word(buffer, "NACC"))
00881   {
00882     NAMD_die("DID NOT FIND NACC AFTER ATOM LIST IN PSF");
00883   }
00884 
00885   /*  Read in the number of hydrogen bond donors and then the donors */
00886   sscanf(buffer, "%d", &numAcceptors);
00887 
00888   if (numAcceptors)
00889     read_acceptors(psf_file);
00890 
00891   /*  look for the explicit non-bonded exclusion section.     */
00892   while (!NAMD_find_word(buffer, "NNB"))
00893   {
00894     ret_code = NAMD_read_line(psf_file, buffer);
00895 
00896     if (ret_code != 0)
00897     {
00898       NAMD_die("EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
00899     }
00900   }
00901 
00902   /*  Read in the number of exclusions and then the exclusions    */
00903   sscanf(buffer, "%d", &numExclusions);
00904 
00905   if (numExclusions)
00906     read_exclusions(psf_file);
00907 
00908   // DRUDE: read lone pair hosts and anisotropic terms from PSF
00909   if (is_drude_psf)
00910   {
00911     while (!NAMD_find_word(buffer, "NUMLP"))
00912     {
00913       ret_code = NAMD_read_line(psf_file, buffer);
00914       if (ret_code != 0)
00915       {
00916         NAMD_die("EOF ENCOUNTERED LOOKING FOR NUMLP IN DRUDE PSF FILE");
00917       }
00918     }
00919     sscanf(buffer, "%d", &numLphosts);
00920     if (numLphosts) read_lphosts(psf_file);
00921 
00922     while (!NAMD_find_word(buffer, "NUMANISO"))
00923     {
00924       ret_code = NAMD_read_line(psf_file, buffer);
00925       if (ret_code != 0)
00926       {
00927         NAMD_die("EOF ENCOUNTERED LOOKING FOR NUMANISO IN DRUDE PSF FILE");
00928       }
00929     }
00930     sscanf(buffer, "%d", &numAnisos);
00931     if (numAnisos) read_anisos(psf_file);
00932 
00933   }
00934   // DRUDE
00935 
00936   /*  look for the cross-term section.     */
00937   int crossterms_present = 1;
00938   while (!NAMD_find_word(buffer, "NCRTERM"))
00939   {
00940     ret_code = NAMD_read_line(psf_file, buffer);
00941 
00942     if (ret_code != 0)
00943     {
00944       // hit EOF before finding cross-term section
00945       crossterms_present = 0;
00946       break;
00947     }
00948   }
00949 
00950   if ( crossterms_present) {
00951 
00952     /*  Read in the number of cross-terms and then the cross-terms*/
00953     sscanf(buffer, "%d", &numCrossterms);
00954 
00955     if (numCrossterms)
00956       read_crossterms(psf_file, params);
00957 
00958   }
00959 
00960   /*  Close the .psf file.  */
00961   Fclose(psf_file);
00962 
00963   //  analyze the data and find the status of each atom
00964   numRealBonds = numBonds;
00965   build_atom_status();
00966 
00967   return;
00968 #endif
00969 }
00970 /*      END OF FUNCTION read_psf_file      */
00971 
00972 /************************************************************************/
00973 /*                  */
00974 /*        FUNCTION read_compressed_psf_file      */
00975 /*                  */
00976 /*   INPUTS:                */
00977 /*  fname - Name of the compressed .psf file to read        */
00978 /*  params - pointer to Parameters object to use to obtain          */
00979 /*     parameters for vdWs, bonds, etc.      */
00980 /*                  */
00981 /*  This function reads a compressed .psf file in.  This is where just about   */
00982 /*   all of the structural information for this class comes from.  The  */
00983 /*   .psf file contains descriptions of the atom, bonds, angles,        */
00984 /*   dihedrals, impropers, and exclusions.  The parameter object is     */
00985 /*   used to look up parameters for each of these entities.    */
00986 /*                  */
00987 /************************************************************************/
00988 void Molecule::read_compressed_psf_file(char *fname, Parameters *params, ConfigList *cfgList){
00989 #ifndef MEM_OPT_VERSION
00990     return;
00991 #else
00992     FILE *psf_file;    //  pointer to .psf file
00993     int ret_code;    //  ret_code from NAMD_read_line calls
00994     char buffer[512];
00995 
00996     /* Try and open the .psf file           */
00997     if ( (psf_file = Fopen(fname, "r")) == NULL)
00998     {
00999         char err_msg[512];
01000         sprintf(err_msg, "UNABLE TO OPEN THE COMPRESSED .psf FILE %s", fname);
01001         NAMD_die(err_msg);
01002     }     
01003 
01004     char strBuf[12];
01005 
01006     //first check compressed psf file format version
01007     NAMD_read_line(psf_file, buffer);
01008     if(!NAMD_find_word(buffer, "FORMAT VERSION")) {
01009         NAMD_die("The compressed psf file format is incorrect, please re-generate!\n");
01010     }
01011     float psfVer = 0.0f;
01012     sscanf(buffer, "FORMAT VERSION: %f\n", &psfVer);
01013     if(fabs(psfVer - COMPRESSED_PSF_VER)>1e-6) {
01014         NAMD_die("The compressed psf file format is incorrect, please re-generate!\n");
01015     }
01016 
01017     //Begin reading constant pools for atoms' basic information
01018     //1. segment names
01019     NAMD_read_line(psf_file, buffer);
01020     if(!NAMD_find_word(buffer, "NSEGMENTNAMES"))
01021         NAMD_die("UNABLE TO FIND NSEGMENTNAMES");
01022     sscanf(buffer, "%d", &segNamePoolSize);
01023     if(segNamePoolSize!=0)
01024         segNamePool = new char *[segNamePoolSize];
01025     for(int i=0; i<segNamePoolSize; i++){
01026         NAMD_read_line(psf_file, buffer);
01027         sscanf(buffer, "%s", strBuf);
01028         segNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
01029         strcpy(segNamePool[i], strBuf);        
01030     }    
01031 
01032     //2. residue names
01033     NAMD_read_line(psf_file, buffer);
01034     if(!NAMD_find_word(buffer, "NRESIDUENAMES"))
01035         NAMD_die("UNABLE TO FIND NRESIDUENAMES");
01036     sscanf(buffer, "%d", &resNamePoolSize);
01037     if(resNamePoolSize!=0)
01038         resNamePool = new char *[resNamePoolSize];
01039     for(int i=0; i<resNamePoolSize; i++){
01040         NAMD_read_line(psf_file, buffer);
01041         sscanf(buffer, "%s", strBuf);
01042         resNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
01043         strcpy(resNamePool[i], strBuf);
01044     }
01045 
01046     //3. atom names
01047     NAMD_read_line(psf_file, buffer);
01048     if(!NAMD_find_word(buffer, "NATOMNAMES"))
01049         NAMD_die("UNABLE TO FIND NATOMNAMES");
01050     sscanf(buffer, "%d", &atomNamePoolSize);
01051     if(atomNamePoolSize!=0)
01052         atomNamePool = new char *[atomNamePoolSize];
01053     for(int i=0; i<atomNamePoolSize; i++){
01054         NAMD_read_line(psf_file, buffer);
01055         sscanf(buffer, "%s", strBuf);
01056         atomNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
01057         strcpy(atomNamePool[i], strBuf);
01058     }
01059 
01060     //4. atom types
01061     NAMD_read_line(psf_file, buffer);
01062     if(!NAMD_find_word(buffer, "NATOMTYPES"))
01063         NAMD_die("UNABLE TO FIND NATOMTYPES");
01064     sscanf(buffer, "%d", &atomTypePoolSize);
01065     if(atomTypePoolSize!=0)
01066         atomTypePool = new char *[atomTypePoolSize];
01067     for(int i=0; i<atomTypePoolSize; i++){
01068         NAMD_read_line(psf_file, buffer);
01069         sscanf(buffer, "%s", strBuf);
01070         atomTypePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
01071         strcpy(atomTypePool[i], strBuf);
01072     }
01073 
01074     //5. charges
01075     NAMD_read_line(psf_file, buffer);
01076     if(!NAMD_find_word(buffer, "NCHARGES"))
01077         NAMD_die("UNABLE TO FIND NCHARGES");
01078     sscanf(buffer, "%d", &chargePoolSize);
01079     if(chargePoolSize!=0)
01080         atomChargePool = new Real[chargePoolSize];
01081     for(int i=0; i<chargePoolSize; i++){
01082         NAMD_read_line(psf_file, buffer);
01083         sscanf(buffer, "%f", atomChargePool+i);        
01084     }
01085 
01086     //6. masses
01087     NAMD_read_line(psf_file, buffer);
01088     if(!NAMD_find_word(buffer, "NMASSES"))
01089         NAMD_die("UNABLE TO FIND NMASSES");
01090     sscanf(buffer, "%d", &massPoolSize);
01091     if(massPoolSize!=0)
01092         atomMassPool = new Real[massPoolSize];
01093     for(int i=0; i<massPoolSize; i++){
01094         NAMD_read_line(psf_file, buffer);
01095         sscanf(buffer, "%f", atomMassPool+i);        
01096     }
01097 
01098     //Begin reading of atom signatures
01099     NAMD_read_line(psf_file, buffer);
01100     if(!NAMD_find_word(buffer, "ATOMSIGS"))
01101         NAMD_die("UNABLE TO FIND ATOMSIGS");
01102     sscanf(buffer, "%d", &atomSigPoolSize);
01103     atomSigPool = new AtomSignature[atomSigPoolSize];
01104     int typeCnt;
01105     int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
01106     int tisReal;
01107     int ttype;
01108     for(int i=0; i<atomSigPoolSize; i++){
01109         //BOND SIGS
01110         NAMD_read_line(psf_file, buffer);
01111         if(!NAMD_find_word(buffer, "NBONDSIGS"))
01112             NAMD_die("UNABLE TO FIND NBONDSIGS");
01113         sscanf(buffer, "%d", &typeCnt);
01114         if(typeCnt!=0){
01115             atomSigPool[i].bondCnt = typeCnt;
01116             atomSigPool[i].bondSigs = new TupleSignature[typeCnt];
01117         }            
01118         for(int j=0; j<typeCnt; j++){
01119             NAMD_read_line(psf_file, buffer);            
01120             sscanf(buffer, "%d | %d | %d", &tmp1, &ttype, &tisReal);
01121             TupleSignature oneSig(1, BOND, (Index)ttype, (char)tisReal);
01122             oneSig.offset[0] = tmp1;            
01123             atomSigPool[i].bondSigs[j]=oneSig;
01124             if(tisReal) numRealBonds++;
01125         }
01126 
01127         //ANGLE SIGS
01128         NAMD_read_line(psf_file, buffer);
01129         if(!NAMD_find_word(buffer, "NTHETASIGS"))
01130             NAMD_die("UNABLE TO FIND NTHETASIGS");
01131         sscanf(buffer, "%d", &typeCnt);
01132         if(typeCnt!=0){
01133             atomSigPool[i].angleCnt = typeCnt;
01134             atomSigPool[i].angleSigs = new TupleSignature[typeCnt];
01135         }            
01136         for(int j=0; j<typeCnt; j++){
01137             NAMD_read_line(psf_file, buffer);            
01138             sscanf(buffer, "%d %d | %d | %d", &tmp1, &tmp2, &ttype, &tisReal);
01139             TupleSignature oneSig(2,ANGLE,(Index)ttype, (char)tisReal);
01140             oneSig.offset[0] = tmp1;
01141             oneSig.offset[1] = tmp2;            
01142             atomSigPool[i].angleSigs[j] = oneSig;
01143         }
01144         //DIHEDRAL SIGS
01145         NAMD_read_line(psf_file, buffer);
01146         if(!NAMD_find_word(buffer, "NPHISIGS"))
01147             NAMD_die("UNABLE TO FIND NPHISIGS");
01148         sscanf(buffer, "%d", &typeCnt);
01149         if(typeCnt!=0){
01150             atomSigPool[i].dihedralCnt = typeCnt;
01151             atomSigPool[i].dihedralSigs = new TupleSignature[typeCnt];
01152         }            
01153         for(int j=0; j<typeCnt; j++){
01154             NAMD_read_line(psf_file, buffer);            
01155             sscanf(buffer, "%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
01156             TupleSignature oneSig(3,DIHEDRAL,(Index)ttype, (char)tisReal);
01157             oneSig.offset[0] = tmp1;
01158             oneSig.offset[1] = tmp2;
01159             oneSig.offset[2] = tmp3;
01160             atomSigPool[i].dihedralSigs[j] = oneSig;
01161         }
01162         //IMPROPER SIGS
01163         NAMD_read_line(psf_file, buffer);
01164         if(!NAMD_find_word(buffer, "NIMPHISIGS"))
01165             NAMD_die("UNABLE TO FIND NIMPHISIGS");
01166         sscanf(buffer, "%d", &typeCnt);
01167         if(typeCnt!=0){
01168             atomSigPool[i].improperCnt = typeCnt;
01169             atomSigPool[i].improperSigs = new TupleSignature[typeCnt];
01170         }            
01171         for(int j=0; j<typeCnt; j++){
01172             NAMD_read_line(psf_file, buffer);            
01173             sscanf(buffer, "%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
01174             TupleSignature oneSig(3,IMPROPER,(Index)ttype, (char)tisReal);
01175             oneSig.offset[0] = tmp1;
01176             oneSig.offset[1] = tmp2;
01177             oneSig.offset[2] = tmp3;
01178             atomSigPool[i].improperSigs[j] = oneSig;
01179         }        
01180         //CROSSTERM SIGS
01181         NAMD_read_line(psf_file, buffer);
01182         if(!NAMD_find_word(buffer, "NCRTERMSIGS"))
01183             NAMD_die("UNABLE TO FIND NCRTERMSIGS");
01184         sscanf(buffer, "%d", &typeCnt);
01185         if(typeCnt!=0){
01186             atomSigPool[i].crosstermCnt = typeCnt;
01187             atomSigPool[i].crosstermSigs = new TupleSignature[typeCnt];
01188         }            
01189         for(int j=0; j<typeCnt; j++){
01190             NAMD_read_line(psf_file, buffer);            
01191             sscanf(buffer, "%d %d %d %d %d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &ttype, &tisReal);
01192             TupleSignature oneSig(7,CROSSTERM,(Index)ttype, (char)tisReal);
01193             oneSig.offset[0] = tmp1;
01194             oneSig.offset[1] = tmp2;
01195             oneSig.offset[2] = tmp3;
01196             oneSig.offset[3] = tmp4;
01197             oneSig.offset[4] = tmp5;
01198             oneSig.offset[5] = tmp6;
01199             oneSig.offset[6] = tmp7;
01200             atomSigPool[i].crosstermSigs[j] = oneSig;
01201         }        
01202     }
01203 
01204     //Reading exclusions' signatures    
01205     NAMD_read_line(psf_file, buffer);
01206     if(!NAMD_find_word(buffer, "NEXCLSIGS")){
01207         NAMD_die("UNABLE TO FIND NEXCLSIGS");
01208     }
01209     sscanf(buffer, "%d", &exclSigPoolSize);
01210     if(exclSigPoolSize>0) exclSigPool = new ExclusionSignature[exclSigPoolSize];
01211     vector<int> fullExcls;
01212     vector<int> modExcls;
01213     for(int i=0; i<exclSigPoolSize; i++){
01214         int fullExclCnt = NAMD_read_int(psf_file, buffer);
01215         for(int j=0; j<fullExclCnt; j++)
01216             fullExcls.push_back(NAMD_read_int(psf_file, buffer));
01217         int modExclCnt = NAMD_read_int(psf_file, buffer);
01218         for(int j=0; j<modExclCnt; j++)
01219             modExcls.push_back(NAMD_read_int(psf_file, buffer));
01220 
01221         //The integers in both vectors should be in the increasing the order
01222         //since they have been sorted during the compression of psf file
01223         exclSigPool[i].setOffsets(fullExcls, modExcls);
01224 
01225         fullExcls.clear();
01226         modExcls.clear();
01227     }
01228 
01229     //Now read the cluster information
01230     NAMD_read_line(psf_file, buffer);
01231     if(!NAMD_find_word(buffer, "NCLUSTERS")) {
01232         NAMD_die("UNABLE TO FIND NCLUSTERS");
01233     }
01234     sscanf(buffer, "%d", &numClusters);
01235 
01236     NAMD_read_line(psf_file, buffer);
01237     if(!NAMD_find_word(buffer, "CLUSTERCONTIGUITY")) {
01238         NAMD_die("UNABLE TO FIND CLUSTERCONTIGUITY");
01239     }
01240     sscanf(buffer, "%d", &isClusterContiguous);
01241     if(!isClusterContiguous) {
01242         clusterSize = new int32[numClusters];
01243         memset(clusterSize, 0, sizeof(int32)*numClusters);
01244     }
01245 
01246     //Now begin to read atoms
01247     //This part could be read in the batch mode. All the above information
01248     //can be stored in memory which only occupies several KBs.
01249     NAMD_read_line(psf_file, buffer);
01250     if(!NAMD_find_word(buffer, "NATOM"))
01251         NAMD_die("UNABLE TO FIND NATOM");
01252     sscanf(buffer, "%d", &numAtoms);
01253 
01254     NAMD_read_line(psf_file, buffer);
01255     if(!NAMD_find_word(buffer, "NHYDROGENGROUP"))
01256         NAMD_die("UNABLE TO FIND NHYDROGENGROUP");
01257     sscanf(buffer, "%d", &numHydrogenGroups);
01258 
01259     int isOccupancyValid, isBFactorValid;
01260     NAMD_read_line(psf_file, buffer);
01261     if(!NAMD_find_word(buffer, "OCCUPANCYVALID"))
01262         NAMD_die("UNABLE TO FIND OCCUPANCYVALID");
01263     sscanf(buffer, "%d", &isOccupancyValid);
01264     NAMD_read_line(psf_file, buffer);
01265     if(!NAMD_find_word(buffer, "TEMPFACTORVALID"))
01266         NAMD_die("UNABLE TO FIND TEMPFACTORVALID");
01267     sscanf(buffer, "%d", &isBFactorValid);
01268 
01269     if(numAtoms>0){
01270         atoms = new AtomCstInfo[numAtoms];
01271         atomNames = new AtomNameIdx[numAtoms];
01272         eachAtomMass = new Index[numAtoms];
01273         eachAtomCharge = new Index[numAtoms];
01274         eachAtomSig = new Index[numAtoms];
01275         eachAtomExclSig = new Index[numAtoms];
01276 
01277         clusterSigs = new int[numAtoms];
01278 
01279         int *tmpHgSize = new int[numAtoms];
01280         int *tmpHgGP = new int[numAtoms];
01281         int *tmpHgSV = new int[numAtoms];
01282         hydrogenGroup.resize(numAtoms);
01283         HydrogenGroupID *hg = hydrogenGroup.begin();
01284 
01285         ResidueLookupElem *tmpResLookup = resLookup;
01286 
01287         if(isOccupancyValid) {
01288             occupancy = new float[numAtoms];
01289         }
01290         if(isBFactorValid) {
01291             bfactor = new float[numAtoms];
01292         }
01293 
01294         int residue_number; //for residue number
01295         char *segment_name;
01296         int read_count;        
01297         float tmpf[2];
01298         
01299 #if BINARY_PERATOM_OUTPUT
01300         //1. open the binary per-atom info file (fname.bin)
01301         char *binFName = new char[strlen(fname)+10];
01302         sprintf(binFName, "%s.bin", fname);
01303         FILE *perAtomFile = Fopen(binFName, "rb");
01304         if(perAtomFile==NULL) {
01305             char err_msg[512];
01306             sprintf(err_msg, "UNABLE TO OPEN THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s", binFName);            
01307             NAMD_die(err_msg);
01308         }                
01309         //2. read the magic number to determine whether the endian is same or not
01310         int needFlip = 0;
01311         int magicNum = COMPRESSED_PSF_MAGICNUM;
01312         int rMagicNum = COMPRESSED_PSF_MAGICNUM;
01313         flipNum((char *)&rMagicNum, sizeof(int), 1);
01314         int fMagicNum;
01315         fread(&fMagicNum, sizeof(int), 1, perAtomFile);
01316         if(fMagicNum==magicNum) {
01317             needFlip = 0;
01318         }else if(fMagicNum==rMagicNum){
01319             needFlip = 1;
01320         }else{
01321             char err_msg[512];
01322             sprintf(err_msg, "THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS CORRUPTED", binFName);            
01323             NAMD_die(err_msg);
01324         }
01325         //3. read the file version number
01326         float verNum =  0.0f;
01327         fread(&verNum, sizeof(float), 1, perAtomFile);
01328         if(needFlip) flipNum((char *)&verNum, sizeof(float), 1);
01329         if(fabs(verNum - COMPRESSED_PSF_VER)>1e-6) {
01330             char err_msg[512];
01331             sprintf(err_msg, "THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n", binFName);
01332             NAMD_die(err_msg);
01333         }
01334         delete[] binFName;
01335         //4. read per-atom info
01336         Index sIdx[8];
01337         int iIdx[8];
01338         for(int i=0; i<numAtoms; i++){   
01339             fread(sIdx, sizeof(Index), 8, perAtomFile);
01340             fread(iIdx, sizeof(int), 7, perAtomFile);
01341             fread(tmpf, sizeof(float), 2, perAtomFile);
01342             if(needFlip) {
01343                 flipNum((char *)sIdx, sizeof(Index), 8);
01344                 flipNum((char *)iIdx, sizeof(int), 7);
01345                 flipNum((char *)tmpf, sizeof(float), 2);
01346             }
01347             segment_name = segNamePool[sIdx[0]];                   
01348             atomNames[i].resnameIdx = sIdx[1];
01349             atomNames[i].atomnameIdx = sIdx[2];
01350             atomNames[i].atomtypeIdx = sIdx[3];
01351             eachAtomCharge[i] = sIdx[4];
01352             eachAtomMass[i] = sIdx[5];
01353             eachAtomSig[i] = sIdx[6];
01354             eachAtomExclSig[i] = sIdx[7];
01355             residue_number = iIdx[0];
01356             clusterSigs[i] = iIdx[1];
01357             if(!isClusterContiguous)
01358                 clusterSize[iIdx[1]]++;
01359 
01360             atoms[i].partner = iIdx[2];
01361             atoms[i].hydrogenList= iIdx[3];
01362 
01363             tmpHgSize[i] = iIdx[4];
01364             tmpHgGP[i] = iIdx[5];
01365             tmpHgSV[i] = iIdx[6];
01366 
01367             //debugExclNum += (exclSigPool[sIdx[7]].fullExclCnt+exclSigPool[sIdx[7]].modExclCnt);
01368 #else
01369         int idx[15];
01370         for(int i=0; i<numAtoms; i++){
01371             NAMD_read_line(psf_file, buffer);
01372             read_count = sscanf(buffer, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %f %f",
01373                                 idx, idx+1, idx+2, idx+3, idx+4,
01374                                 idx+5, idx+6, idx+7, idx+8, idx+9, idx+10, idx+11, idx+12,
01375                                 idx+13, idx+14, tmpf, tmpf+1);
01376             if(read_count!=17){
01377                 char err_msg[128];
01378                 sprintf(err_msg, "BAD ATOM LINE FORMAT IN COMPRESSED PSF FILE IN ATOM LINE %d\nLINE=%s",i+1, buffer);
01379                 NAMD_die(err_msg);
01380             }
01381             segment_name = segNamePool[idx[0]];                   
01382             atomNames[i].resnameIdx = (Index)idx[1];
01383             atomNames[i].atomnameIdx = (Index)idx[2];
01384             atomNames[i].atomtypeIdx = (Index)idx[3];
01385             eachAtomCharge[i] = (Index)idx[4];
01386             eachAtomMass[i] = (Index)idx[5];              
01387             eachAtomSig[i] = (Index)idx[6];
01388             eachAtomExclSig[i] = (Index)idx[7];
01389             residue_number = idx[8];
01390             clusterSigs[i] = idx[9];
01391             if(!isClusterContiguous)
01392                 clusterSize[idx[9]]++;
01393 
01394             atoms[i].partner = idx[10];
01395             atoms[i].hydrogenList= idx[11];
01396 
01397             tmpHgSize[i] = idx[12];
01398             tmpHgGP[i] = idx[13];
01399             tmpHgSV[i] = idx[14];
01400             //debugExclNum += (exclSigPool[idx[7]].fullExclCnt+exclSigPool[idx[7]].modExclCnt);
01401 #endif      
01402 
01403             if(isOccupancyValid)
01404                 occupancy[i] = tmpf[0];
01405             if(isBFactorValid)
01406                 bfactor[i] = tmpf[1];
01407 
01408             //Add this atom to residue lookup table
01409             if(tmpResLookup) tmpResLookup =
01410                 tmpResLookup->append(segment_name, residue_number, i);
01411             //Determine the type of the atom (H or O)
01412             Real thisAtomMass = atomMassPool[eachAtomMass[i]];
01413             if (thisAtomMass <= 0.05) {
01414               atoms[i].status |= LonepairAtom;
01415             } else if (thisAtomMass < 1.0) {
01416               atoms[i].status |= DrudeAtom;
01417             } else if(thisAtomMass <= 3.5){
01418                 atoms[i].status = HydrogenAtom;
01419             }else if(atomNamePool[atomNames[i].atomnameIdx][0]=='O' &&
01420                      (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)){
01421                 atoms[i].status = OxygenAtom;
01422             }
01423 
01424             //Look up the vdw constants for this atom
01425             params->assign_vdw_index(atomTypePool[atomNames[i].atomtypeIdx],
01426                                      &(atoms[i]));
01427         } //end of reading per-atom information
01428 
01429         //re-sort the hydrogenGroup array based on the recorded "hydrogenList"
01430         //value
01431         for(int i=0; i<numAtoms; i++){
01432             int hgIdx = atoms[i].hydrogenList;
01433             hg[hgIdx].atomID = i;
01434             hg[hgIdx].atomsInGroup = tmpHgSize[i];
01435             hg[hgIdx].GPID = tmpHgGP[i];
01436             hg[hgIdx].sortVal = tmpHgSV[i];
01437             hg[hgIdx].isGP = 1;
01438             if(tmpHgSize[i]==0) hg[hgIdx].isGP = 0;
01439         }
01440         
01441         delete [] tmpHgSize;
01442         delete [] tmpHgGP;
01443         delete [] tmpHgSV;      
01444 
01445         //printf("debugExclNum: %d\n", debugExclNum);
01446     }
01447     //build up information for numBonds/numDihedrals.. etc
01448     numBonds=0;
01449     numAngles=0;
01450     numDihedrals=0;
01451     numImpropers=0;
01452     numCrossterms=0;
01453     numTotalExclusions=0;
01454     for(int i=0; i<numAtoms; i++){
01455         AtomSignature *thisSig = &atomSigPool[eachAtomSig[i]];
01456         numBonds += thisSig->bondCnt;
01457         numAngles += thisSig->angleCnt;
01458         numDihedrals += thisSig->dihedralCnt;
01459         numImpropers += thisSig->improperCnt;
01460         numCrossterms += thisSig->crosstermCnt;
01461 
01462         ExclusionSignature *exclSig = &exclSigPool[eachAtomExclSig[i]];
01463         numTotalExclusions += (exclSig->fullExclCnt + exclSig->modExclCnt);
01464     }
01465     
01466     numTotalExclusions /= 2;
01467 
01468     //Just reading for the parameters values; extra Bonds, Dihedrals etc.
01469     //have been taken into account when compressing the molecule object.
01470     //The actual number of Bonds, Dihedrals etc. will be calculated based
01471     //on atom signatures.
01472     if(simParams->extraBondsOn)
01473         build_extra_bonds(params, cfgList->find("extraBondsFile"));
01474 #if 0
01475     //This part has been enabled in build_extra_bonds for memory optimized version
01476     //read extra bond parameters if there is an input of extra bonds (extraBondsOn is true)
01477     int extraDihedralParamNum = 0;
01478     int extraImproperParamNum = 0;
01479     if(simParams->extraBondsOn){
01480         int numExtraParams=0;
01481 
01482         //read extra bond params
01483         NAMD_read_line(psf_file, buffer);
01484         if(!NAMD_find_word(buffer, "NEXTRABONDPARAMS")){
01485             NAMD_die("UNABLE TO FIND NEXTRABONDPARAMS");
01486         }
01487         sscanf(buffer, "%d", &numExtraParams);
01488         if(numExtraParams>0){
01489             BondValue *newParams = new BondValue[params->NumBondParams+numExtraParams];
01490             memcpy(newParams, params->bond_array, params->NumBondParams*sizeof(BondValue));
01491             delete [] params->bond_array;            
01492 
01493             int curNumPs = params->NumBondParams;
01494             for(int i=0; i<numExtraParams; i++){
01495                 Real k, x0;
01496                 NAMD_read_line(psf_file, buffer);
01497                 sscanf(buffer, "%f %f", &k, &x0);
01498                 newParams[curNumPs+i].k = k;
01499                 newParams[curNumPs+i].x0 = x0;
01500             }
01501             params->bond_array = newParams;
01502             params->NumBondParams += numExtraParams;
01503         }
01504 
01505         //read extra angle params
01506         NAMD_read_line(psf_file, buffer);
01507         if(!NAMD_find_word(buffer, "NEXTRAANGLEPARAMS")){
01508             NAMD_die("UNABLE TO FIND NEXTRAANGLEPARAMS");
01509         }
01510         sscanf(buffer, "%d", &numExtraParams);
01511         if(numExtraParams>0){
01512             AngleValue *newParams = new AngleValue[params->NumAngleParams+numExtraParams];
01513             memcpy(newParams, params->angle_array, params->NumAngleParams*sizeof(AngleValue));
01514             delete [] params->angle_array;
01515 
01516             int curNumPs = params->NumAngleParams;
01517             for(int i=0; i<numExtraParams; i++) {
01518                 Real k, theta0, k_ub, r_ub;
01519                 NAMD_read_line(psf_file, buffer);
01520                 sscanf(buffer, "%f %f %f %f", &k, &theta0, &k_ub, &r_ub);
01521                 newParams[curNumPs+i].k = k;
01522                 newParams[curNumPs+i].theta0 = theta0;
01523                 newParams[curNumPs+i].k_ub = k_ub;
01524                 newParams[curNumPs+i].r_ub = r_ub;
01525             }
01526             params->angle_array = newParams;
01527             params->NumAngleParams += numExtraParams;
01528         }
01529         
01530         //read extra diheral params
01531         NAMD_read_line(psf_file, buffer);
01532         if(!NAMD_find_word(buffer, "NEXTRADIHEDRALPARAMS")){
01533             NAMD_die("UNABLE TO FIND NEXTRADIHEDRALPARAMS");
01534         }
01535         sscanf(buffer, "%d", &numExtraParams);
01536         extraDihedralParamNum = numExtraParams;
01537         if(numExtraParams>0){
01538             DihedralValue *newParams = new DihedralValue[params->NumDihedralParams+numExtraParams];
01539             memcpy(newParams, params->dihedral_array, params->NumDihedralParams*sizeof(DihedralValue));
01540             delete [] params->dihedral_array;
01541 
01542             int curNumPs = params->NumDihedralParams;
01543             for(int i=0; i<numExtraParams; i++) {
01544                 int multiplicity;
01545                 Real k, delta;
01546                 int n;
01547                 NAMD_read_line(psf_file, buffer); //read multiplicity
01548                 sscanf(buffer, "%d", &multiplicity);
01549                 newParams[curNumPs+i].multiplicity = multiplicity;
01550                 for(int j=0; j<multiplicity; j++) {
01551                     NAMD_read_line(psf_file, buffer);
01552                     sscanf(buffer, "%f %d %f", &k, &n, &delta);
01553                     newParams[curNumPs+i].values[j].k = k;
01554                     newParams[curNumPs+i].values[j].n = n;
01555                     newParams[curNumPs+i].values[j].delta = delta;
01556                 }
01557             }
01558             params->dihedral_array = newParams;
01559         }
01560 
01561         //read extra improper params
01562         NAMD_read_line(psf_file, buffer);
01563         if(!NAMD_find_word(buffer, "NEXTRAIMPROPERPARAMS")){
01564             NAMD_die("UNABLE TO FIND NEXTRAIMPROPERPARAMS");
01565         }
01566         sscanf(buffer, "%d", &numExtraParams);
01567         extraImproperParamNum = numExtraParams;
01568         if(numExtraParams>0){
01569             ImproperValue *newParams = new ImproperValue[params->NumImproperParams+numExtraParams];
01570             memcpy(newParams, params->improper_array, params->NumImproperParams*sizeof(ImproperValue));
01571             delete [] params->improper_array;
01572 
01573             int curNumPs = params->NumImproperParams;
01574             for(int i=0; i<numExtraParams; i++) {
01575                 int multiplicity;
01576                 Real k, delta;
01577                 int n;
01578                 NAMD_read_line(psf_file, buffer); //read multiplicity
01579                 sscanf(buffer, "%d", &multiplicity);
01580                 newParams[curNumPs+i].multiplicity = multiplicity;
01581                 for(int j=0; j<multiplicity; j++) {
01582                     NAMD_read_line(psf_file, buffer);
01583                     sscanf(buffer, "%f %d %f", &k, &n, &delta);
01584                     newParams[curNumPs+i].values[j].k = k;
01585                     newParams[curNumPs+i].values[j].n = n;
01586                     newParams[curNumPs+i].values[j].delta = delta;
01587                 }
01588             }
01589             params->improper_array = newParams;
01590         }
01591     }
01592 #endif
01593 
01594     //read DIHEDRALPARAMARRAY and IMPROPERPARAMARRAY    
01595     NAMD_read_line(psf_file, buffer);
01596     if(!NAMD_find_word(buffer, "DIHEDRALPARAMARRAY"))
01597         NAMD_die("UNABLE TO FIND DIHEDRALPARAMARRAY");
01598     for(int i=0; i<params->NumDihedralParams; i++){
01599         params->dihedral_array[i].multiplicity = NAMD_read_int(psf_file, buffer);
01600     }
01601     //params->NumDihedralParams += extraDihedralParamNum;
01602 
01603     NAMD_read_line(psf_file, buffer); //to read a simple single '\n' line 
01604     NAMD_read_line(psf_file, buffer);
01605     if(!NAMD_find_word(buffer, "IMPROPERPARAMARRAY"))
01606         NAMD_die("UNABLE TO FIND IMPROPERPARAMARRAY");
01607     for(int i=0; i<params->NumImproperParams; i++){
01608         params->improper_array[i].multiplicity = NAMD_read_int(psf_file, buffer);
01609     }
01610     //params->NumImproperParams += extraImproperParamNum;
01611     Fclose(psf_file);
01612 
01613     //numRealBonds is set when reading bondSignatures from compressed psf file
01614     
01615     build_atom_status();
01616 #endif
01617 }
01618 /*      END OF FUNCTION read_compressed_psf_file      */
01619 
01620 /************************************************************************/
01621 /*                  */
01622 /*        FUNCTION read_atoms      */
01623 /*                  */
01624 /*   INPUTS:                */
01625 /*  fd - file pointer to the .psf file        */
01626 /*  params - Parameters object to use for parameters    */
01627 /*                  */
01628 /*  this function reads in the Atoms section of the .psf file.      */
01629 /*   This section consists of numAtoms lines that are of the form:  */
01630 /*     <atom#> <mol> <seg#> <res> <atomname> <atomtype> <charge> <mass> */
01631 /*   Each line is read into the appropriate entry in the atoms array.   */
01632 /*   The parameters object is then used to determine the vdW constants  */
01633 /*   for this atom.              */
01634 /*                  */
01635 /************************************************************************/
01636 
01637 void Molecule::read_atoms(FILE *fd, Parameters *params)
01638 
01639 {
01640 #ifdef MEM_OPT_VERSION
01641     return;
01642 #else
01643   char buffer[512];  // Buffer for reading from file
01644   int atom_number=0;  // Atom number 
01645   int last_atom_number=0; // Last atom number, used to assure
01646         // atoms are in order
01647   char segment_name[11]; // Segment name
01648   char residue_number[11]; // Residue number
01649   char residue_name[11];  // Residue name
01650   char atom_name[11];  // Atom name
01651   char atom_type[11];  // Atom type
01652   Real charge;    // Charge for the current atom
01653   Real mass;    // Mass for the current atom
01654   int read_count;    // Number of fields read by sscanf
01655 
01656   /*  Allocate the atom arrays          */
01657   atoms     = new Atom[numAtoms];
01658   atomNames = new AtomNameInfo[numAtoms];
01659   if(simParams->genCompressedPsf) {
01660       atomSegResids = new AtomSegResInfo[numAtoms];
01661   }
01662 
01663   // DRUDE: supplement Atom data
01664   if (is_drude_psf) {
01665     drudeConsts = new DrudeConst[numAtoms];
01666   }
01667   // DRUDE
01668 
01669   hydrogenGroup.resize(0);
01670 
01671   if (atoms == NULL || atomNames == NULL )
01672   {
01673     NAMD_die("memory allocation failed in Molecule::read_atoms");
01674   }
01675 
01676   ResidueLookupElem *tmpResLookup = resLookup;
01677 
01678   /*  Loop and read in numAtoms atom lines.      */
01679   while (atom_number < numAtoms)
01680   {
01681     /*  Get the line from the file        */
01682     NAMD_read_line(fd, buffer);
01683 
01684     /*  If its blank or a comment, skip it      */
01685     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') )
01686       continue;
01687 
01688     /*  Parse up the line          */
01689     read_count=sscanf(buffer, "%d %s %s %s %s %s %f %f",
01690        &atom_number, segment_name, residue_number,
01691        residue_name, atom_name, atom_type, &charge, &mass);
01692 
01693     /*  Check to make sure we found what we were expecting  */
01694     if (read_count != 8)
01695     {
01696       char err_msg[128];
01697 
01698       sprintf(err_msg, "BAD ATOM LINE FORMAT IN PSF FILE IN ATOM LINE %d\nLINE=%s",
01699          last_atom_number+1, buffer);
01700       NAMD_die(err_msg);
01701     }
01702 
01703     // DRUDE: read alpha and thole parameters from atom line
01704     if (is_drude_psf)
01705     {
01706       Real alpha, thole;
01707       read_count=sscanf(buffer,
01708 //          "%*d %*s %*s %*s %*s %*s %*f %*f %*d %*f %*f %f %f", &alpha, &thole);
01709                 // the two columns preceding alpha and thole will disappear
01710           "%*d %*s %*s %*s %*s %*s %*f %*f %*d %f %f", &alpha, &thole);
01711       if (read_count != 2)
01712       {
01713         char err_msg[128];
01714 
01715         sprintf(err_msg, "BAD ATOM LINE FORMAT IN PSF FILE "
01716             "IN ATOM LINE %d\nLINE=%s", last_atom_number+1, buffer);
01717         NAMD_die(err_msg);
01718       }
01719       drudeConsts[atom_number-1].alpha = alpha;
01720       drudeConsts[atom_number-1].thole = thole;
01721     }
01722     // DRUDE
01723 
01724     /*  Check if this is in XPLOR format  */
01725     int atom_type_num;
01726     if ( sscanf(atom_type, "%d", &atom_type_num) > 0 )
01727     {
01728       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.");
01729     }
01730 
01731     /*  Make sure the atoms were in sequence    */
01732     if (atom_number != last_atom_number+1)
01733     {
01734       char err_msg[128];
01735 
01736       sprintf(err_msg, "ATOM NUMBERS OUT OF ORDER AT ATOM #%d OF PSF FILE",
01737          last_atom_number+1);
01738       NAMD_die(err_msg);
01739     }
01740 
01741     last_atom_number++;
01742 
01743     /*  Dynamically allocate strings for atom name, atom    */
01744     /*  type, etc so that we only allocate as much space    */
01745     /*  for these strings as we really need      */
01746     int reslength = strlen(residue_name)+1;
01747     int namelength = strlen(atom_name)+1;
01748     int typelength = strlen(atom_type)+1;
01749 
01750     atomNames[atom_number-1].resname = nameArena->getNewArray(reslength);
01751     atomNames[atom_number-1].atomname = nameArena->getNewArray(namelength);
01752     atomNames[atom_number-1].atomtype = nameArena->getNewArray(typelength);
01753   
01754     if (atomNames[atom_number-1].resname == NULL)
01755     {
01756       NAMD_die("memory allocation failed in Molecule::read_atoms");
01757     }
01758 
01759     /*  Put the values from this atom into the atoms array  */
01760     strcpy(atomNames[atom_number-1].resname, residue_name);
01761     strcpy(atomNames[atom_number-1].atomname, atom_name);
01762     strcpy(atomNames[atom_number-1].atomtype, atom_type);
01763     atoms[atom_number-1].mass = mass;
01764     atoms[atom_number-1].charge = charge;
01765     atoms[atom_number-1].status = UnknownAtom;
01766 
01767     /*  Add this atom to residue lookup table */
01768     if ( tmpResLookup ) tmpResLookup =
01769         tmpResLookup->append(segment_name, atoi(residue_number), atom_number-1);
01770 
01771     if(atomSegResids) { //for compressing molecule information
01772         AtomSegResInfo *one = atomSegResids + (atom_number - 1);
01773         memcpy(one->segname, segment_name, strlen(segment_name)+1);
01774         one->resid = atoi(residue_number);
01775     }
01776 
01777     /*  Determine the type of the atom (H or O) */
01778     if (atoms[atom_number-1].mass <= 0.05) {
01779       atoms[atom_number-1].status |= LonepairAtom;
01780     } else if (atoms[atom_number-1].mass < 1.0) {
01781       atoms[atom_number-1].status |= DrudeAtom;
01782     } else if (atoms[atom_number-1].mass <=3.5) {
01783       atoms[atom_number-1].status |= HydrogenAtom;
01784     } else if ((atomNames[atom_number-1].atomname[0] == 'O') && 
01785          (atoms[atom_number-1].mass >= 14.0) && 
01786          (atoms[atom_number-1].mass <= 18.0)) {
01787       atoms[atom_number-1].status |= OxygenAtom;
01788     }
01789 
01790     /*  Look up the vdw constants for this atom    */
01791     params->assign_vdw_index(atomNames[atom_number-1].atomtype, 
01792        &(atoms[atom_number-1]));
01793         }
01794 
01795   return;
01796 #endif
01797 }
01798 /*      END OF FUNCTION read_atoms      */
01799 
01800 /************************************************************************/
01801 /*                  */
01802 /*      FUNCTION read_bonds        */
01803 /*                  */
01804 /*  read_bonds reads in the bond section of the .psf file.  This    */
01805 /*  section contains a list of pairs of numbers where each pair is      */
01806 /*  represents two atoms that are bonded together.  Each atom pair is   */
01807 /*  read in.  Then that parameter object is queried to determine the    */
01808 /*  force constant and rest distance for the bond.      */
01809 /*                  */
01810 /************************************************************************/
01811 
01812 void Molecule::read_bonds(FILE *fd, Parameters *params)
01813 
01814 {
01815 #ifdef MEM_OPT_VERSION
01816     return;
01817 #else
01818   int atom_nums[2];  // Atom indexes for the bonded atoms
01819   char atom1name[11];  // Atom type for atom #1
01820   char atom2name[11];  // Atom type for atom #2
01821   register int j;      // Loop counter
01822   int num_read=0;    // Number of bonds read so far
01823   int origNumBonds = numBonds;   // number of bonds in file header
01824 
01825   /*  Allocate the array to hold the bonds      */
01826   bonds=new Bond[numBonds];
01827 
01828   if (bonds == NULL)
01829   {
01830     NAMD_die("memory allocations failed in Molecule::read_bonds");
01831   }
01832 
01833   /*  Loop through and read in all the bonds      */
01834   while (num_read < numBonds)
01835   {
01836     /*  Loop and read in the two atom indexes    */
01837     for (j=0; j<2; j++)
01838     {
01839       /*  Read the atom number from the file.         */
01840       /*  Subtract 1 to convert the index from the    */
01841       /*  1 to NumAtoms used in the file to the       */
01842       /*  0 to NumAtoms-1 that we need    */
01843       atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
01844 
01845       /*  Check to make sure the index isn't too big  */
01846       if (atom_nums[j] >= numAtoms)
01847       {
01848         char err_msg[128];
01849 
01850         sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01851         NAMD_die(err_msg);
01852       }
01853     }
01854 
01855     /*  Get the atom type for the two atoms.  When we query */
01856     /*  the parameter object, we need to send the atom type */
01857     /*  that is alphabetically first as atom 1.    */
01858     if (strcasecmp(atomNames[atom_nums[0]].atomtype, 
01859          atomNames[atom_nums[1]].atomtype) < 0)
01860     {
01861       strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
01862       strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
01863     }
01864     else
01865     {
01866       strcpy(atom2name, atomNames[atom_nums[0]].atomtype);
01867       strcpy(atom1name, atomNames[atom_nums[1]].atomtype);
01868     }
01869 
01870     /*  Query the parameter object for the constants for    */
01871     /*  this bond            */
01872     Bond *b = &(bonds[num_read]);
01873     params->assign_bond_index(atom1name, atom2name, b);
01874 
01875     /*  Assign the atom indexes to the array element  */
01876     b->atom1=atom_nums[0];
01877     b->atom2=atom_nums[1];
01878 
01879     /*  Make sure this isn't a fake bond meant for shake in x-plor.  */
01880     Real k, x0;
01881     params->get_bond_params(&k,&x0,b->bond_type);
01882     if (simParams->watmodel == WAT_SWM4) {
01883       // need to retain Lonepair bonds for Drude
01884       if ( k == 0. && !is_lp(b->atom1) && !is_lp(b->atom2)) --numBonds;
01885       else ++num_read;
01886     }
01887     else {
01888       if ( k == 0. ) --numBonds;  // fake bond
01889       else ++num_read;  // real bond
01890     }
01891   }
01892 
01893   /*  Tell user about our subterfuge  */
01894   if ( numBonds != origNumBonds ) {
01895     iout << iWARN << "Ignored " << origNumBonds - numBonds <<
01896             " bonds with zero force constants.\n" << endi;
01897     iout << iWARN <<
01898         "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
01899   }
01900 
01901   return;
01902 #endif
01903 }
01904 /*      END OF FUNCTION read_bonds      */
01905 
01906 /************************************************************************/
01907 /*                  */
01908 /*      FUNCTION read_angles        */
01909 /*                  */
01910 /*   INPUTS:                */
01911 /*  fd - File descriptor for .psf file        */
01912 /*  params - Parameters object to query for parameters    */
01913 /*                  */
01914 /*  read_angles reads the angle parameters from the .psf file.      */
01915 /*   This section of the .psf file consists of a list of triplets of    */
01916 /*   atom indexes.  Each triplet represents three atoms connected via   */
01917 /*   an angle bond.  The parameter object is queried to obtain the      */
01918 /*   constants for each bond.            */
01919 /*                  */
01920 /************************************************************************/
01921 
01922 void Molecule::read_angles(FILE *fd, Parameters *params)
01923 
01924 {
01925 #ifdef MEM_OPT_VERSION
01926     return;
01927 #else
01928   int atom_nums[3];  //  Atom numbers for the three atoms
01929   char atom1name[11];  //  Atom type for atom 1
01930   char atom2name[11];  //  Atom type for atom 2
01931   char atom3name[11];  //  Atom type for atom 3
01932   register int j;      //  Loop counter
01933   int num_read=0;    //  Number of angles read so far
01934   int origNumAngles = numAngles;  // Number of angles in file
01935   /*  Alloc the array of angles          */
01936   angles=new Angle[numAngles];
01937 
01938   if (angles == NULL)
01939   {
01940     NAMD_die("memory allocation failed in Molecule::read_angles");
01941   }
01942 
01943   /*  Loop through and read all the angles      */
01944   while (num_read < numAngles)
01945   {
01946     /*  Loop through the 3 atom indexes in the current angle*/
01947     for (j=0; j<3; j++)
01948     {
01949       /*  Read the atom number from the file.         */
01950       /*  Subtract 1 to convert the index from the    */
01951       /*  1 to NumAtoms used in the file to the       */
01952       /*  0 to NumAtoms-1 that we need    */
01953       atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
01954 
01955       /*  Check to make sure the atom index doesn't   */
01956       /*  exceed the Number of Atoms      */
01957       if (atom_nums[j] >= numAtoms)
01958       {
01959         char err_msg[128];
01960 
01961         sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01962         NAMD_die(err_msg);
01963       }
01964     }
01965 
01966     /*  Place the bond name that is alphabetically first  */
01967     /*  in the atom1name.  This is OK since the order of    */
01968     /*  atom1 and atom3 are interchangable.  And to search  */
01969     /*  the tree of angle parameters, we need the order     */
01970     /*  to be predictable.          */
01971     if (strcasecmp(atomNames[atom_nums[0]].atomtype, 
01972          atomNames[atom_nums[2]].atomtype) < 0)
01973     {
01974       strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
01975       strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
01976       strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
01977     }
01978     else
01979     {
01980       strcpy(atom1name, atomNames[atom_nums[2]].atomtype);
01981       strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
01982       strcpy(atom3name, atomNames[atom_nums[0]].atomtype);
01983     }
01984 
01985     /*  Get the constant values for this bond from the  */
01986     /*  parameter object          */
01987     params->assign_angle_index(atom1name, atom2name, 
01988        atom3name, &(angles[num_read]));
01989 
01990     /*  Assign the three atom indices      */
01991     angles[num_read].atom1=atom_nums[0];
01992     angles[num_read].atom2=atom_nums[1];
01993     angles[num_read].atom3=atom_nums[2];
01994 
01995     /*  Make sure this isn't a fake angle meant for shake in x-plor.  */
01996     Real k, t0, k_ub, r_ub;
01997     params->get_angle_params(&k,&t0,&k_ub,&r_ub,angles[num_read].angle_type);
01998     if ( k == 0. && k_ub == 0. ) --numAngles;  // fake angle
01999     else ++num_read;  // real angle
02000   }
02001 
02002   /*  Tell user about our subterfuge  */
02003   if ( numAngles != origNumAngles ) {
02004     iout << iWARN << "Ignored " << origNumAngles - numAngles <<
02005             " angles with zero force constants.\n" << endi;
02006   }
02007 
02008   return;
02009 #endif
02010 }
02011 /*      END OF FUNCTION read_angles      */
02012 
02013 /************************************************************************/
02014 /*                  */
02015 /*        FUNCTION read_dihedrals      */
02016 /*                  */
02017 /*   INPUTS:                */
02018 /*  fd - file descriptor for the .psf file        */
02019 /*  params - pointer to parameter object        */
02020 /*                  */
02021 /*  read_dihedreals reads the dihedral section of the .psf file.    */
02022 /*   This section of the file contains a list of quartets of atom       */
02023 /*   numbers.  Each quartet represents a group of atoms that form a     */
02024 /*   dihedral bond.              */
02025 /*                  */
02026 /************************************************************************/
02027 
02028 void Molecule::read_dihedrals(FILE *fd, Parameters *params)
02029 {
02030 #ifdef MEM_OPT_VERSION
02031     return;
02032 #else
02033   int atom_nums[4];  // The 4 atom indexes
02034   int last_atom_nums[4];  // Atom numbers from previous bond
02035   char atom1name[11];  // Atom type for atom 1
02036   char atom2name[11];  // Atom type for atom 2
02037   char atom3name[11];  // Atom type for atom 3
02038   char atom4name[11];  // Atom type for atom 4
02039   register int j;      // loop counter
02040   int num_read=0;    // number of dihedrals read so far
02041   int multiplicity=1;  // multiplicity of the current bond
02042   Bool duplicate_bond;  // Is this a duplicate of the last bond
02043   int num_unique=0;   // Number of unique dihedral bonds
02044 
02045   //  Initialize the array used to check for duplicate dihedrals
02046   for (j=0; j<4; j++)
02047     last_atom_nums[j] = -1;
02048 
02049   /*  Allocate an array to hold the Dihedrals      */
02050   dihedrals = new Dihedral[numDihedrals];
02051 
02052   if (dihedrals == NULL)
02053   {
02054     NAMD_die("memory allocation failed in Molecule::read_dihedrals");
02055   }
02056 
02057   /*  Loop through and read all the dihedrals      */
02058   while (num_read < numDihedrals)
02059   {
02060     duplicate_bond = TRUE;
02061 
02062     /*  Loop through and read the 4 indexes for this bond   */
02063     for (j=0; j<4; j++)
02064     {
02065       /*  Read the atom number from the file.         */
02066       /*  Subtract 1 to convert the index from the    */
02067       /*  1 to NumAtoms used in the file to the       */
02068       /*  0 to NumAtoms-1 that we need    */
02069       atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
02070 
02071       /*  Check for an atom index that is too large  */
02072       if (atom_nums[j] >= numAtoms)
02073       {
02074         char err_msg[128];
02075 
02076         sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
02077         NAMD_die(err_msg);
02078       }
02079 
02080       //  Check to see if this atom matches the last bond
02081       if (atom_nums[j] != last_atom_nums[j])
02082       {
02083         duplicate_bond = FALSE;
02084       }
02085 
02086       last_atom_nums[j] = atom_nums[j];
02087     }
02088 
02089     /*  Get the atom types for the 4 atoms so we can look  */
02090     /*  up the constants in the parameter object    */
02091     strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
02092     strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
02093     strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
02094     strcpy(atom4name, atomNames[atom_nums[3]].atomtype);
02095 
02096     //  Check to see if this is really a new bond or just
02097     //  a repeat of the last one
02098     if (duplicate_bond)
02099     {
02100       //  This is a duplicate, so increase the multiplicity
02101       multiplicity++;
02102 
02103       if (multiplicity == 2)
02104       {
02105         numMultipleDihedrals++;
02106       }
02107     }
02108     else
02109     {
02110       multiplicity=1;
02111       num_unique++;
02112     }
02113 
02114     /*  Get the constants for this dihedral bond    */
02115     params->assign_dihedral_index(atom1name, atom2name, 
02116        atom3name, atom4name, &(dihedrals[num_unique-1]),
02117        multiplicity);
02118 
02119     /*  Assign the atom indexes        */
02120     dihedrals[num_unique-1].atom1=atom_nums[0];
02121     dihedrals[num_unique-1].atom2=atom_nums[1];
02122     dihedrals[num_unique-1].atom3=atom_nums[2];
02123     dihedrals[num_unique-1].atom4=atom_nums[3];
02124 
02125     num_read++;
02126   }
02127 
02128   numDihedrals = num_unique;
02129 
02130   return;
02131 #endif
02132 }
02133 /*      END OF FUNCTION read_dihedral      */
02134 
02135 /************************************************************************/
02136 /*                  */
02137 /*        FUNCTION read_impropers      */
02138 /*                  */
02139 /*   INPUTS:                */
02140 /*  fd - file descriptor for .psf file        */
02141 /*  params - parameter object          */
02142 /*                  */
02143 /*  read_impropers reads the improper section of the .psf file.  */
02144 /*   This section is identical to the dihedral section in that it is    */
02145 /*   made up of a list of quartets of atom indexes that define the      */
02146 /*   atoms that are bonded together.          */
02147 /*                  */
02148 /************************************************************************/
02149 
02150 void Molecule::read_impropers(FILE *fd, Parameters *params)
02151 {
02152 #ifdef MEM_OPT_VERSION
02153     return;
02154 #else
02155   int atom_nums[4];  //  Atom indexes for the 4 atoms
02156   int last_atom_nums[4];  //  Atom indexes from previous bond
02157   char atom1name[11];  //  Atom type for atom 1
02158   char atom2name[11];  //  Atom type for atom 2
02159   char atom3name[11];  //  Atom type for atom 3
02160   char atom4name[11];  //  Atom type for atom 4
02161   register int j;      //  Loop counter
02162   int num_read=0;    //  Number of impropers read so far
02163   int multiplicity=1;  // multiplicity of the current bond
02164   Bool duplicate_bond;  // Is this a duplicate of the last bond
02165   int num_unique=0;   // Number of unique dihedral bonds
02166 
02167   //  Initialize the array used to look for duplicate improper
02168   //  entries.  Set them all to -1 so we know nothing will match
02169   for (j=0; j<4; j++)
02170     last_atom_nums[j] = -1;
02171 
02172   /*  Allocate the array to hold the impropers      */
02173   impropers=new Improper[numImpropers];
02174 
02175   if (impropers == NULL)
02176   {
02177     NAMD_die("memory allocation failed in Molecule::read_impropers");
02178   }
02179 
02180   /*  Loop through and read all the impropers      */
02181   while (num_read < numImpropers)
02182   {
02183     duplicate_bond = TRUE;
02184 
02185     /*  Loop through the 4 indexes for this improper  */
02186     for (j=0; j<4; j++)
02187     {
02188       /*  Read the atom number from the file.         */
02189       /*  Subtract 1 to convert the index from the    */
02190       /*  1 to NumAtoms used in the file to the       */
02191       /*  0 to NumAtoms-1 that we need    */
02192       atom_nums[j]=NAMD_read_int(fd, "IMPROPERS")-1;
02193 
02194       /*  Check to make sure the index isn't too big  */
02195       if (atom_nums[j] >= numAtoms)
02196       {
02197         char err_msg[128];
02198 
02199         sprintf(err_msg, "IMPROPERS INDEX %d GREATER THAN NATOM %d IN IMPROPERS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
02200         NAMD_die(err_msg);
02201       }
02202 
02203       if (atom_nums[j] != last_atom_nums[j])
02204       {
02205         duplicate_bond = FALSE;
02206       }
02207 
02208       last_atom_nums[j] = atom_nums[j];
02209     }
02210 
02211     /*  Get the atom types so we can look up the parameters */
02212     strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
02213     strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
02214     strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
02215     strcpy(atom4name, atomNames[atom_nums[3]].atomtype);
02216 
02217     //  Check to see if this is a duplicate improper
02218     if (duplicate_bond)
02219     {
02220       //  This is a duplicate improper.  So we don't
02221       //  really count this entry, we just update
02222       //  the parameters object
02223       multiplicity++;
02224 
02225       if (multiplicity == 2)
02226       {
02227         //  Count the number of multiples.
02228         numMultipleImpropers++;
02229       }
02230     }
02231     else
02232     {
02233       //  Not a duplicate
02234       multiplicity = 1;
02235       num_unique++;
02236     }
02237 
02238     /*  Look up the constants for this bond      */
02239     params->assign_improper_index(atom1name, atom2name, 
02240        atom3name, atom4name, &(impropers[num_unique-1]),
02241        multiplicity);
02242 
02243     /*  Assign the atom indexes        */
02244     impropers[num_unique-1].atom1=atom_nums[0];
02245     impropers[num_unique-1].atom2=atom_nums[1];
02246     impropers[num_unique-1].atom3=atom_nums[2];
02247     impropers[num_unique-1].atom4=atom_nums[3];
02248 
02249     num_read++;
02250   }
02251 
02252   //  Now reset the numImpropers value to the number of UNIQUE
02253   //  impropers.  Sure, we waste a few entries in the improper_array
02254   //  on the master node, but it is very little space . . .
02255   numImpropers = num_unique;
02256 
02257   return;
02258 #endif
02259 }
02260 /*      END OF FUNCTION read_impropers      */
02261 
02262 /************************************************************************/
02263 /*                  */
02264 /*        FUNCTION read_crossterms      */
02265 /*                  */
02266 /*   INPUTS:                */
02267 /*  fd - file descriptor for .psf file        */
02268 /*  params - parameter object          */
02269 /*                  */
02270 /*   This section is identical to the dihedral section in that it is    */
02271 /*   made up of a list of quartets of atom indexes that define the      */
02272 /*   atoms that are bonded together.          */
02273 /*                  */
02274 /************************************************************************/
02275 
02276 void Molecule::read_crossterms(FILE *fd, Parameters *params)
02277 
02278 {
02279 #ifdef MEM_OPT_VERSION
02280     return;
02281 #else
02282   int atom_nums[8];  //  Atom indexes for the 4 atoms
02283   int last_atom_nums[8];  //  Atom indexes from previous bond
02284   char atom1name[11];  //  Atom type for atom 1
02285   char atom2name[11];  //  Atom type for atom 2
02286   char atom3name[11];  //  Atom type for atom 3
02287   char atom4name[11];  //  Atom type for atom 4
02288   char atom5name[11];  //  Atom type for atom 5
02289   char atom6name[11];  //  Atom type for atom 6
02290   char atom7name[11];  //  Atom type for atom 7
02291   char atom8name[11];  //  Atom type for atom 8
02292   register int j;      //  Loop counter
02293   int num_read=0;    //  Number of items read so far
02294   Bool duplicate_bond;  // Is this a duplicate of the last bond
02295 
02296   //  Initialize the array used to look for duplicate crossterm
02297   //  entries.  Set them all to -1 so we know nothing will match
02298   for (j=0; j<8; j++)
02299     last_atom_nums[j] = -1;
02300 
02301   /*  Allocate the array to hold the cross-terms */
02302   crossterms=new Crossterm[numCrossterms];
02303 
02304   if (crossterms == NULL)
02305   {
02306     NAMD_die("memory allocation failed in Molecule::read_crossterms");
02307   }
02308 
02309   /*  Loop through and read all the cross-terms      */
02310   while (num_read < numCrossterms)
02311   {
02312     duplicate_bond = TRUE;
02313 
02314     /*  Loop through the 8 indexes for this cross-term */
02315     for (j=0; j<8; j++)
02316     {
02317       /*  Read the atom number from the file.         */
02318       /*  Subtract 1 to convert the index from the    */
02319       /*  1 to NumAtoms used in the file to the       */
02320       /*  0 to NumAtoms-1 that we need    */
02321       atom_nums[j]=NAMD_read_int(fd, "CROSS-TERMS")-1;
02322 
02323       /*  Check to make sure the index isn't too big  */
02324       if (atom_nums[j] >= numAtoms)
02325       {
02326         char err_msg[128];
02327 
02328         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);
02329         NAMD_die(err_msg);
02330       }
02331 
02332       if (atom_nums[j] != last_atom_nums[j])
02333       {
02334         duplicate_bond = FALSE;
02335       }
02336 
02337       last_atom_nums[j] = atom_nums[j];
02338     }
02339 
02340     /*  Get the atom types so we can look up the parameters */
02341     strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
02342     strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
02343     strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
02344     strcpy(atom4name, atomNames[atom_nums[3]].atomtype);
02345     strcpy(atom5name, atomNames[atom_nums[4]].atomtype);
02346     strcpy(atom6name, atomNames[atom_nums[5]].atomtype);
02347     strcpy(atom7name, atomNames[atom_nums[6]].atomtype);
02348     strcpy(atom8name, atomNames[atom_nums[7]].atomtype);
02349 
02350     //  Check to see if this is a duplicate term
02351     if (duplicate_bond)
02352     {
02353       iout << iWARN << "Duplicate cross-term detected.\n" << endi;
02354     }
02355 
02356     /*  Look up the constants for this bond      */
02357     params->assign_crossterm_index(atom1name, atom2name, 
02358        atom3name, atom4name, atom5name, atom6name,
02359        atom7name, atom8name, &(crossterms[num_read]));
02360 
02361     /*  Assign the atom indexes        */
02362     crossterms[num_read].atom1=atom_nums[0];
02363     crossterms[num_read].atom2=atom_nums[1];
02364     crossterms[num_read].atom3=atom_nums[2];
02365     crossterms[num_read].atom4=atom_nums[3];
02366     crossterms[num_read].atom5=atom_nums[4];
02367     crossterms[num_read].atom6=atom_nums[5];
02368     crossterms[num_read].atom7=atom_nums[6];
02369     crossterms[num_read].atom8=atom_nums[7];
02370 
02371     if(!duplicate_bond) num_read++;
02372   }
02373 
02374   numCrossterms = num_read;
02375 
02376   return;
02377 #endif
02378 }
02379 /*      END OF FUNCTION read_impropers      */
02380 
02381 /************************************************************************/
02382 /*                  */
02383 /*      FUNCTION read_donors        */
02384 /*                  */
02385 /*  read_donors reads in the bond section of the .psf file.  This   */
02386 /*  section contains a list of pairs of numbers where each pair is      */
02387 /*  represents two atoms that are part of an H-bond.  Each atom pair is */
02388 /*  read in.                                                            */
02389 /*                  */
02390 /*  Donor atoms are the heavy atoms to which hydrogens are bonded.      */
02391 /*  There will always be a donor atom for each donor pair.  However,    */
02392 /*  for a united-atom model there may not be an explicit hydrogen       */
02393 /*  present, in which case the second atom index in the pair will be    */
02394 /*  given as 0 in the PSF (and stored as -1 in this program's internal  */
02395 /*  storage).                                                           */
02396 /************************************************************************/
02397 
02398 void Molecule::read_donors(FILE *fd)
02399 {
02400 #ifdef MEM_OPT_VERSION
02401     return;
02402 #else
02403   int d[2];               // temporary storage of donor atom index
02404   register int j;      // Loop counter
02405   int num_read=0;    // Number of bonds read so far
02406   int num_no_hydr=0;      // Number of bonds with no hydrogen given
02407 
02408   /*  Allocate the array to hold the bonds      */
02409   donors=new Bond[numDonors];
02410 
02411   if (donors == NULL)
02412   {
02413     NAMD_die("memory allocations failed in Molecule::read_donors");
02414   }
02415 
02416   /*  Loop through and read in all the donors      */
02417   while (num_read < numDonors)
02418   {
02419     /*  Loop and read in the two atom indexes    */
02420     for (j=0; j<2; j++)
02421     {
02422       /*  Read the atom number from the file.         */
02423       /*  Subtract 1 to convert the index from the    */
02424       /*  1 to NumAtoms used in the file to the       */
02425       /*  0 to NumAtoms-1 that we need    */
02426       d[j]=NAMD_read_int(fd, "DONORS")-1;
02427 
02428       /*  Check to make sure the index isn't too big  */
02429       if (d[j] >= numAtoms)
02430       {
02431         char err_msg[128];
02432 
02433         sprintf(err_msg,
02434     "DONOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE",
02435           d[j]+1, numAtoms, num_read+1);
02436         NAMD_die(err_msg);
02437       }
02438 
02439       /*  Check if there is a hydrogen given */
02440       if (d[j] < 0)
02441                           num_no_hydr++;
02442     }
02443 
02444     /*  Assign the atom indexes to the array element  */
02445     Bond *b = &(donors[num_read]);
02446     b->atom1=d[0];
02447     b->atom2=d[1];
02448 
02449     num_read++;
02450   }
02451 
02452   return;
02453 #endif
02454 }
02455 /*      END OF FUNCTION read_donors      */
02456 
02457 
02458 /************************************************************************/
02459 /*                  */
02460 /*      FUNCTION read_acceptors        */
02461 /*                  */
02462 /*  read_acceptors reads in the bond section of the .psf file.      */
02463 /*  This section contains a list of pairs of numbers where each pair is */
02464 /*  represents two atoms that are part of an H-bond.  Each atom pair is */
02465 /*  read in.                                                            */
02466 /*                  */
02467 /*  Acceptor atoms are the heavy atoms to which hydrogens directly      */
02468 /*  orient in a hydrogen bond interaction.  There will always be an     */
02469 /*  acceptor atom for each acceptor pair.  The antecedent atom, to      */
02470 /*  which the acceptor is bound, may not be given in the structure,     */
02471 /*  however, in which case the second atom index in the pair will be    */
02472 /*  given as 0 in the PSF (and stored as -1 in this program's internal  */
02473 /*  storage).                                                           */
02474 /************************************************************************/
02475 
02476 void Molecule::read_acceptors(FILE *fd)
02477 {
02478 #ifdef MEM_OPT_VERSION
02479     return;
02480 #else
02481   int d[2];               // temporary storage of atom index
02482   register int j;      // Loop counter
02483   int num_read=0;    // Number of bonds read so far
02484         int num_no_ante=0;      // number of pairs with no antecedent
02485 
02486   /*  Allocate the array to hold the bonds      */
02487   acceptors=new Bond[numAcceptors];
02488 
02489   if (acceptors == NULL)
02490   {
02491     NAMD_die("memory allocations failed in Molecule::read_acceptors");
02492   }
02493 
02494   /*  Loop through and read in all the acceptors      */
02495   while (num_read < numAcceptors)
02496   {
02497     /*  Loop and read in the two atom indexes    */
02498     for (j=0; j<2; j++)
02499     {
02500       /*  Read the atom number from the file.         */
02501       /*  Subtract 1 to convert the index from the    */
02502       /*  1 to NumAtoms used in the file to the       */
02503       /*  0 to NumAtoms-1 that we need    */
02504       d[j]=NAMD_read_int(fd, "ACCEPTORS")-1;
02505 
02506       /*  Check to make sure the index isn't too big  */
02507       if (d[j] >= numAtoms)
02508       {
02509         char err_msg[128];
02510 
02511         sprintf(err_msg, "ACCEPTOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE", d[j]+1, numAtoms, num_read+1);
02512         NAMD_die(err_msg);
02513       }
02514 
02515       /*  Check if there is an antecedent given */
02516       if (d[j] < 0)
02517                           num_no_ante++;
02518     }
02519 
02520     /*  Assign the atom indexes to the array element  */
02521     Bond *b = &(acceptors[num_read]);
02522     b->atom1=d[0];
02523     b->atom2=d[1];
02524 
02525     num_read++;
02526   }
02527 
02528   return;
02529 #endif
02530 }
02531 /*      END OF FUNCTION read_acceptors      */
02532 
02533 
02534 /************************************************************************/
02535 /*                  */
02536 /*      FUNCTION read_exclusions      */
02537 /*                  */
02538 /*   INPUTS:                */
02539 /*  fd - file descriptor for .psf file        */
02540 /*                  */
02541 /*  read_exclusions reads in the explicit non-bonded exclusions     */
02542 /*  from the .psf file.  This section is a little funky, so hang on.    */
02543 /*  Ok, first there is a list of atom indexes that is NumExclusions     */
02544 /*  long.  These are in some sense the atoms that will be exlcuded.     */
02545 /*  Following this list is a list of NumAtoms length that is a list     */
02546 /*  of indexes into the list of excluded atoms.  So an example.  Suppose*/
02547 /*  we have a 5 atom simulation with 3 explicit exclusions.  The .psf   */
02548 /*  file could look like:            */
02549 /*                  */
02550 /*  3!NNB                */
02551 /*  3 4 5                */
02552 /*  0 1 3 3 3              */
02553 /*                  */
02554 /*  This would mean that atom 1 has no explicit exclusions.  Atom 2     */
02555 /*  has an explicit exclusion with atom 3.  Atom 3 has an explicit      */
02556 /*  exclusion with atoms 4 AND 5.  And atoms 4 and 5 have no explicit   */
02557 /*  exclusions.  Got it!?!  I'm not sure who dreamed this up . . .      */
02558 /*                  */
02559 /************************************************************************/
02560 
02561 void Molecule::read_exclusions(FILE *fd)
02562 
02563 {
02564 #ifdef MEM_OPT_VERSION
02565     return;
02566 #else
02567   int *exclusion_atoms;  //  Array of indexes of excluded atoms
02568   register int num_read=0;    //  Number fo exclusions read in
02569   int current_index;  //  Current index value
02570   int last_index;    //  the previous index value
02571   register int insert_index=0;  //  index of where we are in exlcusions array
02572 
02573   /*  Allocate the array of exclusion structures and the array of */
02574   /*  exlcuded atom indexes          */
02575   exclusions      = new Exclusion[numExclusions];
02576   exclusion_atoms = new int[numExclusions];
02577 
02578   if ( (exclusions == NULL) || (exclusion_atoms == NULL) )
02579   {
02580     NAMD_die("memory allocation failed in Molecule::read_exclusions");
02581   }
02582 
02583   /*  First, read in the excluded atoms list      */
02584   for (num_read=0; num_read<numExclusions; num_read++)
02585   {
02586     /*  Read the atom number from the file. Subtract 1 to   */
02587     /*  convert the index from the 1 to NumAtoms used in the*/
02588     /*  file to the  0 to NumAtoms-1 that we need    */
02589     exclusion_atoms[num_read]=NAMD_read_int(fd, "IMPROPERS")-1;
02590 
02591     /*  Check for an illegal index        */
02592     if (exclusion_atoms[num_read] >= numAtoms)
02593     {
02594       char err_msg[128];
02595 
02596       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);
02597       NAMD_die(err_msg);
02598     }
02599   }
02600 
02601   /*  Now, go through and read the list of NumAtoms pointers into */
02602   /*  the array that we just read in        */
02603   last_index=0;
02604 
02605   for (num_read=0; num_read<numAtoms; num_read++)
02606   {
02607     /*  Read in the current index value      */
02608     current_index=NAMD_read_int(fd, "EXCLUSIONS");
02609 
02610     /*  Check for an illegal pointer      */
02611     if (current_index>numExclusions)
02612     {
02613       char err_msg[128];
02614 
02615       sprintf(err_msg, "EXCLUSION INDEX %d LARGER THAN NUMBER OF EXLCUSIONS %d IN PSF FILE, EXCLUSION #%d\n", 
02616          current_index+1, numExclusions, num_read);
02617       NAMD_die(err_msg);
02618     }
02619 
02620     /*  Check to see if it matches the last index.  If so   */
02621     /*  than this atom has no exclusions.  If not, then     */
02622     /*  we have to build some exclusions      */
02623     if (current_index != last_index)
02624     {
02625       /*  This atom has some exlcusions.  Loop from   */
02626       /*  the last_index to the current index.  This  */
02627       /*  will include how ever many exclusions this  */
02628       /*  atom has          */
02629       for (insert_index=last_index; 
02630            insert_index<current_index; insert_index++)
02631       {
02632         /*  Assign the two atoms involved.      */
02633         /*  The first one is our position in    */
02634         /*  the list, the second is based on    */
02635         /*  the pointer into the index list     */
02636         int a1 = num_read;
02637         int a2 = exclusion_atoms[insert_index];
02638         if ( a1 < a2 ) {
02639           exclusions[insert_index].atom1 = a1;
02640           exclusions[insert_index].atom2 = a2;
02641         } else if ( a2 < a1 ) {
02642           exclusions[insert_index].atom1 = a2;
02643           exclusions[insert_index].atom2 = a1;
02644         } else {
02645           char err_msg[128];
02646           sprintf(err_msg, "ATOM %d EXCLUDED FROM ITSELF IN PSF FILE\n", a1+1);
02647           NAMD_die(err_msg);
02648         }
02649       }
02650 
02651       last_index=current_index;
02652     }
02653   }
02654 
02655   /*  Free our temporary list of indexes        */
02656   delete [] exclusion_atoms;
02657 
02658   return;
02659 #endif
02660 }
02661 /*      END OF FUNCTION read_exclusions      */
02662 
02663 /************************************************************************/
02664 /*                  */
02665 /*        FUNCTION read_lphosts    */
02666 /*                  */
02667 /*   INPUTS:                */
02668 /*  fd - file pointer to the .psf file        */
02669 /*                  */
02670 /*  this function reads in the lone pair host section of the .psf file. */
02671 /*                  */
02672 void Molecule::read_lphosts(FILE *fd)
02673 {
02674   char buffer[512];  // Buffer for reading from file
02675   char lptype[8];
02676   int numhosts, index, i, read_count;
02677   Real distance, angle, dihedral;
02678 
02679   lphosts = new Lphost[numLphosts];
02680   if (lphosts == NULL)
02681   {
02682     NAMD_die("memory allocation failed in Molecule::read_lphosts");
02683   }
02684   for (i = 0;  i < numLphosts;  i++)
02685   {
02686     NAMD_read_line(fd, buffer);
02687     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
02688     read_count=sscanf(buffer, "%d %d %6s %f %f %f",
02689         &numhosts, &index, lptype, &distance, &angle, &dihedral);
02690     if (read_count != 6 || numhosts != 3 || index != 4*i + 1
02691         || strcmp(lptype,"F") != 0)
02692     {
02693       char err_msg[128];
02694       sprintf(err_msg, "BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n"
02695           "LINE=%s\n", i+1, buffer);
02696       NAMD_die(err_msg);
02697     }
02698     lphosts[i].distance = distance;
02699     lphosts[i].angle = angle * (M_PI/180);
02700     lphosts[i].dihedral = dihedral;
02701   }
02702   for (i = 0;  i < numLphosts;  i++) {
02703     lphosts[i].atom1 = NAMD_read_int(fd, "LPHOSTS")-1;
02704     lphosts[i].atom2 = NAMD_read_int(fd, "LPHOSTS")-1;
02705     lphosts[i].atom3 = NAMD_read_int(fd, "LPHOSTS")-1;
02706     lphosts[i].atom4 = NAMD_read_int(fd, "LPHOSTS")-1;
02707   }
02708 }
02709 /*      END OF FUNCTION read_lphosts    */
02710 
02711 /************************************************************************/
02712 /*                  */
02713 /*        FUNCTION read_anisos     */
02714 /*                  */
02715 /*   INPUTS:                */
02716 /*  fd - file pointer to the .psf file        */
02717 /*                  */
02718 /*  this function reads in the anisotropic terms section of .psf file. */
02719 /*                  */
02720 void Molecule::read_anisos(FILE *fd)
02721 {
02722   char buffer[512];  // Buffer for reading from file
02723   int numhosts, index, i, read_count;
02724   Real k11, k22, k33;
02725 
02726   anisos = new Aniso[numAnisos];
02727   if (anisos == NULL)
02728   {
02729     NAMD_die("memory allocation failed in Molecule::read_anisos");
02730   }
02731   for (i = 0;  i < numAnisos;  i++)
02732   {
02733     NAMD_read_line(fd, buffer);
02734     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
02735     read_count=sscanf(buffer, "%f %f %f", &k11, &k22, &k33);
02736     if (read_count != 3)
02737     {
02738       char err_msg[128];
02739       sprintf(err_msg, "BAD FORMAT FOR ANISO LINE %d IN PSF FILE LINE\n"
02740           "LINE=%s\n", i+1, buffer);
02741       NAMD_die(err_msg);
02742     }
02743     anisos[i].k11 = k11;
02744     anisos[i].k22 = k22;
02745     anisos[i].k33 = k33;
02746   }
02747   for (i = 0;  i < numAnisos;  i++) {
02748     anisos[i].atom1 = NAMD_read_int(fd, "ANISOS")-1;
02749     anisos[i].atom2 = NAMD_read_int(fd, "ANISOS")-1;
02750     anisos[i].atom3 = NAMD_read_int(fd, "ANISOS")-1;
02751     anisos[i].atom4 = NAMD_read_int(fd, "ANISOS")-1;
02752   }
02753 }
02754 /*      END OF FUNCTION read_anisos     */
02755 
02756 void Molecule::setOccupancyData(molfile_atom_t *atomarray){
02757     occupancy = new float[numAtoms];
02758     for(int i=0; i<numAtoms; i++) {
02759         occupancy[i] = atomarray[i].occupancy;
02760     }
02761 }
02762 
02763 void Molecule::setBFactorData(molfile_atom_t *atomarray){
02764     bfactor = new float[numAtoms];
02765     for(int i=0; i<numAtoms; i++) {
02766         bfactor[i] = atomarray[i].bfactor;
02767     }
02768 }
02769 
02770 void Molecule::plgLoadAtomBasics(molfile_atom_t *atomarray){
02771 #ifndef MEM_OPT_VERSION
02772     atoms = new Atom[numAtoms];
02773     atomNames = new AtomNameInfo[numAtoms];
02774     if(simParams->genCompressedPsf) {
02775         atomSegResids = new AtomSegResInfo[numAtoms];
02776     }    
02777     hydrogenGroup.resize(0);
02778 
02779     ResidueLookupElem *tmpResLookup = resLookup;
02780 
02781     for(int i=0; i<numAtoms; i++) {
02782         int reslength = strlen(atomarray[i].resname)+1;
02783         int namelength = strlen(atomarray[i].name)+1;
02784         int typelength = strlen(atomarray[i].type)+1;
02785         atomNames[i].resname = nameArena->getNewArray(reslength);
02786         atomNames[i].atomname = nameArena->getNewArray(namelength);
02787         atomNames[i].atomtype = nameArena->getNewArray(typelength);
02788         strcpy(atomNames[i].resname, atomarray[i].resname);
02789         strcpy(atomNames[i].atomname, atomarray[i].name);
02790         strcpy(atomNames[i].atomtype, atomarray[i].type);
02791 
02792         atoms[i].mass = atomarray[i].mass;
02793         atoms[i].charge = atomarray[i].charge;
02794         atoms[i].status = UnknownAtom;
02795 
02796         //add this atom to residue lookup table
02797         if(tmpResLookup) {
02798             tmpResLookup = tmpResLookup->append(atomarray[i].segid, atomarray[i].resid, i);
02799         }
02800 
02801         if(atomSegResids) { //for compressing molecule information
02802             AtomSegResInfo *one = atomSegResids + i;
02803             memcpy(one->segname, atomarray[i].segid, strlen(atomarray[i].segid)+1);
02804             one->resid = atomarray[i].resid;
02805         }
02806         //Determine the type of the atom
02807         if(atoms[i].mass <= 0.05) {
02808             atoms[i].status |= LonepairAtom;
02809         }else if(atoms[i].mass < 1.0) {
02810             atoms[i].status |= DrudeAtom;
02811         }else if(atoms[i].mass <= 3.5) {
02812             atoms[i].status |= HydrogenAtom;
02813         }else if((atomNames[i].atomname[0] == 'O') &&
02814                  (atoms[i].mass>=14.0) && (atoms[i].mass<=18.0)){
02815             atoms[i].status |= OxygenAtom;
02816         }
02817         //Lookk up the vdw constants for this atom
02818         params->assign_vdw_index(atomNames[i].atomtype, &atoms[i]);
02819     }
02820 #endif
02821 }
02822 
02823 void Molecule::plgLoadBonds(int *from, int *to){
02824 #ifndef MEM_OPT_VERSION
02825     bonds = new Bond[numBonds];
02826     char atom1name[11];
02827     char atom2name[11];
02828     int realNumBonds = 0;
02829     for(int i=0; i<numBonds; i++) {
02830         Bond *thisBond = bonds+realNumBonds;
02831         thisBond->atom1 = from[i]-1;
02832         thisBond->atom2 = to[i]-1;
02833         /* Get the atom type for the two atoms.
02834          * When we query the parameter object, we
02835          * need to send the atom type that is alphabetically
02836          * first as atom 1.
02837          */
02838         if(strcasecmp(atomNames[thisBond->atom1].atomtype,
02839                       atomNames[thisBond->atom2].atomtype)<0) {
02840             strcpy(atom1name, atomNames[thisBond->atom1].atomtype);
02841             strcpy(atom2name, atomNames[thisBond->atom2].atomtype);
02842         }else{
02843             strcpy(atom2name, atomNames[thisBond->atom1].atomtype);
02844             strcpy(atom1name, atomNames[thisBond->atom2].atomtype);
02845         }
02846         params->assign_bond_index(atom1name, atom2name, thisBond);
02847 
02848         //Make sure this isn't a fake bond meant for shake in x-plor
02849         Real k, x0;
02850         params->get_bond_params(&k, &x0, thisBond->bond_type);
02851         if(simParams->watmodel == WAT_SWM4) {
02852             //need to retain Lonepair bonds for Drude
02853             if(k!=0. || is_lp(thisBond->atom1) || 
02854                is_lp(thisBond->atom2)) {               
02855                 realNumBonds++;
02856             }
02857         }else{
02858             if(k != 0.) realNumBonds++;
02859         }
02860     }
02861 
02862     if(numBonds != realNumBonds) {
02863         iout << iWARN << "Ignored" << numBonds-realNumBonds <<
02864             "bonds with zero force constants.\n" <<endi;
02865         iout << iWARN << "Will get H-H distance in rigid H20 from H-O-H angle.\n" <<endi;
02866     }
02867     numBonds = realNumBonds;
02868 #endif
02869 }
02870 
02871 void Molecule::plgLoadAngles(int *plgAngles)
02872 {    
02873 #ifndef MEM_OPT_VERSION
02874     char atom1name[11];
02875     char atom2name[11];
02876     char atom3name[11];
02877 
02878     angles=new Angle[numAngles];
02879     int *atomid = plgAngles;
02880     int numRealAngles = 0;
02881     for(int i=0; i<numAngles; i++) {
02882         Angle *thisAngle = angles+numRealAngles;
02883         thisAngle->atom1 = atomid[0]-1;
02884         thisAngle->atom2 = atomid[1]-1;
02885         thisAngle->atom3 = atomid[2]-1;
02886         atomid += 3;
02887 
02888         if(strcasecmp(atomNames[thisAngle->atom1].atomtype,
02889                       atomNames[thisAngle->atom2].atomtype)<0) {
02890             strcpy(atom1name, atomNames[thisAngle->atom1].atomtype);
02891             strcpy(atom2name, atomNames[thisAngle->atom2].atomtype);
02892             strcpy(atom3name, atomNames[thisAngle->atom3].atomtype);
02893         }else{
02894             strcpy(atom1name, atomNames[thisAngle->atom3].atomtype);
02895             strcpy(atom2name, atomNames[thisAngle->atom2].atomtype);
02896             strcpy(atom3name, atomNames[thisAngle->atom1].atomtype);
02897         }
02898 
02899         params->assign_angle_index(atom1name, atom2name, atom3name, thisAngle);
02900 
02901         Real k, t0, k_ub, r_ub;
02902         params->get_angle_params(&k, &t0, &k_ub, &r_ub, thisAngle->angle_type);
02903         if(k!=0. || k_ub!=0.) numRealAngles++;
02904     }
02905 
02906     if(numAngles != numRealAngles) {
02907         iout << iWARN << "Ignored" << numAngles-numRealAngles << 
02908             " angles with zero force constants.\n" << endi; 
02909     }
02910     numAngles = numRealAngles;
02911 #endif
02912 }
02913 
02914 void Molecule::plgLoadDihedrals(int *plgDihedrals)
02915 {
02916 #ifndef MEM_OPT_VERSION
02917     char atom1name[11];
02918     char atom2name[11];
02919     char atom3name[11];
02920     char atom4name[11];
02921     int lastAtomIds[4];
02922     int multiplicity = 1; //multiplicity of the current bond
02923 
02924     lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
02925     dihedrals = new Dihedral[numDihedrals];
02926     int numRealDihedrals = 0;
02927     int *atomid = plgDihedrals;
02928     for(int i=0; i<numDihedrals; i++, atomid+=4) {
02929         Dihedral *thisDihedral = dihedrals + numRealDihedrals;
02930         Bool duplicate_bond = TRUE;
02931         for(int j=0; j<4; j++) {
02932             if(atomid[j] != lastAtomIds[j]) {
02933                 duplicate_bond = FALSE;
02934             }
02935             lastAtomIds[j] = atomid[j];
02936         }
02937 
02938         strcpy(atom1name, atomNames[atomid[0]-1].atomtype);
02939         strcpy(atom2name, atomNames[atomid[1]-1].atomtype);
02940         strcpy(atom3name, atomNames[atomid[2]-1].atomtype);
02941         strcpy(atom4name, atomNames[atomid[3]-1].atomtype);
02942 
02943         if(duplicate_bond) {
02944             multiplicity++;
02945             if(multiplicity==2) {
02946                 numMultipleDihedrals++;
02947             }
02948         }else{
02949             multiplicity=1;
02950             numRealDihedrals++;
02951         }
02952 
02953         params->assign_dihedral_index(atom1name, atom2name,
02954                                       atom3name, atom4name, thisDihedral,
02955                                       multiplicity);
02956         thisDihedral->atom1 = atomid[0]-1;
02957         thisDihedral->atom2 = atomid[1]-1;
02958         thisDihedral->atom3 = atomid[2]-1;
02959         thisDihedral->atom4 = atomid[3]-1;
02960     }
02961 
02962     numDihedrals = numRealDihedrals;
02963 #endif
02964 }
02965 
02966 void Molecule::plgLoadImpropers(int *plgImpropers)
02967 {
02968 #ifndef MEM_OPT_VERSION
02969     char atom1name[11];
02970     char atom2name[11];
02971     char atom3name[11];
02972     char atom4name[11];
02973     int lastAtomIds[4];
02974     int multiplicity = 1; //multiplicity of the current bond
02975 
02976     lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
02977     impropers = new Improper[numImpropers];
02978     int numRealImpropers = 0;
02979     int *atomid = plgImpropers;
02980     for(int i=0; i<numImpropers; i++, atomid+=4) {
02981         Improper *thisImproper = impropers + numRealImpropers;
02982         Bool duplicate_bond = TRUE;
02983         for(int j=0; j<4; j++) {
02984             if(atomid[j] != lastAtomIds[j]) {
02985                 duplicate_bond = FALSE;
02986             }
02987             lastAtomIds[j] = atomid[j];
02988         }
02989 
02990         strcpy(atom1name, atomNames[atomid[0]-1].atomtype);
02991         strcpy(atom2name, atomNames[atomid[1]-1].atomtype);
02992         strcpy(atom3name, atomNames[atomid[2]-1].atomtype);
02993         strcpy(atom4name, atomNames[atomid[3]-1].atomtype);
02994 
02995         if(duplicate_bond) {
02996             multiplicity++;
02997             if(multiplicity==2) {
02998                 numMultipleImpropers++;
02999             }
03000         }else{
03001             multiplicity=1;
03002             numRealImpropers++;
03003         }
03004 
03005         params->assign_improper_index(atom1name, atom2name,
03006                                       atom3name, atom4name, thisImproper,
03007                                       multiplicity);
03008         thisImproper->atom1 = atomid[0]-1;
03009         thisImproper->atom2 = atomid[1]-1;
03010         thisImproper->atom3 = atomid[2]-1;
03011         thisImproper->atom4 = atomid[3]-1;
03012     }
03013 
03014     numImpropers = numRealImpropers;
03015 #endif
03016 }
03017 
03018 void Molecule::plgLoadCrossterms(int *plgCterms)
03019 {
03020 #ifndef MEM_OPT_VERSION
03021     char atom1name[11];
03022     char atom2name[11];
03023     char atom3name[11];
03024     char atom4name[11];
03025     char atom5name[11];
03026     char atom6name[11];
03027     char atom7name[11];
03028     char atom8name[11];
03029     int lastAtomIds[8];    
03030 
03031     for(int i=0; i<8; i++)
03032         lastAtomIds[i]=-1;
03033     
03034     crossterms = new Crossterm[numCrossterms];
03035     int numRealCrossterms = 0;
03036     int *atomid = plgCterms;
03037     for(int i=0; i<numCrossterms; i++, atomid+=8) {
03038         Crossterm *thisCrossterm = crossterms + numRealCrossterms;
03039         Bool duplicate_bond = TRUE;
03040         for(int j=0; j<8; j++) {
03041             if(atomid[j] != lastAtomIds[j]) {
03042                 duplicate_bond = FALSE;
03043             }
03044             lastAtomIds[j] = atomid[j];
03045         }
03046 
03047         strcpy(atom1name, atomNames[atomid[0]-1].atomtype);
03048         strcpy(atom2name, atomNames[atomid[1]-1].atomtype);
03049         strcpy(atom3name, atomNames[atomid[2]-1].atomtype);
03050         strcpy(atom4name, atomNames[atomid[3]-1].atomtype);
03051         strcpy(atom5name, atomNames[atomid[4]-1].atomtype);
03052         strcpy(atom6name, atomNames[atomid[5]-1].atomtype);
03053         strcpy(atom7name, atomNames[atomid[6]-1].atomtype);
03054         strcpy(atom8name, atomNames[atomid[7]-1].atomtype);
03055 
03056         if(duplicate_bond) {
03057             iout << iWARN <<"Duplicate cross-term detected.\n" << endi;
03058         }else
03059             numRealCrossterms++;
03060 
03061         params->assign_crossterm_index(atom1name, atom2name,
03062                                        atom3name, atom4name, atom5name,
03063                                        atom6name, atom7name, atom8name,
03064                                        thisCrossterm);
03065 
03066         thisCrossterm->atom1 = atomid[0]-1;
03067         thisCrossterm->atom2 = atomid[1]-1;
03068         thisCrossterm->atom3 = atomid[2]-1;
03069         thisCrossterm->atom4 = atomid[3]-1;
03070         thisCrossterm->atom5 = atomid[4]-1;
03071         thisCrossterm->atom6 = atomid[5]-1;
03072         thisCrossterm->atom7 = atomid[6]-1;
03073         thisCrossterm->atom8 = atomid[7]-1;
03074     }
03075 
03076     numCrossterms = numRealCrossterms;
03077 #endif
03078 }
03079 
03080 
03081 /************************************************************************/
03082 /*                  */
03083 /*      FUNCTION print_atoms        */
03084 /*                  */
03085 /*  print_atoms prints out the list of atoms stored in this object. */
03086 /*  It is inteded mainly for debugging purposes.      */
03087 /*                  */
03088 /************************************************************************/
03089 
03090 void Molecule::print_atoms(Parameters *params)
03091 
03092 {
03093   register int i;
03094   Real sigma;
03095   Real epsilon;
03096   Real sigma14;
03097   Real epsilon14;
03098 
03099   DebugM(2,"ATOM LIST\n" \
03100       << "******************************************\n" \
03101                   << "NUM  NAME TYPE RES  MASS    CHARGE CHARGE   FEP-CHARGE"  \
03102       << "SIGMA   EPSILON SIGMA14 EPSILON14\n" \
03103         << endi);
03104 
03105   for (i=0; i<numAtoms; i++)
03106   {
03107     params->get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, 
03108         atoms[i].vdw_type);
03109 
03110 #ifdef MEM_OPT_VERSION
03111      DebugM(2,i+1 << " " << atomNamePool[atomNames[i].atomnameIdx]  \
03112               << " " << atomTypePool[atomNames[i].atomtypeIdx] << " " \
03113               << resNamePool[atomNames[i].resnameIdx]  << " "  \
03114             << atommass(i)  \
03115         << " " << atomcharge(i) << " " << sigma \
03116         << " " << epsilon << " " << sigma14 \
03117         << " " << epsilon14 << "\n" \
03118         << endi);
03119 #else
03120     DebugM(2,i+1 << " " << atomNames[i].atomname  \
03121               << " " << atomNames[i].atomtype << " " \
03122               << atomNames[i].resname  << " " << atoms[i].mass  \
03123         << " " << atoms[i].charge << " " << sigma \
03124         << " " << epsilon << " " << sigma14 \
03125         << " " << epsilon14 << "\n" \
03126         << endi);
03127 #endif
03128   }
03129 }
03130 /*      END OF FUNCTION print_atoms      */
03131 
03132 /************************************************************************/
03133 /*                  */
03134 /*      FUNCTION print_bonds        */
03135 /*                  */
03136 /*  print_bonds prints out the list of bonds stored in this object. */
03137 /*  It is inteded mainly for debugging purposes.      */
03138 /*                  */
03139 /************************************************************************/
03140 
03141 void Molecule::print_bonds(Parameters *params)
03142 
03143 {
03144   register int i;
03145   Real k;
03146   Real x0;
03147 
03148   DebugM(2,"BOND LIST\n" << "********************************\n" \
03149       << "ATOM1 ATOM2 TYPE1 TYPE2      k        x0" \
03150       << endi);
03151 
03152   #ifdef MEM_OPT_VERSION
03153   for(i=0; i<numAtoms; i++){
03154       AtomSignature *sig = &atomSigPool[eachAtomSig[i]];
03155       int bCnt = sig->bondCnt;
03156       TupleSignature *bSigs = sig->bondSigs;
03157       for(int j=0; j<bCnt; j++){
03158           Bond aBond;
03159           aBond.atom1 = i;
03160           aBond.atom2 = i+bSigs[j].offset[0];
03161           aBond.bond_type = bSigs[j].tupleParamType;
03162           params->get_bond_params(&k, &x0, aBond.bond_type);
03163 
03164           DebugM(2,aBond.atom1+1 << " " \
03165              << aBond.atom2+1 << " "   \
03166              << get_atomtype(aBond.atom1) << " "  \
03167              << get_atomtype(aBond.atom2) << " " << k \
03168              << " " << x0 << endi);
03169       }
03170   }
03171   #else
03172   for (i=0; i<numBonds; i++)
03173   {
03174     params->get_bond_params(&k, &x0, bonds[i].bond_type);
03175 
03176     DebugM(2,bonds[i].atom1+1 << " " \
03177        << bonds[i].atom2+1 << " "   \
03178        << atomNames[bonds[i].atom1].atomtype << " "  \
03179        << atomNames[bonds[i].atom2].atomtype << " " << k \
03180        << " " << x0 << endi);
03181   }
03182   #endif
03183 }
03184 /*      END OF FUNCTION print_bonds      */
03185 
03186 /************************************************************************/
03187 /*                  */
03188 /*      FUNCTION print_exclusions      */
03189 /*                  */
03190 /*  print_exlcusions prints out the list of exlcusions stored in    */
03191 /*  this object.  It is inteded mainly for debugging purposes.    */
03192 /*                  */
03193 /************************************************************************/
03194 
03195 void Molecule::print_exclusions()
03196 {
03197 #ifdef MEM_OPT_VERSION
03198     DebugM(2, "WARNING: this function is not availabe in memory optimized version!\n" << endi);
03199 #else
03200   register int i;
03201 
03202   DebugM(2,"EXPLICIT EXCLUSION LIST\n" \
03203       << "********************************\n" \
03204             << "ATOM1 ATOM2 " \
03205       << endi);
03206 
03207   for (i=0; i<numExclusions; i++)
03208   {
03209     DebugM(2,exclusions[i].atom1+1 << "  " \
03210        << exclusions[i].atom2+1 << endi);
03211   }
03212 #endif
03213 }
03214 /*      END OF FUNCTION print_exclusions    */
03215 
03216 /************************************************************************/
03217 /*                  */
03218 /*      FUNCTION send_Molecule        */
03219 /*                  */
03220 /*  send_Molecule is used by the Master node to distribute the      */
03221 /*   structural information to all the client nodes.  It is NEVER called*/
03222 /*   by the client nodes.              */
03223 /*                  */
03224 /************************************************************************/
03225 
03226 
03227 void Molecule::send_Molecule(MOStream *msg)
03228 
03229 {
03230     #ifdef MEM_OPT_VERSION
03231     //  Now build arrays of indexes into these arrays by atom      
03232     //generating all the new atom signatures and exclusion signatures only on pe0
03233     build_lists_by_atom();
03234     #endif
03235       
03236   /*//  Message to send to clients
03237   int bufSize = BUFSIZE;
03238   // When the simulation system is very large, then the buffer size should be expanded to reduce the number of one-to-all broadcasts.
03239   if(numAtoms>=1000000) bufSize=16*BUFSIZE;
03240   MOStream *msg=com_obj->newOutputStream(ALLBUTME, MOLECULETAG, bufSize);
03241   if ( msg == NULL )
03242   {
03243     NAMD_die("Memory allocation failed in Molecule::send_Molecule");
03244   }*/
03245 
03246   #ifdef MEM_OPT_VERSION
03247       msg->put(numAtoms);
03248       //mass and charge pool needed to be sent to other processors
03249       //for the sake of function call: build_atom_status
03250       msg->put(massPoolSize);
03251       msg->put(massPoolSize, atomMassPool);
03252       msg->put(numAtoms*sizeof(Index), (char *)eachAtomMass);
03253 
03254       msg->put(chargePoolSize);
03255       msg->put(chargePoolSize, atomChargePool);
03256       msg->put(numAtoms*sizeof(Index), (char *)eachAtomCharge);
03257 
03258       //vdw_type, partner etc.
03259       msg->put(numAtoms*sizeof(AtomCstInfo), (char *)atoms);
03260 
03261       //put atoms' signatures
03262       msg->put(atomSigPoolSize);
03263       for(int i=0; i<atomSigPoolSize; i++)
03264           atomSigPool[i].pack(msg);
03265 
03266       //put atom's exclusion signatures
03267       msg->put(exclSigPoolSize);
03268       for(int i=0; i<exclSigPoolSize; i++)
03269           exclSigPool[i].pack(msg);
03270 
03271       //put eachAtomSig and eachAtomExclSig
03272       msg->put(numAtoms*sizeof(Index), (char *)eachAtomSig);
03273       msg->put(numAtoms*sizeof(Index), (char *)eachAtomExclSig);
03274 
03275       msg->put(numAtoms*sizeof(int32), (char *)clusterSigs);
03276   #else
03277       msg->put(numAtoms);
03278       msg->put(numAtoms*sizeof(Atom), (char*)atoms);
03279   #endif
03280       //  Send the bond information
03281       msg->put(numRealBonds);
03282       msg->put(numBonds);
03283 
03284       #ifndef MEM_OPT_VERSION
03285       if (numBonds)
03286       {
03287         msg->put(numBonds*sizeof(Bond), (char*)bonds);
03288       }
03289       #endif
03290 
03291       //  Send the angle information
03292       msg->put(numAngles);
03293 
03294       #ifndef MEM_OPT_VERSION
03295       if (numAngles)
03296       {
03297         msg->put(numAngles*sizeof(Angle), (char*)angles);
03298       }
03299       #endif
03300 
03301       //  Send the dihedral information
03302       msg->put(numDihedrals);
03303 
03304       #ifndef MEM_OPT_VERSION
03305       if (numDihedrals)
03306       {
03307         msg->put(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
03308       }
03309       #endif
03310 
03311       //  Send the improper information
03312       msg->put(numImpropers);
03313 
03314       #ifndef MEM_OPT_VERSION
03315       if (numImpropers)
03316       {
03317         msg->put(numImpropers*sizeof(Improper), (char*)impropers);
03318       }
03319       #endif
03320 
03321       //  Send the crossterm information
03322       msg->put(numCrossterms);
03323 
03324       #ifndef MEM_OPT_VERSION
03325       if (numCrossterms)
03326       {
03327         msg->put(numCrossterms*sizeof(Crossterm), (char*)crossterms);
03328       }
03329       #endif
03330 
03331       // send the hydrogen bond donor information
03332       msg->put(numDonors);
03333 
03334       if(numDonors)
03335       {
03336         msg->put(numDonors*sizeof(Bond), (char*)donors);
03337       }
03338 
03339       // send the hydrogen bond acceptor information
03340       msg->put(numAcceptors);
03341 
03342       if(numAcceptors)
03343       {
03344         msg->put(numAcceptors*sizeof(Bond), (char*)acceptors);
03345       }
03346 
03347       //  Send the exclusion information
03348       #ifndef MEM_OPT_VERSION
03349       msg->put(numExclusions);
03350 
03351       if (numExclusions)
03352       {
03353         msg->put(numExclusions*sizeof(Exclusion), (char*)exclusions);
03354       }
03355       #endif
03356 
03357       //hydrogen group info is calculated when generating the compressed molecule
03358       //information, so this has to be distributed to other nodes instead of
03359       //being recalculated in build_atom_status in receive_Molecule function
03360       #ifdef MEM_OPT_VERSION
03361       msg->put(numHydrogenGroups);      
03362       msg->put(numAtoms*sizeof(HydrogenGroupID), (char *)hydrogenGroup.begin());      
03363       #endif
03364       
03365       //  Send the constraint information, if used
03366       if (simParams->constraintsOn)
03367       {
03368          msg->put(numConstraints);
03369          
03370          msg->put(numAtoms, consIndexes);
03371          
03372          if (numConstraints)
03373          {
03374            msg->put(numConstraints*sizeof(ConstraintParams), (char*)consParams);
03375          }
03376       }
03377       
03378       /* BEGIN gf */
03379       // Send the gridforce information, if used
03380       if (simParams->mgridforceOn)
03381       {
03382          DebugM(3, "Sending gridforce info\n");
03383          msg->put(numGridforceGrids);
03384          
03385          for (int grid = 0; grid < numGridforceGrids; grid++) {
03386              msg->put(numGridforces[grid]);
03387              msg->put(numAtoms, gridfrcIndexes[grid]);
03388              if (numGridforces[grid])
03389              {
03390                  msg->put(numGridforces[grid]*sizeof(GridforceParams), (char*)gridfrcParams[grid]);
03391              }
03392              gridfrcGrid[grid]->pack(msg);      // grid object writes its private data to message itself
03393          }
03394       }
03395       /* END gf */
03396 
03397       //  Send the stirring information, if used
03398       if (simParams->stirOn)
03399       {
03400                //CkPrintf ("DEBUG: putting numStirredAtoms..\n");
03401          msg->put(numStirredAtoms);
03402                //CkPrintf ("DEBUG: putting numAtoms,stirIndexes.. numAtoms=%d\n",numStirredAtoms);
03403          msg->put(numAtoms, stirIndexes);
03404                //CkPrintf ("DEBUG: if numStirredAtoms..\n");
03405          if (numStirredAtoms)
03406          {
03407            //CkPrintf ("DEBUG: big put, with (char*)stirParams\n");
03408            msg->put(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
03409          }
03410       }
03411       
03412       
03413       //  Send the moving drag information, if used
03414       if (simParams->movDragOn) {
03415          msg->put(numMovDrag);
03416          msg->put(numAtoms, movDragIndexes);
03417          if (numMovDrag)
03418          {
03419            msg->put(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
03420          }
03421       }
03422       
03423       //  Send the rotating drag information, if used
03424       if (simParams->rotDragOn) {
03425          msg->put(numRotDrag);
03426          msg->put(numAtoms, rotDragIndexes);
03427          if (numRotDrag)
03428          {
03429            msg->put(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
03430          }
03431       }
03432       
03433       //  Send the "constant" torque information, if used
03434       if (simParams->consTorqueOn) {
03435          msg->put(numConsTorque);
03436          msg->put(numAtoms, consTorqueIndexes);
03437          if (numConsTorque)
03438          {
03439            msg->put(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
03440          }
03441       }
03442       
03443       // Send the constant force information, if used
03444       if (simParams->consForceOn)
03445       { msg->put(numConsForce);
03446         msg->put(numAtoms, consForceIndexes);
03447         if (numConsForce)
03448           msg->put(numConsForce*sizeof(Vector), (char*)consForce);
03449       }
03450 
03451       //  Send the langevin parameters, if active
03452       if (simParams->langevinOn || simParams->tCoupleOn)
03453       {
03454         msg->put(numAtoms, langevinParams);
03455       }
03456 
03457       //  Send fixed atoms, if active
03458       if (simParams->fixedAtomsOn)
03459       {
03460         msg->put(numFixedAtoms);
03461         msg->put(numAtoms, fixedAtomFlags);
03462         msg->put(numFixedRigidBonds);
03463       }
03464 
03465       if (simParams->excludeFromPressure) {
03466         msg->put(numExPressureAtoms);
03467         msg->put(numAtoms, exPressureAtomFlags);
03468       }
03469 
03470 //fepb
03471       // send fep atom info
03472       if (simParams->alchFepOn || simParams->alchThermIntOn || simParams->lesOn || simParams->pairInteractionOn) {
03473         msg->put(numFepInitial);
03474         msg->put(numFepFinal);
03475         msg->put(numAtoms*sizeof(char), (char*)fepAtomFlags);
03476       }
03477 //fepe
03478 
03479       // Broadcast the message to the other nodes
03480       msg->end();
03481       delete msg;
03482 
03483       #ifdef MEM_OPT_VERSION
03484       build_excl_check_signatures();
03485       #else      
03486       //  Now build arrays of indexes into these arrays by atom      
03487       build_lists_by_atom();
03488       #endif
03489 
03490 }
03491     /*      END OF FUNCTION send_Molecule      */
03492 
03493     /************************************************************************/
03494     /*                  */
03495     /*      FUNCTION receive_Molecule      */
03496     /*                  */
03497     /*  receive_Molecule is used by all the clients to receive the  */
03498     /*   structural data sent out by the master node.  It is NEVER called   */
03499     /*   by the Master node.            */
03500     /*                  */
03501     /************************************************************************/
03502 
03503 void Molecule::receive_Molecule(MIStream *msg)
03504 {
03505       //  Get the atom information
03506       msg->get(numAtoms);
03507 
03508   #ifdef MEM_OPT_VERSION   
03509       //are mass and charge pool needed to be sent to other processors???
03510       msg->get(massPoolSize);
03511       if(atomMassPool) delete [] atomMassPool;
03512       atomMassPool = new Real[massPoolSize];
03513       msg->get(massPoolSize, atomMassPool);
03514       if(eachAtomMass) delete [] eachAtomMass;
03515       eachAtomMass = new Index[numAtoms];
03516       msg->get(numAtoms*sizeof(Index), (char *)eachAtomMass);
03517       
03518       msg->get(chargePoolSize);
03519       if(atomChargePool) delete [] atomChargePool;
03520       atomChargePool = new Real[chargePoolSize];
03521       msg->get(chargePoolSize, atomChargePool);
03522       if(eachAtomCharge) delete [] eachAtomCharge;
03523       eachAtomCharge = new Index[numAtoms];
03524       msg->get(numAtoms*sizeof(Index), (char *)eachAtomCharge);
03525 
03526       //vdw_type, partner etc.
03527       if(atoms) delete [] atoms;
03528       atoms = new AtomCstInfo[numAtoms];
03529       msg->get(numAtoms*sizeof(AtomCstInfo), (char *)atoms);
03530 
03531       //get atoms' signatures
03532       msg->get(atomSigPoolSize);
03533       if(atomSigPool) delete [] atomSigPool;
03534       atomSigPool = new AtomSignature[atomSigPoolSize];
03535       for(int i=0; i<atomSigPoolSize; i++)
03536           atomSigPool[i].unpack(msg);
03537 
03538       //get exclusions' signatures
03539       msg->get(exclSigPoolSize);
03540       if(exclSigPool) delete [] exclSigPool;
03541       exclSigPool = new ExclusionSignature[exclSigPoolSize];
03542       for(int i=0; i<exclSigPoolSize; i++)
03543           exclSigPool[i].unpack(msg);
03544 
03545       //get eachAtomSig and eachAtomExclSig
03546       if(eachAtomSig) delete [] eachAtomSig;
03547       eachAtomSig = new Index[numAtoms];
03548       msg->get(numAtoms*sizeof(Index), (char *)eachAtomSig);
03549       if(eachAtomExclSig) delete [] eachAtomExclSig;
03550       eachAtomExclSig = new Index[numAtoms];
03551       msg->get(numAtoms*sizeof(Index), (char *)eachAtomExclSig);
03552 
03553       if(clusterSigs) delete [] clusterSigs;
03554       clusterSigs = new int32[numAtoms];
03555       msg->get(numAtoms*sizeof(int32), (char *)clusterSigs); 
03556   #else
03557       delete [] atoms;
03558       atoms= new Atom[numAtoms];
03559 
03560       msg->get(numAtoms*sizeof(Atom), (char*)atoms);
03561   #endif
03562 
03563       //  Get the bond information
03564       msg->get(numRealBonds);
03565       msg->get(numBonds);
03566 
03567       #ifndef MEM_OPT_VERSION
03568       if (numBonds)
03569       {
03570         delete [] bonds;
03571         bonds=new Bond[numBonds];
03572 
03573         msg->get(numBonds*sizeof(Bond), (char*)bonds);
03574       }
03575       #endif
03576 
03577       //  Get the angle information
03578       msg->get(numAngles);
03579 
03580       #ifndef MEM_OPT_VERSION
03581       if (numAngles)
03582       {
03583         delete [] angles;
03584         angles=new Angle[numAngles];
03585 
03586         msg->get(numAngles*sizeof(Angle), (char*)angles);
03587       }
03588       #endif
03589 
03590       //  Get the dihedral information
03591       msg->get(numDihedrals);
03592 
03593       #ifndef MEM_OPT_VERSION
03594       if (numDihedrals)
03595       {
03596         delete [] dihedrals;
03597         dihedrals=new Dihedral[numDihedrals];
03598 
03599         msg->get(numDihedrals*sizeof(Dihedral), (char*)dihedrals);
03600       }
03601       #endif
03602 
03603       //  Get the improper information
03604       msg->get(numImpropers);
03605 
03606       #ifndef MEM_OPT_VERSION
03607       if (numImpropers)
03608       {
03609         delete [] impropers;
03610         impropers=new Improper[numImpropers];
03611 
03612         msg->get(numImpropers*sizeof(Improper), (char*)impropers);
03613       }
03614       #endif
03615 
03616       //  Get the crossterm information
03617       msg->get(numCrossterms);
03618 
03619       #ifndef MEM_OPT_VERSION
03620       if (numCrossterms)
03621       {
03622         delete [] crossterms;
03623         crossterms=new Crossterm[numCrossterms];
03624 
03625         msg->get(numCrossterms*sizeof(Crossterm), (char*)crossterms);
03626       }
03627       #endif
03628 
03629       //  Get the hydrogen bond donors
03630       msg->get(numDonors);
03631 
03632       if (numDonors)
03633       {
03634         delete [] donors;
03635         donors=new Bond[numDonors];
03636 
03637         msg->get(numDonors*sizeof(Bond), (char*)donors);
03638       }
03639 
03640       //  Get the hydrogen bond acceptors
03641       msg->get(numAcceptors);
03642 
03643       if (numAcceptors)
03644       {
03645         delete [] acceptors;
03646         acceptors=new Bond[numAcceptors];
03647 
03648         msg->get(numAcceptors*sizeof(Bond), (char*)acceptors);
03649       }
03650 
03651       //  Get the exclusion information
03652       #ifndef MEM_OPT_VERSION
03653       msg->get(numExclusions);
03654 
03655       if (numExclusions)
03656       {
03657         delete [] exclusions;
03658         exclusions=new Exclusion[numExclusions];
03659 
03660         msg->get(numExclusions*sizeof(Exclusion), (char*)exclusions);
03661       }
03662       #endif
03663 
03664       #ifdef MEM_OPT_VERSION
03665       msg->get(numHydrogenGroups);      
03666       hydrogenGroup.resize(numAtoms);
03667       msg->get(numAtoms*sizeof(HydrogenGroupID), (char *)hydrogenGroup.begin());     
03668       #endif
03669       
03670       //  Get the constraint information, if they are active
03671       if (simParams->constraintsOn)
03672       {
03673          msg->get(numConstraints);
03674 
03675          delete [] consIndexes;
03676          consIndexes = new int32[numAtoms];
03677          
03678          msg->get(numAtoms, consIndexes);
03679          
03680          if (numConstraints)
03681          {
03682            delete [] consParams;
03683            consParams = new ConstraintParams[numConstraints];
03684       
03685            msg->get(numConstraints*sizeof(ConstraintParams), (char*)consParams);
03686          }
03687       }
03688 
03689       /* BEGIN gf */
03690       if (simParams->mgridforceOn)
03691       {
03692          DebugM(3, "Receiving gridforce info\n");
03693          
03694          msg->get(numGridforceGrids);
03695          
03696          delete [] numGridforces;
03697          numGridforces = new int[numGridforceGrids];
03698          
03699          delete [] gridfrcIndexes;      // Should I be deleting elements of these first?
03700          delete [] gridfrcParams;
03701          delete [] gridfrcGrid;
03702          gridfrcIndexes = new int32*[numGridforceGrids];
03703          gridfrcParams = new GridforceParams*[numGridforceGrids];
03704          gridfrcGrid = new GridforceGrid*[numGridforceGrids];
03705          
03706          for (int grid = 0; grid < numGridforceGrids; grid++) {
03707              msg->get(numGridforces[grid]);
03708              
03709              gridfrcIndexes[grid] = new int32[numAtoms];
03710              msg->get(numAtoms, gridfrcIndexes[grid]);
03711          
03712              if (numGridforces[grid])
03713              {
03714                  gridfrcParams[grid] = new GridforceParams[numGridforces[grid]];
03715                  msg->get(numGridforces[grid]*sizeof(GridforceParams), (char*)gridfrcParams[grid]);
03716              }
03717              
03718              gridfrcGrid[grid] = new GridforceGrid(grid);
03719              gridfrcGrid[grid]->unpack(msg);
03720              CProxy_ComputeGridForceNodeMgr 
03721                mgr(CkpvAccess(BOCclass_group).computeGridForceNodeMgr);
03722              GridDepositMsg *outmsg = new GridDepositMsg;
03723              outmsg->gridnum = grid;
03724              outmsg->grid = gridfrcGrid[grid];
03725              outmsg->num_grids = numGridforceGrids;
03726              mgr[CkMyNode()].depositInitialGrid(outmsg);
03727          }
03728       }
03729       /* END gf */
03730       
03731       //  Get the stirring information, if stirring is  active
03732       if (simParams->stirOn)
03733       {
03734          msg->get(numStirredAtoms);
03735 
03736          delete [] stirIndexes;
03737          stirIndexes = new int32[numAtoms];
03738          
03739          msg->get(numAtoms, stirIndexes);
03740          
03741          if (numStirredAtoms)
03742          {
03743            delete [] stirParams;
03744            stirParams = new StirParams[numStirredAtoms];
03745       
03746            msg->get(numStirredAtoms*sizeof(StirParams), (char*)stirParams);
03747          }
03748       }
03749       
03750       //  Get the moving drag information, if it is active
03751       if (simParams->movDragOn) {
03752          msg->get(numMovDrag);
03753          delete [] movDragIndexes;
03754          movDragIndexes = new int32[numAtoms];
03755          msg->get(numAtoms, movDragIndexes);
03756          if (numMovDrag)
03757          {
03758            delete [] movDragParams;
03759            movDragParams = new MovDragParams[numMovDrag];
03760            msg->get(numMovDrag*sizeof(MovDragParams), (char*)movDragParams);
03761          }
03762       }
03763       
03764       //  Get the rotating drag information, if it is active
03765       if (simParams->rotDragOn) {
03766          msg->get(numRotDrag);
03767          delete [] rotDragIndexes;
03768          rotDragIndexes = new int32[numAtoms];
03769          msg->get(numAtoms, rotDragIndexes);
03770          if (numRotDrag)
03771          {
03772            delete [] rotDragParams;
03773            rotDragParams = new RotDragParams[numRotDrag];
03774            msg->get(numRotDrag*sizeof(RotDragParams), (char*)rotDragParams);
03775          }
03776       }
03777       
03778       //  Get the "constant" torque information, if it is active
03779       if (simParams->consTorqueOn) {
03780          msg->get(numConsTorque);
03781          delete [] consTorqueIndexes;
03782          consTorqueIndexes = new int32[numAtoms];
03783          msg->get(numAtoms, consTorqueIndexes);
03784          if (numConsTorque)
03785          {
03786            delete [] consTorqueParams;
03787            consTorqueParams = new ConsTorqueParams[numConsTorque];
03788            msg->get(numConsTorque*sizeof(ConsTorqueParams), (char*)consTorqueParams);
03789          }
03790       }
03791       
03792       // Get the constant force information, if it's active
03793       if (simParams->consForceOn)
03794       { msg->get(numConsForce);
03795         delete [] consForceIndexes;
03796         consForceIndexes = new int32[numAtoms];
03797         msg->get(numAtoms, consForceIndexes);
03798         if (numConsForce)
03799         { delete [] consForce;
03800           consForce = new Vector[numConsForce];
03801           msg->get(numConsForce*sizeof(Vector), (char*)consForce);
03802         }
03803       }
03804 
03805       //  Get the langevin parameters, if they are active
03806       if (simParams->langevinOn || simParams->tCoupleOn)
03807       {
03808         delete [] langevinParams;
03809         langevinParams = new Real[numAtoms];
03810 
03811         msg->get(numAtoms, langevinParams);
03812       }
03813 
03814       //  Get the fixed atoms, if they are active
03815       if (simParams->fixedAtomsOn)
03816       {
03817         delete [] fixedAtomFlags;
03818         fixedAtomFlags = new int32[numAtoms];
03819 
03820         msg->get(numFixedAtoms);
03821         msg->get(numAtoms, fixedAtomFlags);
03822         msg->get(numFixedRigidBonds);
03823       }
03824 
03825       if (simParams->excludeFromPressure) {
03826         exPressureAtomFlags = new int32[numAtoms];
03827         msg->get(numExPressureAtoms);
03828         msg->get(numAtoms, exPressureAtomFlags);
03829       }
03830 
03831 //fepb
03832       //receive fep atom info
03833       if (simParams->alchFepOn || simParams->lesOn || simParams->alchThermIntOn || simParams->pairInteractionOn) {
03834         delete [] fepAtomFlags;
03835         fepAtomFlags = new unsigned char[numAtoms];
03836 
03837         msg->get(numFepInitial);
03838         msg->get(numFepFinal);
03839         msg->get(numAtoms*sizeof(unsigned char), (char*)fepAtomFlags);
03840       }
03841 //fepe
03842 
03843       //  Now free the message 
03844       delete msg;
03845       
03846       //  analyze the data and find the status of each atom
03847       build_atom_status();
03848 
03849       //  Now build arrays of indexes into these arrays by atom
03850       #ifdef MEM_OPT_VERSION
03851       build_excl_check_signatures();
03852       #else
03853       build_lists_by_atom();      
03854       #endif
03855 
03856       #ifdef MEM_OPT_VERSION
03857       delEachAtomSigs();
03858       delChargeSpace();
03859       delMassSpace();
03860       delOtherEachAtomStructs();
03861       #endif
03862     }
03863     /*      END OF FUNCTION receive_Molecule    */
03864 
03865 
03866 #ifdef MEM_OPT_VERSION
03867     //Well, the exclusion check signatures could also done on PE0 and
03868     //sent to other processors through send_Molecule/receive_Molecule 
03869     //two procedures.
03870     void Molecule::build_excl_check_signatures(){
03871        exclChkSigPool = new ExclusionCheck[exclSigPoolSize];
03872        for(int i=0; i<exclSigPoolSize; i++){
03873            ExclusionSignature *sig = &exclSigPool[i];
03874            ExclusionCheck *sigChk = &exclChkSigPool[i];
03875            if(sig->fullExclCnt){
03876                if(!sig->modExclCnt){ //only having fullExclusion
03877                    sigChk->min = sig->fullOffset[0];
03878                    sigChk->max = sig->fullOffset[sig->fullExclCnt-1];
03879                }else{ //have both full and modified exclusion
03880                    int fullMin, fullMax, modMin, modMax;
03881                    
03882                    fullMin = sig->fullOffset[0];
03883                    fullMax = sig->fullOffset[sig->fullExclCnt-1];
03884                
03885                    modMin = sig->modOffset[0];
03886                    modMax = sig->modOffset[sig->modExclCnt-1];
03887                    
03888                    if(fullMin < modMin)
03889                        sigChk->min = fullMin;
03890                    else
03891                        sigChk->min = modMin;
03892                    if(fullMax < modMax)
03893                        sigChk->max = modMax;
03894                    else
03895                        sigChk->max = fullMax;
03896                }        
03897            }else{
03898                if(sig->modExclCnt){
03899                    sigChk->min = sig->modOffset[0];
03900                    sigChk->max = sig->modOffset[sig->modExclCnt-1];
03901                }else{ //both count are 0
03902                    if(CkMyPe()==0)
03903                        printf("Warning: an empty exclusion signature with index %d!\n", i);               
03904                    continue;
03905                }
03906            }           
03907 
03908            if((sigChk->max-sigChk->min) > simParams->maxExclusionFlags){
03909                printf("The distance of a exclusion check %d exceeds the value %d simulation parameter maxExclusionFlags specifies! Please increase the value\n", sigChk->max-sigChk->min, simParams->maxExclusionFlags);
03910                NAMD_die("Currently not supporting building exclusion check on the fly for memory optimized version!\n");
03911            }
03912 
03913            sigChk->flags = new char[sigChk->max-sigChk->min+1];
03914            memset(sigChk->flags, 0, sizeof(char)*(sigChk->max-sigChk->min+1));
03915            for(int j=0; j<sig->fullExclCnt; j++){
03916                int dist = sig->fullOffset[j] - sigChk->min;
03917                sigChk->flags[dist] = EXCHCK_FULL;
03918            }
03919            for(int j=0; j<sig->modExclCnt; j++){
03920                int dist = sig->modOffset[j] - sigChk->min;
03921                sigChk->flags[dist] = EXCHCK_MOD;
03922            }
03923        }
03924     }
03925 #endif
03926 
03927     /************************************************************************/
03928     /*                  */
03929     /*      FUNCTION build_lists_by_atom      */
03930     /*                  */
03931     /*  This function builds O(NumAtoms) arrays that store the bonds,   */
03932     /*  angles, dihedrals, and impropers, that each atom is involved in.    */
03933     /*  This is a space hog, but VERY fast.  This will certainly have to    */
03934     /*  change to make things scalable in memory, but for now, speed is the */
03935     /*  thing!                */
03936     /*                  */
03937     /************************************************************************/
03938 #ifdef MEM_OPT_VERSION
03939 
03940     //In memory optimized version, this function is only performed on PE0 since other nodes
03941     //may have less memory than the master node.
03942     //3 tasks is performed in the function: 1) cluster building 2) atom signatures changed due to
03943     // fixed atoms and alch atoms 3) exclusion signatures (in the format of offset) changed due to
03944     //hydrogen group/alch/fixed atoms.
03945     void Molecule::build_lists_by_atom()
03946        
03947     {
03948        register int i;      //  Loop counter
03949 
03950        register int numFixedAtoms = this->numFixedAtoms;
03951        // if we want forces on fixed atoms then just pretend
03952        // there are none for the purposes of this routine;
03953        if ( simParams->fixedAtomsForces ) numFixedAtoms = 0;
03954                      
03955        const int pair_self = 
03956          simParams->pairInteractionOn ? simParams->pairInteractionSelf : 0;
03957 
03958        //Cluster building is done when compressing the psf file
03959        
03960        //The safety check for bonds has been already done
03961        //in the period of generating compressed psf file
03962        DebugM(3, "Building the atom list for each signature.\n");
03963 
03964        //build up tuplesByAtom structure where each atom's signature may be updated
03965        //because of fixed atoms, atom signature pool will be changed.       
03966        //The method to update the current signature pool:
03967 
03968        //For each atom signature, go through those atoms that belong to this signature
03969        //and split the current signature correspondingly. 
03970        // Note: 1.The splitted signatures originating from the same atom signature will always
03971        // be different. But two splitted signatures may be same if they originate from two different
03972        // atom signatures. Currently, this redundancy is not handled for the sake of convenience
03973        // and simplicity. I believe practically this will not blow up memory since
03974        // the space for atom signatures will always remain very low (~10K) 
03975        //       2. This method is used under the assumption that chance of splitting a current
03976        // atom signature is small
03977 
03978        // There is another way to do this: (The following code is an implementation of this
03979        // 1. Finding all the tuples that are not needed to be considered either for 
03980        // fixed number or for pair_self
03981        // 2. Build atom signatures for them
03982        // 3. Add these new signatures to the existing atom signature pool
03983        // -- Chao Mei       
03984 
03985        numCalcBonds = numBonds;
03986        numCalcAngles = numAngles;
03987        numCalcDihedrals = numDihedrals;
03988        numCalcImpropers = numImpropers;
03989        numCalcCrossterms = numCrossterms;
03990        if(numFixedAtoms || pair_self){ //either condition holds                            
03991            vector<AtomSignature> newAtomSigPool;
03992            //we assume there are always one empty atom signature which is indexed
03993            //at 0 in the newAtomSigPool but with atomSigPoolSize in the final combined
03994            //atom signature pool
03995            AtomSignature emptyAtomSig;
03996            newAtomSigPool.push_back(emptyAtomSig);
03997 
03998            //begin to building new atom signatures for those atoms
03999            //affected by fixed atoms and fepAtomFlags. Besides,
04000            //the affected atoms' signature index should be updated too
04001            AtomSignature newAtomSig;           
04002            for(i=0; i<numAtoms; i++){ //"i" is the id of atom
04003                AtomSignature *curAtomSig = &atomSigPool[eachAtomSig[i]];
04004 
04005                //deal with fep atoms
04006                if(pair_self && fepAtomFlags[i]!=1){
04007                    //all tuples associated with this atom are gone                   
04008                    //update the number of calculated tuples
04009                    numCalcBonds -= curAtomSig->bondCnt;
04010                    numCalcAngles -= curAtomSig->angleCnt;
04011                    numCalcDihedrals -= curAtomSig->dihedralCnt;
04012                    numCalcImpropers -= curAtomSig->improperCnt;
04013                    numCalcCrossterms -= curAtomSig->crosstermCnt;
04014 
04015                    //update this atom's signature
04016                    eachAtomSig[i] = (Index)atomSigPoolSize;                   
04017                    
04018                    continue;
04019                }
04020 
04021                //deal with fixed atoms
04022                int sigChanged = 0;
04023                if(numFixedAtoms && fixedAtomFlags[i]){
04024                    //examine all the tuples associated with this atom                   
04025                    newAtomSig = *curAtomSig;
04026 
04027                    //1. bonds
04028                    int tupleCnt = curAtomSig->bondCnt;
04029                    TupleSignature *tupleSigs = curAtomSig->bondSigs;
04030                    for(int j=0; j<tupleCnt; j++){
04031                        int atom2 = i+tupleSigs[j].offset[0];
04032                        if(fixedAtomFlags[atom2]){ // find a fixed tuple                           
04033                            newAtomSig.bondSigs[j].setEmpty();
04034                            numCalcBonds--;
04035                            sigChanged = 1;
04036                        }
04037                    }
04038 
04039                    //2. angles
04040                    tupleCnt = curAtomSig->angleCnt;
04041                    tupleSigs = curAtomSig->angleSigs;
04042                    for(int j=0; j<tupleCnt; j++){
04043                        int atom2 = i+tupleSigs[j].offset[0];
04044                        int atom3 = i+tupleSigs[j].offset[1];
04045                        if(fixedAtomFlags[atom2] && fixedAtomFlags[atom3]){ // find a fixed tuple                           
04046                            newAtomSig.angleSigs[j].setEmpty();
04047                            numCalcAngles--;
04048                            sigChanged = 1;
04049                        }
04050                    }
04051 
04052                    //3. dihedrals
04053                    tupleCnt = curAtomSig->dihedralCnt;
04054                    tupleSigs = curAtomSig->dihedralSigs;
04055                    for(int j=0; j<tupleCnt; j++){
04056                        int atom2 = i+tupleSigs[j].offset[0];
04057                        int atom3 = i+tupleSigs[j].offset[1];
04058                        int atom4 = i+tupleSigs[j].offset[2];
04059                        if(fixedAtomFlags[atom2] && fixedAtomFlags[atom3]
04060                           && fixedAtomFlags[atom4]){ // find a fixed tuple                           
04061                            newAtomSig.dihedralSigs[j].setEmpty();
04062                            numCalcDihedrals--;
04063                            sigChanged = 1;
04064                        }
04065                    }
04066 
04067                    //4. impropers
04068                    tupleCnt = curAtomSig->improperCnt;
04069                    tupleSigs = curAtomSig->improperSigs;
04070                    for(int j=0; j<tupleCnt; j++){
04071                        int atom2 = i+tupleSigs[j].offset[0];
04072                        int atom3 = i+tupleSigs[j].offset[1];
04073                        int atom4 = i+tupleSigs[j].offset[2];
04074                        if(fixedAtomFlags[atom2] && fixedAtomFlags[atom3]
04075                           && fixedAtomFlags[atom4]){ // find a fixed tuple                           
04076                            newAtomSig.improperSigs[j].setEmpty();
04077                            numCalcImpropers--;
04078                            sigChanged = 1;
04079                        }
04080                    }
04081 
04082                    //5. crossterms
04083                    tupleCnt = curAtomSig->crosstermCnt;
04084                    tupleSigs = curAtomSig->crosstermSigs;
04085                    for(int j=0; j<tupleCnt; j++){
04086                        int atom2 = i+tupleSigs[j].offset[0];
04087                        int atom3 = i+tupleSigs[j].offset[1];
04088                        int atom4 = i+tupleSigs[j].offset[2];
04089                        int atom5 = i+tupleSigs[j].offset[3];
04090                        int atom6 = i+tupleSigs[j].offset[4];
04091                        int atom7 = i+tupleSigs[j].offset[5];
04092                        int atom8 = i+tupleSigs[j].offset[6];
04093 
04094                        if(fixedAtomFlags[atom2] && fixedAtomFlags[atom3]
04095                           && fixedAtomFlags[atom4] && fixedAtomFlags[atom5]
04096                           && fixedAtomFlags[atom6] && fixedAtomFlags[atom7]
04097                           && fixedAtomFlags[atom8]){ // find a fixed tuple                           
04098                            newAtomSig.crosstermSigs[j].setEmpty();
04099                            numCalcCrossterms--;
04100                            sigChanged = 1;
04101                        }
04102                    }
04103 
04104                    //update the newAtomSig to remove all empty tuple signatures
04105                    if(sigChanged){
04106                        newAtomSig.removeEmptyTupleSigs();
04107 
04108                        //add the new atom signature to the newAtomSigPool if it is not in it
04109                        int poolIndex = lookupCstPool(newAtomSigPool, newAtomSig);
04110                        if(poolIndex==-1){
04111                            poolIndex = newAtomSigPool.size();
04112                            newAtomSigPool.push_back(newAtomSig);
04113                        }
04114                        eachAtomSig[i] = (Index)(poolIndex+atomSigPoolSize);
04115                    }
04116                }
04117            }
04118            //update the atom signature pool from the new one                      
04119            AtomSignature *tmpAtomSigPool = new AtomSignature[atomSigPoolSize+newAtomSigPool.size()];
04120            for(i=0; i<atomSigPoolSize; i++)
04121                tmpAtomSigPool[i] = atomSigPool[i];           
04122            for(i=0; i<newAtomSigPool.size(); i++)
04123                tmpAtomSigPool[atomSigPoolSize+i] = newAtomSigPool[i];
04124 
04125            delete [] atomSigPool;
04126            atomSigPoolSize += newAtomSigPool.size();
04127            atomSigPool = tmpAtomSigPool;
04128            newAtomSigPool.clear();
04129 
04130            if(CkMyPe()==0){
04131                printf("Current number of atom signatures: %d\n", atomSigPoolSize);
04132            }
04133        }
04134           
04135        DebugM(3,"Building exclusion data.\n");
04136 
04137 
04138        //  Build the arrays of exclusions for each atom       
04139        //inside this function call, unnecessary exclusions are eliminated in the
04140        //case of fep atoms and fixed atoms
04141        if(CkMyPe()==0){
04142            CkPrintf("Current exclusion signatures: %d\n", exclSigPoolSize);
04143        }
04144        build_exclusions();       
04145        
04146        if(CkMyPe()==0){
04147            CkPrintf("After stripping, exclusion signatures: %d\n", exclSigPoolSize);
04148        }
04149        
04150        //Finally build the exclChkSigPool of type ExclusionCheck which is separated
04151        //from this function for saving memory
04152        //build_excl_check_signature();
04153     }
04154 #else
04155     void Molecule::build_lists_by_atom()
04156        
04157     {
04158        register int i;      //  Loop counter
04159 
04160        register int numFixedAtoms = this->numFixedAtoms;
04161        // if we want forces on fixed atoms then just pretend
04162        // there are none for the purposes of this routine;
04163        if ( simParams->fixedAtomsForces ) numFixedAtoms = 0;
04164 
04165 //fepb
04166 //     int numFepInitial = this->numFepInitial;
04167 //     int numFepFinal = this->numFepFinal;
04168 //fepe
04169        tmpArena = new ObjectArena<int32>;
04170        bondsWithAtom = new int32 *[numAtoms];
04171        cluster = new int32 [numAtoms];
04172        clusterSize = new int32 [numAtoms];
04173 
04174        bondsByAtom = new int32 *[numAtoms];
04175        anglesByAtom = new int32 *[numAtoms];
04176        dihedralsByAtom = new int32 *[numAtoms];
04177        impropersByAtom = new int32 *[numAtoms];
04178        crosstermsByAtom = new int32 *[numAtoms];
04179 
04180        exclusionsByAtom = new int32 *[numAtoms];
04181        fullExclusionsByAtom = new int32 *[numAtoms];
04182        modExclusionsByAtom = new int32 *[numAtoms];
04183 
04184        int32 *byAtomSize = new int32[numAtoms];
04185 
04186        const int pair_self = 
04187          simParams->pairInteractionOn ? simParams->pairInteractionSelf : 0;
04188 
04189        DebugM(3,"Building bond lists.\n");
04190     
04191        //  Build the bond lists
04192        for (i=0; i<numAtoms; i++)
04193        {
04194          byAtomSize[i] = 0;
04195        }
04196        for (i=0; i<numRealBonds; i++)
04197        {
04198          byAtomSize[bonds[i].atom1]++;
04199          byAtomSize[bonds[i].atom2]++;
04200        }
04201        for (i=0; i<numAtoms; i++)
04202        {
04203          bondsWithAtom[i] = tmpArena->getNewArray(byAtomSize[i]+1);
04204          bondsWithAtom[i][byAtomSize[i]] = -1;
04205          byAtomSize[i] = 0;
04206        }
04207        for (i=0; i<numRealBonds; i++)
04208        {
04209          int a1 = bonds[i].atom1;
04210          int a2 = bonds[i].atom2;
04211          bondsWithAtom[a1][byAtomSize[a1]++] = i;
04212          bondsWithAtom[a2][byAtomSize[a2]++] = i;
04213        }
04214 
04215        //  Build cluster information (contiguous clusters)
04216        for (i=0; i<numAtoms; i++) {
04217          cluster[i] = i;
04218        }
04219        for (i=0; i<numAtoms; i++) {
04220          int ci = i;
04221          while ( cluster[ci] != ci ) ci = cluster[ci];
04222          for ( int32 *b = bondsWithAtom[i]; *b != -1; ++b ) {
04223            int a = bonds[*b].atom1;
04224            if ( a == i ) a = bonds[*b].atom2;
04225            if ( a > i ) {
04226              int ca = a;
04227              while ( cluster[ca] != ca ) ca = cluster[ca];
04228              if ( ca > ci ) cluster[ca] = cluster[ci];
04229              else cluster[ci] = cluster[ca];
04230            }
04231          }
04232        }
04233        while ( 1 ) {
04234          int allok = 1;
04235          for (i=0; i<numAtoms; i++) {
04236            int ci = cluster[i];
04237            if ( cluster[ci] != ci ) {
04238              allok = 0;
04239              cluster[i] = cluster[ci];
04240            }
04241          }
04242          if ( allok ) break;
04243        }
04244        
04245        for (i=0; i<numAtoms; i++) {
04246          clusterSize[i] = 0;
04247        }       
04248        for (i=0; i<numAtoms; i++) {           
04249          clusterSize[cluster[i]] += 1;
04250        }
04251 
04252 /*
04253        //Getting number of clusters for debugging
04254        int numClusters=0;
04255        for(int i=0; i<numAtoms; i++){
04256            if(clusterSize[i]!=0) numClusters++;
04257        }
04258        printf("Num of clusters: %d\n", numClusters);
04259 */
04260 
04261        //  Build the bond lists
04262        for (i=0; i<numAtoms; i++)
04263        {
04264          byAtomSize[i] = 0;
04265        }
04266        numCalcBonds = 0;
04267        for (i=0; i<numBonds; i++)
04268        {
04269          if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
04270                             && fixedAtomFlags[bonds[i].atom2] ) continue;
04271    
04272          if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
04273          byAtomSize[bonds[i].atom1]++;
04274          numCalcBonds++;
04275        }
04276        for (i=0; i<numAtoms; i++)
04277        {
04278          bondsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
04279          bondsByAtom[i][byAtomSize[i]] = -1;
04280          byAtomSize[i] = 0;
04281        }
04282        for (i=0; i<numBonds; i++)
04283        {
04284          if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
04285                             && fixedAtomFlags[bonds[i].atom2] ) continue;
04286          if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
04287          int a1 = bonds[i].atom1;
04288          bondsByAtom[a1][byAtomSize[a1]++] = i;
04289        }
04290        for (i=0; i<numBonds; ++i) {
04291          int a1 = bonds[i].atom1;
04292          int a2 = bonds[i].atom2;
04293          int j;
04294          if ( a1 == a2 ) {
04295            char buff[512];
04296            sprintf(buff,"Atom %d is bonded to itself", a1+1);
04297            NAMD_die(buff);
04298          }
04299          for ( j = 0; j < byAtomSize[a1]; ++j ) {
04300            int b = bondsByAtom[a1][j];
04301            int ba1 = bonds[b].atom1;
04302            int ba2 = bonds[b].atom2;
04303            if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
04304              char buff[512];
04305              sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
04306              NAMD_die(buff);
04307            }
04308          }
04309          for ( j = 0; j < byAtomSize[a2]; ++j ) {
04310            int b = bondsByAtom[a2][j];
04311            int ba1 = bonds[b].atom1;
04312            int ba2 = bonds[b].atom2;
04313            if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
04314              char buff[512];
04315              sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
04316              NAMD_die(buff);
04317            }
04318          }
04319        }
04320        
04321        DebugM(3,"Building angle lists.\n");
04322     
04323        //  Build the angle lists
04324        for (i=0; i<numAtoms; i++)
04325        {
04326          byAtomSize[i] = 0;
04327        }
04328        numCalcAngles = 0;
04329        for (i=0; i<numAngles; i++)
04330        {
04331          if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
04332                             && fixedAtomFlags[angles[i].atom2]
04333                             && fixedAtomFlags[angles[i].atom3] ) continue;
04334          if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
04335          byAtomSize[angles[i].atom1]++;
04336          numCalcAngles++;
04337        }
04338        for (i=0; i<numAtoms; i++)
04339        {
04340          anglesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
04341          anglesByAtom[i][byAtomSize[i]] = -1;
04342          byAtomSize[i] = 0;
04343        }
04344        for (i=0; i<numAngles; i++)
04345        {
04346          if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
04347          if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
04348                             && fixedAtomFlags[angles[i].atom2]
04349                             && fixedAtomFlags[angles[i].atom3] ) continue;
04350          int a1 = angles[i].atom1;
04351          anglesByAtom[a1][byAtomSize[a1]++] = i;
04352        }
04353        
04354        DebugM(3,"Building improper lists.\n");
04355     
04356        //  Build the improper lists
04357        for (i=0; i<numAtoms; i++)
04358        {
04359          byAtomSize[i] = 0;
04360        }
04361        numCalcImpropers = 0;
04362        for (i=0; i<numImpropers; i++)
04363        {
04364          if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
04365                             && fixedAtomFlags[impropers[i].atom2]
04366                             && fixedAtomFlags[impropers[i].atom3]
04367                             && fixedAtomFlags[impropers[i].atom4] ) continue;
04368          if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
04369          byAtomSize[impropers[i].atom1]++;
04370          numCalcImpropers++;
04371        }
04372        for (i=0; i<numAtoms; i++)
04373        {
04374          impropersByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
04375          impropersByAtom[i][byAtomSize[i]] = -1;
04376          byAtomSize[i] = 0;
04377        }
04378        for (i=0; i<numImpropers; i++)
04379        {
04380          if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
04381                             && fixedAtomFlags[impropers[i].atom2]
04382                             && fixedAtomFlags[impropers[i].atom3]
04383                             && fixedAtomFlags[impropers[i].atom4] ) continue;
04384          if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
04385          int a1 = impropers[i].atom1;
04386          impropersByAtom[a1][byAtomSize[a1]++] = i;
04387        }
04388        
04389        DebugM(3,"Building dihedral lists.\n");
04390     
04391        //  Build the dihedral lists
04392        for (i=0; i<numAtoms; i++)
04393        {
04394          byAtomSize[i] = 0;
04395        }
04396        numCalcDihedrals = 0;
04397        for (i=0; i<numDihedrals; i++)
04398        {
04399          if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
04400                             && fixedAtomFlags[dihedrals[i].atom2]
04401                             && fixedAtomFlags[dihedrals[i].atom3]
04402                             && fixedAtomFlags[dihedrals[i].atom4] ) continue;
04403          if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
04404          byAtomSize[dihedrals[i].atom1]++;
04405          numCalcDihedrals++;
04406        }
04407        for (i=0; i<numAtoms; i++)
04408        {
04409          dihedralsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
04410          dihedralsByAtom[i][byAtomSize[i]] = -1;
04411          byAtomSize[i] = 0;
04412        }
04413        for (i=0; i<numDihedrals; i++)
04414        {
04415          if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
04416                             && fixedAtomFlags[dihedrals[i].atom2]
04417                             && fixedAtomFlags[dihedrals[i].atom3]
04418                             && fixedAtomFlags[dihedrals[i].atom4] ) continue;
04419          if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
04420          int a1 = dihedrals[i].atom1;
04421          dihedralsByAtom[a1][byAtomSize[a1]++] = i;
04422        }
04423     
04424        DebugM(3,"Building crossterm lists.\n");
04425     
04426        //  Build the crossterm lists
04427        for (i=0; i<numAtoms; i++)
04428        {
04429          byAtomSize[i] = 0;
04430        }
04431        numCalcCrossterms = 0;
04432        for (i=0; i<numCrossterms; i++)
04433        {
04434          if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
04435                             && fixedAtomFlags[crossterms[i].atom2]
04436                             && fixedAtomFlags[crossterms[i].atom3]
04437                             && fixedAtomFlags[crossterms[i].atom4]
04438                             && fixedAtomFlags[crossterms[i].atom5]
04439                             && fixedAtomFlags[crossterms[i].atom6]
04440                             && fixedAtomFlags[crossterms[i].atom7]
04441                             && fixedAtomFlags[crossterms[i].atom8] ) continue;
04442          if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
04443          byAtomSize[crossterms[i].atom1]++;
04444          numCalcCrossterms++;
04445        }
04446        for (i=0; i<numAtoms; i++)
04447        {
04448          crosstermsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
04449          crosstermsByAtom[i][byAtomSize[i]] = -1;
04450          byAtomSize[i] = 0;
04451        }
04452        for (i=0; i<numCrossterms; i++)
04453        {
04454          if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
04455                             && fixedAtomFlags[crossterms[i].atom2]
04456                             && fixedAtomFlags[crossterms[i].atom3]
04457                             && fixedAtomFlags[crossterms[i].atom4]
04458                             && fixedAtomFlags[crossterms[i].atom5]
04459                             && fixedAtomFlags[crossterms[i].atom6]
04460                             && fixedAtomFlags[crossterms[i].atom7]
04461                             && fixedAtomFlags[crossterms[i].atom8] ) continue;
04462          if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
04463          int a1 = crossterms[i].atom1;
04464          crosstermsByAtom[a1][byAtomSize[a1]++] = i;
04465        }
04466 
04467        // DRUDE: init lphostIndexes array
04468        DebugM(3,"Initializing lone pair host index array.\n");
04469        lphostIndexes = new int32[numAtoms];
04470        for (i = 0;  i < numAtoms;  i++) {
04471          lphostIndexes[i] = -1;
04472        }
04473        for (i = 0;  i < numLphosts;  i++) {
04474          int32 index = lphosts[i].atom1;
04475          lphostIndexes[index] = i;
04476        }
04477        // DRUDE
04478     
04479        DebugM(3,"Building exclusion data.\n");
04480     
04481        //  Build the arrays of exclusions for each atom
04482        build_exclusions();
04483 
04484        //  Remove temporary structures
04485        delete [] bondsWithAtom;  bondsWithAtom = 0;
04486        delete tmpArena;  tmpArena = 0;
04487 
04488        if (exclusions != NULL)
04489       delete [] exclusions;
04490 
04491        // 1-4 exclusions which are also fully excluded were eliminated by hash table
04492        numTotalExclusions = exclusionSet.size();
04493        exclusions = new Exclusion[numTotalExclusions];
04494        UniqueSetIter<Exclusion> exclIter(exclusionSet);
04495        for ( exclIter=exclIter.begin(),i=0; exclIter != exclIter.end(); exclIter++,i++ )
04496        {
04497          exclusions[i] = *exclIter;
04498        }
04499        // Free exclusionSet storage
04500        // exclusionSet.clear(1);
04501        exclusionSet.clear();
04502 
04503        DebugM(3,"Building exclusion lists.\n");
04504     
04505        for (i=0; i<numAtoms; i++)
04506        {
04507          byAtomSize[i] = 0;
04508        }
04509        numCalcExclusions = 0;
04510        for (i=0; i<numTotalExclusions; i++)
04511        {
04512          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
04513                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
04514          byAtomSize[exclusions[i].atom1]++;
04515          numCalcExclusions++;
04516        }
04517        for (i=0; i<numAtoms; i++)
04518        {
04519          exclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
04520          exclusionsByAtom[i][byAtomSize[i]] = -1;
04521          byAtomSize[i] = 0;
04522        }
04523        for (i=0; i<numTotalExclusions; i++)
04524        {
04525          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
04526                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
04527          int a1 = exclusions[i].atom1;
04528          exclusionsByAtom[a1][byAtomSize[a1]++] = i;
04529        }
04530 
04531        int32 *byAtomSize2 = new int32[numAtoms];
04532 
04533        for (i=0; i<numAtoms; i++)
04534        {
04535          byAtomSize[i] = 0;
04536          byAtomSize2[i] = 0;
04537        }
04538 
04539        for (i=0; i<numTotalExclusions; i++)
04540        {
04541          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
04542                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
04543          if ( exclusions[i].modified ) {
04544            byAtomSize2[exclusions[i].atom1]++;
04545            byAtomSize2[exclusions[i].atom2]++;
04546          } else {
04547            byAtomSize[exclusions[i].atom1]++;
04548            byAtomSize[exclusions[i].atom2]++;
04549          }
04550        }
04551 
04552        for (i=0; i<numAtoms; i++)
04553        {
04554          fullExclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
04555          fullExclusionsByAtom[i][0] = 0;
04556          modExclusionsByAtom[i] = arena->getNewArray(byAtomSize2[i]+1);
04557          modExclusionsByAtom[i][0] = 0;
04558        }
04559 
04560        for (i=0; i<numTotalExclusions; i++)
04561        {
04562          int a1 = exclusions[i].atom1;
04563          int a2 = exclusions[i].atom2;
04564          if ( numFixedAtoms && fixedAtomFlags[a1]
04565                             && fixedAtomFlags[a2] ) continue;
04566          int32 *l1, *l2;
04567          if ( exclusions[i].modified ) {
04568            l1 = modExclusionsByAtom[a1];
04569            l2 = modExclusionsByAtom[a2];
04570          } else {
04571            l1 = fullExclusionsByAtom[a1];
04572            l2 = fullExclusionsByAtom[a2];
04573          }
04574          l1[++(*l1)] = a2;
04575          l2[++(*l2)] = a1;
04576        }
04577 
04578        delete [] byAtomSize;  byAtomSize = 0;
04579        delete [] byAtomSize2;  byAtomSize2 = 0;
04580 
04581 
04582        //  Allocate an array to hold the exclusions for each atom
04583        all_exclusions = new ExclusionCheck[numAtoms];
04584 
04585        for (i=0; i<numAtoms; i++)
04586        {
04587          all_exclusions[i].min = numAtoms;
04588          all_exclusions[i].max = -1;
04589        }
04590        for (i=0; i<numTotalExclusions; i++)
04591        {
04592          // first atom should alway have lower number!
04593          int a1 = exclusions[i].atom1;
04594          int a2 = exclusions[i].atom2;
04595          if ( numFixedAtoms && fixedAtomFlags[a1]
04596                             && fixedAtomFlags[a2] ) continue;
04597          if ( all_exclusions[a1].min > a2 ) all_exclusions[a1].min = a2;
04598          if ( all_exclusions[a2].min > a1 ) all_exclusions[a2].min = a1;
04599          if ( a2 > all_exclusions[a1].max ) all_exclusions[a1].max = a2;
04600          if ( a1 > all_exclusions[a2].max ) all_exclusions[a2].max = a1;
04601        }
04602        int exclmem = 0;
04603        int maxExclusionFlags = simParams->maxExclusionFlags;
04604        for (i=0; i<numAtoms; i++)
04605        {
04606          int s = all_exclusions[i].max - all_exclusions[i].min + 1;
04607          if ( all_exclusions[i].max != -1 ) {
04608            if ( s < maxExclusionFlags ) {
04609              char *f = all_exclusions[i].flags = exclArena->getNewArray(s);
04610              for ( int k=0; k<s; ++k ) f[k] = 0;
04611              exclmem += s;
04612            } else {
04613              all_exclusions[i].flags = 0;  // need to build on the fly
04614            }
04615          } else {
04616            all_exclusions[i].flags = (char*)-1; // should never dereference
04617          }
04618        }
04619        if ( 0 ) {
04620          iout << iINFO << numTotalExclusions << " exclusions consume "
04621             << exclmem << " bytes.\n" << endi;
04622        }
04623        for (i=0; i<numTotalExclusions; i++)
04624        {
04625          int a1 = exclusions[i].atom1;
04626          int a2 = exclusions[i].atom2;
04627          if ( numFixedAtoms && fixedAtomFlags[a1]
04628                             && fixedAtomFlags[a2] ) continue;
04629          if ( exclusions[i].modified ) {
04630            if ( all_exclusions[a1].flags )
04631              all_exclusions[a1].flags[a2-all_exclusions[a1].min] = EXCHCK_MOD;
04632            if ( all_exclusions[a2].flags )
04633              all_exclusions[a2].flags[a1-all_exclusions[a2].min] = EXCHCK_MOD;
04634          } else {
04635            if ( all_exclusions[a1].flags )
04636              all_exclusions[a1].flags[a2-all_exclusions[a1].min] = EXCHCK_FULL;
04637            if ( all_exclusions[a2].flags )
04638              all_exclusions[a2].flags[a1-all_exclusions[a2].min] = EXCHCK_FULL;
04639          }
04640        }
04641 
04642     }
04643 #endif
04644 
04645     /*    END OF FUNCTION build_lists_by_atom    */
04646 
04647     /****************************************************************/
04648     /*                */
04649     /*      FUNCTION build_exclusions    */
04650     /*                */
04651     /*  This function builds a list of all the exlcusions       */
04652     /*  atoms.  These lists include explicit exclusions as well as  */
04653     /*  exclusions that are calculated based on the bonded structure*/
04654     /*  and the exclusion flag.  For each pair of atoms that are    */
04655     /*  excluded, the larger of the 2 atom indexes is stored in the */
04656     /*  array of the smaller index.  All the arrays are not sorted. */
04657     /*  Then to determine if two atoms have an exclusion, a linear  */
04658     /*  search is done on the array of the atom with the smaller    */
04659     /*  index for the larger index.          */
04660     /*  If the exclusion policy is set to scaled1-4, there are  */
04661     /*  actually two lists built.  One contains the pairs of atoms  */
04662     /*  that are to be exlcuded (i.e., explicit exclusions, 1-2,    */
04663     /*  and 1-3 interactions) and the other contains just the 1-4   */
04664     /*  interactions, since they will need to be determined   */
04665     /*  independantly of the other exclusions.      */
04666     /*                */
04667     /****************************************************************/
04668 
04669     void Molecule::build_exclusions()
04670     {
04671 #ifdef MEM_OPT_VERSION
04672         numCalcExclusions = numTotalExclusions*2;        
04673 
04674         //stripHGroupExcl() is not needed because the function purpose 
04675         //is to reduce the memory usage of exclusions. Since now the exclusion
04676         //signature is used, such function is not critical now
04677         //What about the following function? Need to be confirmed --Chao Mei
04678         
04679         ExclusionSettings exclude_flag;    //  Exclusion policy
04680                 exclude_flag = simParams->exclude;
04681                 int stripHGroupExclFlag = (simParams->splitPatch == SPLIT_PATCH_HYDROGEN);
04682         if (!simParams->amberOn || !simParams->readExclusions){                 
04683                 switch (exclude_flag)
04684                 {
04685                  case NONE:                
04686                  case ONETWO:           
04687                    break;
04688                   case ONETHREE:                                
04689                   case ONEFOUR:                 
04690                   case SCALED14:            
04691                         if ( stripHGroupExclFlag ) stripHGroupExcl();
04692                     break;
04693                 }
04694                 }
04695                 else if (stripHGroupExclFlag && exclude_flag!=NONE && exclude_flag!=ONETWO)
04696                         stripHGroupExcl();               
04697 
04698         stripFepFixedExcl();
04699         return;
04700 #else
04701       register int i;          //  Loop counter
04702       ExclusionSettings exclude_flag;    //  Exclusion policy
04703 
04704       exclude_flag = simParams->exclude;
04705       int stripHGroupExclFlag = (simParams->splitPatch == SPLIT_PATCH_HYDROGEN);
04706 
04707       //  Go through the explicit exclusions and add them to the arrays
04708       for (i=0; i<numExclusions; i++)
04709       {
04710         exclusionSet.add(exclusions[i]);
04711       }
04712 
04713       // If this is AMBER force field, and readExclusions is TRUE,
04714       // then all the exclusions were read from parm file, and we
04715       // shouldn't generate any of them.
04716       if (!simParams->amberOn || !simParams->readExclusions)
04717       { //  Now calculate the bonded exlcusions based on the exclusion policy
04718         switch (exclude_flag)
04719         {
04720          case NONE:
04721            break;
04722          case ONETWO:
04723            build12excl();
04724            break;
04725           case ONETHREE:
04726             build12excl();
04727             build13excl();
04728             if ( stripHGroupExclFlag ) stripHGroupExcl();
04729             break;
04730           case ONEFOUR:
04731             build12excl();
04732             build13excl();
04733             build14excl(0);
04734             if ( stripHGroupExclFlag ) stripHGroupExcl();
04735             break;
04736           case SCALED14:
04737             build12excl();
04738             build13excl();
04739             build14excl(1);
04740             if ( stripHGroupExclFlag ) stripHGroupExcl();
04741             break;
04742         }
04743       }
04744       else if (stripHGroupExclFlag && exclude_flag!=NONE && exclude_flag!=ONETWO)
04745         stripHGroupExcl();
04746 
04747       stripFepExcl();
04748 
04749       // DRUDE
04750       if (is_drude_psf) build_inherited_excl();
04751 #endif
04752     }
04753     /*      END OF FUNCTION build_exclusions    */
04754 
04755 
04756     // DRUDE: extend exclusions for Drude and LP
04757     // Drude and LP particles "inherit" exclusions from their parents
04758     void Molecule::build_inherited_excl(void) {
04759 #ifdef MEM_OPT_VERSION
04760       NAMD_die("Drude and LP particles not supported in memopt version.");
04761 #else
04762       ExclusionSettings exclude_flag = simParams->exclude;
04763       int32 *bond1, *bond2, *bond3;
04764       int32 i, mid1, mid2;
04765 
04766       if (exclude_flag == ONEFOUR || exclude_flag == SCALED14) {
04767         NAMD_die("DRUDE MODEL SUPPORTS ONLY UP TO 1-3 EXCLUSION POLICY");
04768       }
04769 
04770       if (exclude_flag == NONE) return;
04771 
04772       // validate that each Drude or lone pair particle
04773       // has a unique parent that is a heavy atom
04774       for (i = 0;  i < numAtoms;  i++) {
04775 
04776         if (is_drude(i) || is_lp(i)) {
04777           // find parent (heavy) atom of particle i
04778           bond1 = bondsWithAtom[i];
04779 
04780           if (-1 == *bond1) {  // i must have one bond
04781             char err_msg[512];
04782             const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
04783             sprintf(err_msg, "FOUND ISOLATED %s PARTICLE %d", idescrip, i+1);
04784             NAMD_die(err_msg);
04785           }
04786           if (-1 != *(bond1+1)) {  // and only one bond
04787             char err_msg[512];
04788             const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
04789             sprintf(err_msg, "FOUND MULTIPLY LINKED %s PARTICLE %d",
04790                 idescrip, i+1);
04791             NAMD_die(err_msg);
04792           }
04793 
04794           // mid1 is parent of particle i
04795           mid1 = bonds[*bond1].atom1;
04796           if (mid1 == i) mid1 = bonds[*bond1].atom2;
04797 
04798           // make sure that mid1 is a heavy atom
04799           if (is_drude(mid1) || is_lp(mid1) || is_hydrogen(mid1)) {
04800             char err_msg[512];
04801             const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
04802             sprintf(err_msg, "PARENT ATOM %d of %s PARTICLE %d "
04803                 "IS NOT HEAVY ATOM", mid1+1, idescrip, i+1);
04804             NAMD_die(err_msg);
04805           }
04806 
04807           // follow build14excl() code
04808           bond2 = bondsWithAtom[mid1];
04809 
04810           // loop through all the bonds connected to atom mid1
04811           while (*bond2 != -1) {
04812             if (bonds[*bond2].atom1 == mid1) {
04813               mid2 = bonds[*bond2].atom2;
04814             }
04815             else {
04816               mid2 = bonds[*bond2].atom1;
04817             }
04818 
04819             // Make sure that we don't double back to where we started from.
04820             // Doing so causes strange behavior.
04821             if (mid2 == i) {
04822               bond2++;
04823               continue;
04824             }
04825 
04826             if (exclude_flag == ONETWO) {
04827               // add (i,mid2) as an exclusion
04828               if (i < mid2) {
04829                 exclusionSet.add(Exclusion(i, mid2));
04830               }
04831               else {
04832                 exclusionSet.add(Exclusion(mid2, i));
04833               }
04834             }
04835             else {  // exclude_flag == ONETHREE
04836 
04837               bond3 = bondsWithAtom[mid2];
04838 
04839               // loop through all the bonds connected to mid2
04840               while (*bond3 != -1) {
04841 
04842                 if (bonds[*bond3].atom1 == mid2) {
04843                   // Make sure we don't double back to where we started.
04844                   // Doing so causes strange behavior.
04845                   if (bonds[*bond3].atom2 != mid1) {
04846                     if (i < bonds[*bond3].atom2) {
04847                       exclusionSet.add(Exclusion(i, bonds[*bond3].atom2));
04848                     }
04849                     else if (bonds[*bond3].atom2 < i) {
04850                       exclusionSet.add(Exclusion(bonds[*bond3].atom2, i));
04851                     }
04852                   }
04853                 }
04854                 else {
04855                   // Make sure we don't double back to where we started.
04856                   // Doing so causes strange behavior.
04857                   if (bonds[*bond3].atom1 != mid1) {
04858                     if (i < bonds[*bond3].atom1) {
04859                       exclusionSet.add(Exclusion(i, bonds[*bond3].atom1));
04860                     }
04861                     else if (bonds[*bond3].atom1 < i) {
04862                       exclusionSet.add(Exclusion(bonds[*bond3].atom1, i));
04863                     }
04864                   }
04865                 }
04866 
04867                 ++bond3;
04868               } // while bond3
04869             } // else ONETHREE
04870            
04871             ++bond2;
04872           } // while bond2
04873 
04874         } // if i is Drude or LP
04875 
04876       } // for i
04877 
04878 #endif
04879     } 
04880     // DRUDE
04881 
04882 
04883     /************************************************************************/
04884     /*                  */
04885     /*      FUNCTION build12excl        */
04886     /*                  */
04887     /************************************************************************/
04888 
04889     void Molecule::build12excl(void)
04890        
04891     {
04892 #ifndef MEM_OPT_VERSION
04893        int32 *current_val;  //  Current value to check
04894        register int i;    //  Loop counter to loop through all atoms
04895        
04896        //  Loop through all the atoms marking the bonded interactions for each one
04897        for (i=0; i<numAtoms; i++)
04898        {
04899       current_val = bondsWithAtom[i];
04900        
04901       //  Loop through all the bonds for this atom
04902       while (*current_val != -1)
04903       {
04904          if (bonds[*current_val].atom1 == i)
04905          {
04906       if (i<bonds[*current_val].atom2)
04907       {
04908          exclusionSet.add(Exclusion(i,bonds[*current_val].atom2));
04909       }
04910          }
04911          else
04912          {
04913       if (i<bonds[*current_val].atom1)
04914       {
04915          exclusionSet.add(Exclusion(i,bonds[*current_val].atom1));
04916       }
04917          }
04918     
04919          ++current_val;
04920       }
04921        }
04922 #endif
04923     }
04924     /*      END OF FUNCTION build12excl      */
04925 
04926     /************************************************************************/
04927     /*                  */
04928     /*      FUNCTION build13excl        */
04929     /*                  */
04930     /************************************************************************/
04931 
04932     void Molecule::build13excl(void)
04933        
04934     {
04935 #ifndef MEM_OPT_VERSION
04936        int32 *bond1, *bond2;  //  The two bonds being checked
04937        int middle_atom;  //  Common third atom
04938        register int i;    //  Loop counter to loop through all atoms
04939        
04940        //  Loop through all the atoms looking at the bonded connections
04941        //  for each one
04942        for (i=0; i<numAtoms; i++)
04943        {
04944        bond1 = bondsWithAtom[i];
04945        
04946        //  Loop through all the bonds directly connect to atom i
04947        while (*bond1 != -1)
04948        {
04949         if (bonds[*bond1].atom1 == i)
04950         {
04951           middle_atom=bonds[*bond1].atom2;
04952         }
04953         else
04954         {
04955           middle_atom=bonds[*bond1].atom1;
04956         }
04957 
04958         bond2 = bondsWithAtom[middle_atom];
04959 
04960         //  Now loop through all the bonds connect to the
04961         //  middle atom
04962         while (*bond2 != -1)
04963         {
04964           if (bonds[*bond2].atom1 == middle_atom)
04965           {
04966             if (i < bonds[*bond2].atom2)
04967             {
04968               exclusionSet.add(Exclusion(i,bonds[*bond2].atom2));
04969             }
04970           }
04971           else
04972           {
04973             if (i < bonds[*bond2].atom1)
04974             {
04975               exclusionSet.add(Exclusion(i,bonds[*bond2].atom1));
04976             }
04977           }
04978 
04979           ++bond2;
04980         }
04981 
04982         ++bond1;
04983       }
04984        }
04985 #endif
04986     }
04987     /*      END OF FUNCTION build13excl      */
04988 
04989     /************************************************************************/
04990     /*                  */
04991     /*        FUNCTION build14excl      */
04992     /*                  */
04993     /************************************************************************/
04994 
04995 
04996     void Molecule::build14excl(int modified)
04997        
04998     {
04999 #ifndef MEM_OPT_VERSION
05000        int32 *bond1, *bond2, *bond3;  //  The two bonds being checked
05001        int mid1, mid2;    //  Middle atoms
05002        register int i;      //  Counter to loop through all atoms
05003        
05004        //  Loop through all the atoms
05005        for (i=0; i<numAtoms; i++)
05006        {  
05007       // Get all the bonds connect directly to atom i
05008       bond1 = bondsWithAtom[i];
05009        
05010       while (*bond1 != -1)
05011       {
05012         if (bonds[*bond1].atom1 == i)
05013         {
05014           mid1=bonds[*bond1].atom2;
05015         }
05016         else
05017         {
05018           mid1=bonds[*bond1].atom1;
05019         }
05020 
05021         bond2 = bondsWithAtom[mid1];
05022 
05023         //  Loop through all the bonds connected to atom mid1
05024         while (*bond2 != -1)
05025         {
05026           if (bonds[*bond2].atom1 == mid1)
05027           {
05028             mid2 = bonds[*bond2].atom2;
05029           }
05030           else
05031           {
05032             mid2 = bonds[*bond2].atom1;
05033           }
05034 
05035           //  Make sure that we don't double back to where
05036           //  we started from.  This causes strange behavior.
05037           //  Trust me, I've been there . . .
05038           if (mid2 == i)
05039           {
05040             ++bond2;
05041             continue;
05042           }
05043 
05044           bond3=bondsWithAtom[mid2];
05045 
05046           //  Loop through all the bonds connected to mid2
05047           while (*bond3 != -1)
05048           {
05049             if (bonds[*bond3].atom1 == mid2)
05050             {
05051               //  Make sure that we don't double back to where
05052               //  we started from.  This causes strange behavior.
05053               //  Trust me, I've been there . . .
05054               //  I added this!!!  Why wasn't it there before?  -JCP
05055               if (bonds[*bond3].atom2 != mid1)
05056               if (i<bonds[*bond3].atom2)
05057               {
05058                  exclusionSet.add(Exclusion(i,bonds[*bond3].atom2,modified));
05059               }
05060             }
05061             else
05062             {
05063               //  Make sure that we don't double back to where
05064               //  we started from.  This causes strange behavior.
05065               //  Trust me, I've been there . . .
05066               //  I added this!!!  Why wasn't it there before?  -JCP
05067               if (bonds[*bond3].atom1 != mid1)
05068               if (i<bonds[*bond3].atom1)
05069               {
05070                  exclusionSet.add(Exclusion(i,bonds[*bond3].atom1,modified));
05071               }
05072             }
05073 
05074             ++bond3;
05075           }
05076 
05077           ++bond2;
05078         }
05079     
05080         ++bond1;
05081       }
05082        }
05083 #endif
05084     }
05085     /*      END OF FUNCTION build14excl      */
05086 
05087 
05088     /************************************************************************/
05089     /*                                                                      */
05090     /*        FUNCTION stripHGroupExcl                                      */
05091     /*                                                                      */
05092     /*  This function removes all exclusions which are entirely             */
05093     /*  within a single hydrogen group.  This assumes that they all         */
05094     /*  exist, which should be true for exclusion policy 1-3 or higher.     */
05095     /*                                                                      */
05096     /************************************************************************/
05097 
05098   void Molecule::stripHGroupExcl(void)
05099   {
05100 #ifdef NAMD_CUDA
05101     return;
05102 #endif
05103 
05104 #ifdef MEM_OPT_VERSION
05105    //don't eliminate the exclusion signatures but only adjusting numCalcExclusions
05106    //which is related to the checksum in reductions. HGroupExcl are not calculated
05107    //at all
05108     int delExclCnt=0;
05109     UniqueSet<Exclusion> strippedExcls;
05110     HydrogenGroup::iterator h_i, h_e, h_j;
05111     h_i = hydrogenGroup.begin();  h_e = hydrogenGroup.end();
05112     for( ; h_i != h_e; ++h_i ) {
05113       for ( h_j = h_i + 1; h_j != h_e && ! h_j->isGP; ++h_j ) {
05114                 if ( h_i->atomID < h_j->atomID )
05115                         strippedExcls.add(Exclusion(h_i->atomID,h_j->atomID));
05116                 else
05117                         strippedExcls.add(Exclusion(h_j->atomID,h_i->atomID));
05118       }
05119     }
05120     
05121     UniqueSetIter<Exclusion> iter(strippedExcls);
05122     for(iter=iter.begin(); iter!=iter.end(); iter++){
05123         int atom1 = iter->atom1;
05124         int atom2 = iter->atom2;
05125         ExclusionSignature *a1Sig = &exclSigPool[eachAtomExclSig[atom1]];
05126         int fullOrMod, idx;
05127         idx = a1Sig->findOffset(atom2-atom1, &fullOrMod);
05128         if(idx==-1) continue;
05129         delExclCnt++;
05130         if(!exclStrippedByFepOrFixedAtoms(atom1, atom2))
05131             numCalcExclusions -= 2;
05132     }
05133 
05134     numTotalExclusions -= delExclCnt;    
05135 
05136     return;
05137 
05138 #if 0
05139         UniqueSet<Exclusion> strippedExcls;
05140         HydrogenGroup::iterator h_i, h_e, h_j;
05141     h_i = hydrogenGroup.begin();  h_e = hydrogenGroup.end();
05142     for( ; h_i != h_e; ++h_i ) {
05143       for ( h_j = h_i + 1; h_j != h_e && ! h_j->isGP; ++h_j ) {
05144                 if ( h_i->atomID < h_j->atomID )
05145                         strippedExcls.add(Exclusion(h_i->atomID,h_j->atomID));
05146                 else
05147                         strippedExcls.add(Exclusion(h_j->atomID,h_i->atomID));
05148       }
05149     }
05150     
05151     vector<int> *strippedList = new vector<int>[numAtoms];
05152     UniqueSetIter<Exclusion> iter(strippedExcls);
05153     for(iter=iter.begin(); iter!=iter.end(); iter++){
05154         int atom1 = iter->atom1;
05155         int atom2 = iter->atom2;
05156         strippedList[atom1].push_back(atom2);
05157         strippedList[atom2].push_back(atom1);
05158     }
05159     
05160     vector<ExclusionSignature> newExclSigPool;
05161     for(int i=0; i<numAtoms; i++){
05162         ExclusionSignature a1Sig = exclSigPool[eachAtomExclSig[i]];
05163         sort(strippedList[i].begin(), strippedList[i].end());           
05164         //eliminate the atoms in the strippedList from the atom signature               
05165         int fullOrMod = 0; //0: full; 1: mod            
05166         vector<int> whichOffs;
05167         vector<int> changedOffs;
05168         for(int j=0; j<strippedList[i].size(); j++){
05169                 int atom2 = strippedList[i][j];
05170                 int offset = atom2-i;
05171             //fullOffset and modOffset should be in increasing order
05172             //for findOffset to be correct. So that assignment has
05173             //to be done later
05174                 int idx = a1Sig.findOffset(offset, &fullOrMod);
05175                 if(idx==-1) continue;
05176                         //NAMD_die("Something wrong in Molecule::stripHGroupExcl!\n");                  
05177                         whichOffs.push_back(fullOrMod);
05178             changedOffs.push_back(idx);
05179         }
05180 
05181         for(int j=0; j<whichOffs.size(); j++){
05182             if(whichOffs[j]==0){ //full exclusion             
05183                                 a1Sig.fullOffset[changedOffs[j]] = 0;
05184                         }else{
05185                                 a1Sig.modOffset[changedOffs[j]] = 0;
05186                         }
05187         }
05188         whichOffs.clear();
05189         changedOffs.clear();
05190 
05191                 if(strippedList[i].size()){
05192                         //indicating a1Sig has changed
05193                         a1Sig.removeEmptyOffset();
05194                         int poolIndex = lookupCstPool(newExclSigPool, a1Sig);
05195                         if(poolIndex==-1){
05196                                 poolIndex = newExclSigPool.size();
05197                                 newExclSigPool.push_back(a1Sig);
05198                         }
05199                         eachAtomExclSig[i] = poolIndex + exclSigPoolSize;                       
05200                 }                                                                   
05201         strippedList[i].clear();
05202     }
05203     
05204     addNewExclSigPool(newExclSigPool);
05205 
05206     newExclSigPool.clear();
05207            
05208     delete [] strippedList;     
05209 #endif   
05210 #else
05211     HydrogenGroup::iterator h_i, h_e, h_j;
05212     h_i = hydrogenGroup.begin();  h_e = hydrogenGroup.end();
05213     for( ; h_i != h_e; ++h_i ) {
05214       for ( h_j = h_i + 1; h_j != h_e && ! h_j->isGP; ++h_j ) {
05215         if ( h_i->atomID < h_j->atomID )
05216           exclusionSet.del(Exclusion(h_i->atomID,h_j->atomID));
05217         else
05218           exclusionSet.del(Exclusion(h_j->atomID,h_i->atomID));
05219       }
05220     }
05221 #endif
05222   }
05223     /*      END OF FUNCTION stripHGroupExcl      */
05224 
05225     /************************************************************************/
05226     /*                                                                      */
05227     /*        FUNCTION stripFepExcl                                         */
05228     /*                                                                      */
05229     /************************************************************************/
05230 
05231 #ifdef MEM_OPT_VERSION
05232   void Molecule::stripFepFixedExcl(void)
05233   {
05234       int delExclCnt=0;
05235 
05236       vector<ExclusionSignature> newExclSigPool;
05237       if(simParams->alchFepOn || simParams->alchThermIntOn || simParams->lesOn){
05238           for(int i=0; i<numAtoms; i++){
05239               int atom1 = i;
05240               int t1 = get_fep_type(atom1);
05241               if(t1==0) continue;
05242               ExclusionSignature a1Sig = exclSigPool[eachAtomExclSig[atom1]];
05243               int sigChanged = 0;
05244               for(int j=0; j<a1Sig.fullExclCnt; j++){
05245                   int atom2 = atom1 + a1Sig.fullOffset[j];
05246                   int t2 = get_fep_type(atom2);
05247                   if(t2 && t1!=t2){
05248                       a1Sig.fullOffset[j] = 0;
05249                       sigChanged = 1;
05250                       delExclCnt++;
05251                   }                  
05252               }
05253               for(int j=0; j<a1Sig.modExclCnt; j++){
05254                   int atom2 = atom1 + a1Sig.modOffset[j];
05255                   int t2 = get_fep_type(atom2);
05256                   if(t2 && t1!=t2){
05257                       a1Sig.modOffset[j] = 0;
05258                       sigChanged = 1;
05259                       delExclCnt++;
05260                   }
05261               }
05262 
05263               //deal with the new signature
05264               if(sigChanged){
05265                   a1Sig.removeEmptyOffset();
05266                   int poolIndex = lookupCstPool(newExclSigPool, a1Sig);
05267                   if(poolIndex==-1){
05268                       poolIndex = newExclSigPool.size();
05269                       newExclSigPool.push_back(a1Sig);
05270                   }
05271                   eachAtomExclSig[atom1] = poolIndex + exclSigPoolSize;
05272               }                            
05273           }
05274           //integerate new exclusion signatures into the old one
05275           addNewExclSigPool(newExclSigPool);
05276           newExclSigPool.clear();
05277       }else if(simParams->pairInteractionOn){
05278           for(int i=0; i<numAtoms; i++){
05279               int atom1 = i;
05280               int t1 = get_fep_type(atom1);              
05281               ExclusionSignature a1Sig = exclSigPool[eachAtomExclSig[atom1]];
05282               int sigChanged = 0;
05283               for(int j=0; j<a1Sig.fullExclCnt; j++){
05284                   int atom2 = atom1 + a1Sig.fullOffset[j];
05285                   int t2 = get_fep_type(atom2);
05286 
05287                   int cond1 = simParams->pairInteractionSelf && (t1!=1 || t2!=1);
05288                   int cond2 = (t1!=1 || t2!=2) && (t1!=2 || t2!=1);
05289                   if(cond1 || cond2){                                   
05290                       a1Sig.fullOffset[j] = 0;
05291                       sigChanged = 1;
05292                       delExclCnt++;
05293                   }                  
05294               }
05295               for(int j=0; j<a1Sig.modExclCnt; j++){
05296                   int atom2 = atom1 + a1Sig.modOffset[j];
05297                   int t2 = get_fep_type(atom2);
05298                   int cond1 = simParams->pairInteractionSelf && (t1!=1 || t2!=1);
05299                   int cond2 = (t1!=1 || t2!=2) && (t1!=2 || t2!=1);
05300                   if(cond1 || cond2){                  
05301                       a1Sig.modOffset[j] = 0;
05302                       sigChanged = 1;
05303                       delExclCnt++;
05304                   }
05305               }
05306 
05307               //deal with the new signature
05308               if(sigChanged){
05309                   a1Sig.removeEmptyOffset();
05310                   int poolIndex = lookupCstPool(newExclSigPool, a1Sig);
05311                   if(poolIndex==-1){
05312                       poolIndex = newExclSigPool.size();
05313                       newExclSigPool.push_back(a1Sig);
05314                   }
05315                   eachAtomExclSig[atom1] = poolIndex + exclSigPoolSize;
05316               }                            
05317           }
05318           //integerate new exclusion signatures into the old one
05319           addNewExclSigPool(newExclSigPool);
05320           newExclSigPool.clear();
05321       }
05322 
05323       numTotalExclusions -= (delExclCnt/2); 
05324 
05325       //deal with fixed atoms
05326     //numCalcExclusions has to be calculated in stripHGroupExcl
05327       //numCalcExclusions = 0;
05328       if(!numFixedAtoms){
05329         for(int i=0; i<numAtoms; i++){
05330             ExclusionSignature *a1Sig = &exclSigPool[eachAtomExclSig[i]];
05331             //numCalcExclusions += (a1Sig->fullExclCnt + a1Sig->modExclCnt);
05332         }
05333       }else{
05334           for(int i=0; i<numAtoms; i++){
05335               int atom1 = i;
05336               int t1 = fixedAtomFlags[atom1];
05337               ExclusionSignature a1Sig = exclSigPool[eachAtomExclSig[atom1]];
05338               if(t1==0){
05339                   //numCalcExclusions += (a1Sig.fullExclCnt + a1Sig.modExclCnt);
05340                   continue;
05341               }              
05342               int sigChanged = 0;
05343               for(int j=0; j<a1Sig.fullExclCnt; j++){
05344                   int atom2 = atom1 + a1Sig.fullOffset[j];
05345                   int t2 = fixedAtomFlags[atom2];
05346                   if(t2){
05347                       a1Sig.fullOffset[j] = 0;
05348                       sigChanged = 1;
05349                       delExclCnt++;
05350                   }                  
05351               }
05352               for(int j=0; j<a1Sig.modExclCnt; j++){
05353                   int atom2 = atom1 + a1Sig.modOffset[j];
05354                   int t2 = fixedAtomFlags[atom2];
05355                   if(t2){
05356                       a1Sig.modOffset[j] = 0;
05357                       sigChanged = 1;
05358                       delExclCnt++;
05359                   }
05360               }
05361 
05362               //deal with the new signature
05363               if(sigChanged){
05364                   a1Sig.removeEmptyOffset();
05365                   int poolIndex = lookupCstPool(newExclSigPool, a1Sig);
05366                   if(poolIndex==-1){
05367                       poolIndex = newExclSigPool.size();
05368                       newExclSigPool.push_back(a1Sig);
05369                   }
05370                   eachAtomExclSig[atom1] = poolIndex + exclSigPoolSize;
05371               }              
05372               //numCalcExclusions += (a1Sig.fullExclCnt + a1Sig.modExclCnt);                            
05373           }
05374           //integerate new exclusion signatures into the old one
05375           addNewExclSigPool(newExclSigPool);
05376           newExclSigPool.clear();
05377       }      
05378 
05379       numCalcExclusions = (numCalcExclusions - delExclCnt)/2;      
05380   }
05381 #else
05382   void Molecule::stripFepExcl(void)
05383   {      
05384     UniqueSet<Exclusion> fepExclusionSet;
05385     UniqueSetIter<Exclusion> exclIter(exclusionSet);
05386 
05387     if ( simParams->alchFepOn || simParams->alchThermIntOn || simParams->lesOn ) {
05388        for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
05389        {
05390          int t1 = get_fep_type(exclIter->atom1);
05391          int t2 = get_fep_type(exclIter->atom2);
05392          if ( t1 && t2 && t1 != t2 ) {
05393            fepExclusionSet.add(*exclIter);
05394          }
05395        }
05396     } else if ( simParams->pairInteractionOn ) {
05397       for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
05398       {
05399         int ifep_type = get_fep_type(exclIter->atom1);
05400         int jfep_type = get_fep_type(exclIter->atom2);
05401         if ( simParams->pairInteractionSelf ) {
05402           // for pair-self, both atoms must be in group 1.
05403           if (ifep_type != 1 || jfep_type != 1) {
05404             fepExclusionSet.add(*exclIter);
05405           }
05406         } else {
05407 
05408           // for pair, must have one from each group.
05409           if (!(ifep_type == 1 && jfep_type == 2) &&
05410               !(ifep_type == 2 && jfep_type == 1)) {
05411             fepExclusionSet.add(*exclIter);
05412           }
05413         }
05414        }
05415     }
05416 
05417     UniqueSetIter<Exclusion> fepIter(fepExclusionSet);
05418     for ( fepIter=fepIter.begin(); fepIter != fepIter.end(); fepIter++ )
05419     {
05420       exclusionSet.del(*fepIter);
05421     }
05422   }
05423 #endif
05424     /*      END OF FUNCTION stripFepExcl      */
05425 
05426 
05427 #endif  // MOLECULE2_C undefined = first object file
05428 #ifdef MOLECULE2_C  // second object file
05429 
05430 /* BEGIN gf */
05431     /************************************************************************/
05432     /*                                                                      */
05433     /*      FUNCTION build_gridforce_params                                 */
05434     /*                                                                      */
05435     /*   INPUTS:                                                            */
05436     /*  gridfrcfile - Value of gridforcefile from config file               */
05437     /*  gridfrccol - Value of gridforcecol from config file                 */
05438     /*  gridfrcchrgcol - Value of gridforcechargecol from config file       */
05439     /*  potfile - Value of gridforcepotfile from config file                  */
05440     /*  initial_pdb - PDB object that contains initial positions            */
05441     /*  cwd - Current working directory                                     */
05442     /*                                                                      */
05443     // This function builds all the parameters that are necessary to
05444     // do gridforcing. This involves looking through a PDB object to
05445     // determine which atoms are to be gridforced, and what the force
05446     // multiplier is for each atom.  This information is then stored
05447     // in the arrays gridfrcIndexes and gridfrcParams.
05448     /************************************************************************/
05449 
05450 void Molecule::build_gridforce_params(StringList *gridfrcfile,
05451                                       StringList *gridfrccol,
05452                                       StringList *gridfrcchrgcol,
05453                                       StringList *potfile,
05454                                       PDB *initial_pdb,
05455                                       char *cwd)
05456 {
05457     PDB *kPDB;
05458     register int i;             //  Loop counters
05459     register int j;
05460     register int k;
05461 
05462     DebugM(3,  "Entered build_gridforce_params multi...\n");
05463 //     DebugM(3, "\tgridfrcfile = " << gridfrcfile->data << endi);
05464 //     DebugM(3, "\tgridfrccol = " << gridfrccol->data << endi);
05465     
05466     MGridforceParams* mgridParams = simParams->mgridforcelist.get_first();
05467     numGridforceGrids = 0;
05468     while (mgridParams != NULL) {
05469         numGridforceGrids++;
05470         mgridParams = mgridParams->next;
05471     }
05472     
05473     gridfrcIndexes = new int32*[numGridforceGrids];
05474     gridfrcParams = new GridforceParams*[numGridforceGrids];
05475     gridfrcGrid = new GridforceGrid*[numGridforceGrids];
05476     numGridforces = new int[numGridforceGrids];
05477     
05478     mgridParams = simParams->mgridforcelist.get_first();
05479     for (int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
05480         int current_index=0;    //  Index into values used
05481         int kcol = 5;           //  Column to look for force constant in
05482         int qcol = 0;           //  Column for charge (default 0: use electric charge)
05483         Real kval = 0;          //  Force constant value retreived
05484         char filename[129];     //  PDB filename
05485         char potfilename[129];  //  Potential file name
05486         
05487         if (mgridParams == NULL) {
05488             NAMD_die("Problem with mgridParams!");
05489         }
05490         
05491         // Now load values from mgridforcelist object
05492         if (mgridParams->gridforceFile == NULL)
05493         {
05494             kPDB = initial_pdb;
05495         }
05496         else
05497         {
05498             DebugM(4, "mgridParams->gridforceFile = " << mgridParams->gridforceFile << "\n" << endi);
05499             
05500             if ( (cwd == NULL) || (mgridParams->gridforceFile[0] == '/') )
05501             {
05502                 strcpy(filename, mgridParams->gridforceFile);
05503             }
05504             else
05505             {
05506                 strcpy(filename, cwd);
05507                 strcat(filename, mgridParams->gridforceFile);
05508             }
05509         
05510             kPDB = new PDB(filename);
05511             if ( kPDB == NULL )
05512             {
05513                 NAMD_die("Memory allocation failed in Molecule::build_gridforce_params");
05514             }
05515            
05516             if (kPDB->num_atoms() != numAtoms)
05517             {
05518                 NAMD_die("Number of atoms in grid force PDB doesn't match coordinate PDB");
05519             }
05520         }
05521 
05522         //  Get the column that the force constant is going to be in.  It
05523         //  can be in any of the 5 floating point fields in the PDB, according
05524         //  to what the user wants.  The allowable fields are X, Y, Z, O, or
05525         //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
05526         //  The default is the 5th field, which is beta (temperature factor)
05527         if (mgridParams->gridforceCol == NULL)
05528         {
05529             kcol = 5;
05530         }
05531         else
05532         {
05533             if (strcasecmp(mgridParams->gridforceCol, "X") == 0)
05534             {
05535                 kcol=1;
05536             }
05537             else if (strcasecmp(mgridParams->gridforceCol, "Y") == 0)
05538             {
05539                 kcol=2;
05540             }
05541             else if (strcasecmp(mgridParams->gridforceCol, "Z") == 0)
05542             {
05543                 kcol=3;
05544             }
05545             else if (strcasecmp(mgridParams->gridforceCol, "O") == 0)
05546             {
05547                 kcol=4;
05548             }
05549             else if (strcasecmp(mgridParams->gridforceCol, "B") == 0)
05550             {
05551                 kcol=5;
05552             }
05553             else
05554             {
05555                 NAMD_die("gridforcecol must have value of X, Y, Z, O, or B");
05556             }
05557         }
05558     
05559         //  Get the column that the charge is going to be in.
05560         if (mgridParams->gridforceQcol == NULL)
05561         {
05562             qcol = 0;   // Default: don't read charge from file, use electric charge
05563         }
05564         else
05565         {
05566             if (strcasecmp(mgridParams->gridforceQcol, "X") == 0)
05567             {
05568                 qcol=1;
05569             }
05570             else if (strcasecmp(mgridParams->gridforceQcol, "Y") == 0)
05571             {
05572                 qcol=2;
05573             }
05574             else if (strcasecmp(mgridParams->gridforceQcol, "Z") == 0)
05575             {
05576                 qcol=3;
05577             }
05578             else if (strcasecmp(mgridParams->gridforceQcol, "O") == 0)
05579             {
05580                 qcol=4;
05581             }
05582             else if (strcasecmp(mgridParams->gridforceQcol, "B") == 0)
05583             {
05584                 qcol=5;
05585             }
05586             else
05587             {
05588                 NAMD_die("gridforcechargecol must have value of X, Y, Z, O, or B");
05589             }
05590         }
05591     
05592         if (kcol == qcol) {
05593             NAMD_die("gridforcecol and gridforcechargecol cannot have same value");
05594         }
05595 
05596     
05597         //  Allocate an array that will store an index into the constraint
05598         //  parameters for each atom.  If the atom is not constrained, its
05599         //  value will be set to -1 in this array.
05600         DebugM(4, "Here!\n" << endi);
05601         gridfrcIndexes[gridnum] = new int32[numAtoms];
05602         DebugM(4, "There!\n" << endi);
05603        
05604         if (gridfrcIndexes[gridnum] == NULL)
05605         {
05606             NAMD_die("memory allocation failed in Molecule::build_gridforce_params()");
05607         }
05608 
05609         //  Loop through all the atoms and find out which ones are constrained
05610         for (i=0; i<numAtoms; i++)
05611         {
05612             //  Get the k value based on where we were told to find it
05613             switch (kcol)
05614             {
05615             case 1:
05616                 kval = (kPDB->atom(i))->xcoor();
05617                 break;
05618             case 2:
05619                 kval = (kPDB->atom(i))->ycoor();
05620                 break;
05621             case 3:
05622                 kval = (kPDB->atom(i))->zcoor();
05623                 break;
05624             case 4:
05625                 kval = (kPDB->atom(i))->occupancy();
05626                 break;
05627             case 5:
05628                 kval = (kPDB->atom(i))->temperaturefactor();
05629                 break;
05630             }
05631            
05632             if (kval > 0.0)
05633             {
05634                 //  This atom is constrained
05635                 gridfrcIndexes[gridnum][i] = current_index;
05636                 current_index++;
05637             }
05638             else
05639             {
05640                 //  This atom is not constrained
05641                 gridfrcIndexes[gridnum][i] = -1;
05642             }
05643         }
05644     
05645         if (current_index == 0)
05646         {
05647             //  Constraints were turned on, but there weren't really any constrained
05648             iout << iWARN << "NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" << endi;
05649         }
05650         else
05651         {
05652             //  Allocate an array to hold the constraint parameters
05653             gridfrcParams[gridnum] = new GridforceParams[current_index];
05654             if (gridfrcParams[gridnum] == NULL)
05655             {
05656                 NAMD_die("memory allocation failed in Molecule::build_gridforce_params");
05657             }
05658         }
05659     
05660         numGridforces[gridnum] = current_index;
05661 
05662         //  Loop through all the atoms and assign the parameters for those
05663         //  that are constrained
05664         for (i=0; i<numAtoms; i++)
05665         {
05666             if (gridfrcIndexes[gridnum][i] != -1)
05667             {
05668                 //  This atom has grid force, so get the k value again
05669                 switch (kcol)
05670                 {
05671                 case 1:
05672                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->xcoor();
05673                     break;
05674                 case 2:
05675                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->ycoor();
05676                     break;
05677                 case 3:
05678                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->zcoor();
05679                     break;
05680                 case 4:
05681                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->occupancy();
05682                     break;
05683                 case 5:
05684                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->atom(i))->temperaturefactor();
05685                     break;
05686                 }
05687             
05688                 //  Also get charge column
05689                 switch (qcol)
05690                 {
05691                 case 0:
05692 #ifdef MEM_OPT_VERSION
05693                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
05694 #else
05695                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
05696 #endif
05697                     break;
05698                 case 1:
05699                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->xcoor();
05700                     break;
05701                 case 2:
05702                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->ycoor();
05703                     break;
05704                 case 3:
05705                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->zcoor();
05706                     break;
05707                 case 4:
05708                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->occupancy();
05709                     break;
05710                 case 5:
05711                     gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->atom(i))->temperaturefactor();
05712                     break;
05713                 }
05714             }
05715         }
05716        
05717         //  If we had to create new PDB objects, delete them now
05718         if (mgridParams->gridforceFile != NULL)
05719         {
05720             delete kPDB;
05721         }
05722     
05723         //  Now we fill in our grid information
05724     
05725         // Open potential file
05726         if ( (cwd == NULL) || (mgridParams->gridforceVfile[0] == '/') )
05727         {
05728             strcpy(potfilename, mgridParams->gridforceVfile);
05729         }
05730         else
05731         {
05732             strcpy(potfilename, cwd);
05733             strcat(potfilename, mgridParams->gridforceVfile);
05734         }
05735     
05736 //        iout << iINFO << "Allocating grid " << gridnum
05737 //             << "\n" << endi;
05738              
05739         gridfrcGrid[gridnum] = new GridforceGrid(gridnum);
05740         //      gridfrcGrid[gridnum]->initialize(potfilename, simParams, mgridParams);
05741         gridfrcGrid[gridnum]->init1(potfilename, simParams, mgridParams);
05742 //        ComputeGridForceNodeMgr* mgr = CProxy_ComputeGridForceNodeMgr::
05743 //          ckLocalBranch(CkpvAccess(BOCclass_group).
05744 //          computeGridForceNodeMgr);
05745 //      mgr->depositInitialGrid(gridnum,gridfrcGrid[gridnum],numGridforceGrids);
05746 
05747         CProxy_ComputeGridForceNodeMgr 
05748           mgr(CkpvAccess(BOCclass_group).computeGridForceNodeMgr);
05749           
05750         GridDepositMsg *outmsg = new GridDepositMsg;
05751         outmsg->gridnum = gridnum;
05752         outmsg->grid = gridfrcGrid[gridnum];
05753         outmsg->num_grids = numGridforceGrids;
05754         mgr[CkMyNode()].depositInitialGrid(outmsg);        
05755  
05756         // Finally, get next mgridParams pointer
05757         mgridParams = mgridParams->next;
05758     }
05759 }
05760 /* END gf */
05761 
05762 
05763 
05764     /************************************************************************/
05765     /*                  */
05766     /*      FUNCTION build_constraint_params    */
05767     /*                  */
05768     /*   INPUTS:                */
05769     /*  consref - Value of consref parameter from config file    */
05770     /*  conskfile - Value of conskfile from config file      */
05771     /*  conskcol - Value of conskcol from config file      */
05772     /*  initial_pdb - PDB object that contains initial positions  */
05773     /*  cwd - Current working directory          */
05774     /*                  */
05775     /*  This function builds all the parameters that are necessary  */
05776     /*   to do harmonic constraints.  This involves looking through    */
05777     /*   one or more PDB objects to determine which atoms are constrained,  */
05778     /*   and what the force constant and reference position is force each   */
05779     /*   atom that is constrained.  This information is then stored    */
05780     /*   in the arrays consIndexes and consParams.        */
05781     /*                  */
05782     /************************************************************************/
05783 
05784     void Molecule::build_constraint_params(StringList *consref, 
05785              StringList *conskfile, 
05786              StringList *conskcol, 
05787              PDB *initial_pdb,
05788              char *cwd)
05789        
05790     {
05791        PDB *refPDB, *kPDB;    //  Pointer to other PDB's if used
05792        register int i;      //  Loop counter
05793        int current_index=0;    //  Index into values used
05794        int kcol = 4;      //  Column to look for force constant in
05795        Real kval = 0;      //  Force constant value retreived
05796        char filename[129];    //  PDB filename
05797        
05798        //  Get the PDB object that contains the reference positions.  If
05799        //  the user gave another file name, use it.  Otherwise, just use
05800        //  the PDB file that has the initial coordinates.  i.e., constrain
05801        //  the atoms around their initial position.  This is the most likely
05802        //  case anyway
05803        if (consref == NULL)
05804        {
05805     refPDB = initial_pdb;
05806        }
05807        else
05808        {
05809     if (consref->next != NULL)
05810     {
05811        NAMD_die("Multiple definitions of constraint reference file in configruation file");
05812     }
05813 
05814     if ( (cwd == NULL) || (consref->data[0] == '/') )
05815     {
05816          strcpy(filename, consref->data);
05817     }
05818     else
05819     {
05820          strcpy(filename, cwd);
05821          strcat(filename, consref->data);
05822     }
05823     
05824     refPDB = new PDB(filename);
05825     if ( refPDB == NULL )
05826     {
05827       NAMD_die("Memory allocation failed in Molecule::build_constraint_params");
05828     }
05829     
05830     if (refPDB->num_atoms() != numAtoms)
05831     {
05832        NAMD_die("Number of atoms in constraint reference PDB doesn't match coordinate PDB");
05833     }
05834        }
05835        
05836        //  Get the PDB to read the force constants from.  Again, if the user
05837        //  gave us another file name, open that one.  Otherwise, just use
05838        //  the PDB with the initial coordinates
05839        if (conskfile == NULL)
05840        {
05841     kPDB = initial_pdb;
05842        }
05843        else
05844        {
05845     if (conskfile->next != NULL)
05846     {
05847        NAMD_die("Multiple definitions of constraint constant file in configuration file");
05848     }
05849 
05850     if ( (consref != NULL) && (strcasecmp(consref->data, conskfile->data) == 0) )
05851     {
05852        //  Same PDB used for reference positions and force constants
05853        kPDB = refPDB; 
05854     }
05855     else
05856     {
05857       if ( (cwd == NULL) || (conskfile->data[0] == '/') )
05858       {
05859         strcpy(filename, conskfile->data);
05860       }
05861       else
05862       {
05863         strcpy(filename, cwd);
05864         strcat(filename, conskfile->data);
05865       }
05866 
05867       kPDB = new PDB(filename);
05868       if ( kPDB == NULL )
05869       {
05870         NAMD_die("Memory allocation failed in Molecule::build_constraint_params");
05871       }
05872     
05873       if (kPDB->num_atoms() != numAtoms)
05874       {
05875          NAMD_die("Number of atoms in constraint constant PDB doesn't match coordinate PDB");
05876       }
05877     }
05878        }
05879        
05880        //  Get the column that the force constant is going to be in.  It
05881        //  can be in any of the 5 floating point fields in the PDB, according
05882        //  to what the user wants.  The allowable fields are X, Y, Z, O, or
05883        //  B which correspond to the 1st, 2nd, ... 5th floating point fields.
05884        //  The default is the 4th field, which is the occupancy
05885        if (conskcol == NULL)
05886        {
05887     kcol = 4;
05888        }
05889        else
05890        {
05891     if (conskcol->next != NULL)
05892     {
05893        NAMD_die("Multiple definitions of harmonic constraint column in config file");
05894     }
05895     
05896     if (strcasecmp(conskcol->data, "X") == 0)
05897     {
05898        kcol=1;
05899     }
05900     else if (strcasecmp(conskcol->data, "Y") == 0)
05901     {
05902        kcol=2;
05903     }
05904     else if (strcasecmp(conskcol->data, "Z") == 0)
05905     {
05906        kcol=3;
05907     }
05908     else if (strcasecmp(conskcol->data, "O") == 0)
05909     {
05910        kcol=4;
05911     }
05912     else if (strcasecmp(conskcol->data, "B") == 0)
05913     {
05914        kcol=5;
05915     }
05916     else
05917     {
05918        NAMD_die("conskcol must have value of X, Y, Z, O, or B");
05919     }
05920        }
05921        
05922        //  Allocate an array that will store an index into the constraint
05923        //  parameters for each atom.  If the atom is not constrained, its
05924        //  value will be set to -1 in this array.
05925        consIndexes = new int32[numAtoms];
05926        
05927        if (consIndexes == NULL)
05928        {
05929     NAMD_die("memory allocation failed in Molecule::build_constraint_params()");
05930        }
05931        
05932        //  Loop through all the atoms and find out which ones are constrained
05933        for (i=0; i<numAtoms; i++)
05934        {
05935     //  Get the k value based on where we were told to find it
05936     switch (kcol)
05937     {
05938        case 1:
05939     kval = (kPDB->atom(i))->xcoor();
05940     break;
05941        case 2:
05942     kval = (kPDB->atom(i))->ycoor();
05943     break;
05944        case 3:
05945     kval = (kPDB->atom(i))->zcoor();
05946     break;
05947        case 4:
05948     kval = (kPDB->atom(i))->occupancy();
05949     break;
05950        case 5:
05951     kval = (kPDB->atom(i))->temperaturefactor();
05952     break;
05953     }
05954     
05955     if (kval > 0.0)
05956     {
05957        //  This atom is constrained
05958        consIndexes[i] = current_index;
05959        current_index++;
05960     }
05961     else
05962     {
05963        //  This atom is not constrained
05964        consIndexes[i] = -1;
05965     }
05966        }
05967        
05968        if (current_index == 0)
05969        {
05970     //  Constraints were turned on, but there weren't really any constrained
05971     iout << iWARN << "NO CONSTRAINED ATOMS WERE FOUND, BUT CONSTRAINTS ARE ON . . .\n" << endi;
05972        }
05973        else
05974        {
05975     //  Allocate an array to hold the constraint parameters
05976     consParams = new ConstraintParams[current_index];
05977     
05978     if (consParams == NULL)
05979     {
05980        NAMD_die("memory allocation failed in Molecule::build_constraint_params");
05981     }
05982        }
05983        
05984        numConstraints = current_index;
05985        
05986        //  Loop through all the atoms and assign the parameters for those
05987        //  that are constrained
05988        for (i=0; i<numAtoms; i++)
05989        {
05990     if (consIndexes[i] != -1)
05991     {
05992        //  This atom is constrained, so get the k value again
05993        switch (kcol)
05994        {
05995           case 1:
05996        consParams[consIndexes[i]].k = (kPDB->atom(i))->xcoor();
05997        break;
05998           case 2:
05999        consParams[consIndexes[i]].k = (kPDB->atom(i))->ycoor();
06000        break;
06001           case 3:
06002        consParams[consIndexes[i]].k = (kPDB->atom(i))->zcoor();
06003        break;
06004           case 4:
06005        consParams[consIndexes[i]].k = (kPDB->atom(i))->occupancy();
06006        break;
06007           case 5:
06008        consParams[consIndexes[i]].k = (kPDB->atom(i))->temperaturefactor();
06009        break;
06010        }
06011        
06012        //  Get the reference position
06013        consParams[consIndexes[i]].refPos.x = (refPDB->atom(i))->xcoor();
06014        consParams[consIndexes[i]].refPos.y = (refPDB->atom(i))->ycoor();
06015        consParams[consIndexes[i]].refPos.z = (refPDB->atom(i))->zcoor();
06016     }
06017        }
06018        
06019        //  If we had to create new PDB objects, delete them now
06020        if (consref != NULL)
06021        {
06022     delete refPDB;
06023        }
06024        
06025        if ((conskfile != NULL) &&
06026      !((consref != NULL) && 
06027        (strcasecmp(consref->data, conskfile->data) == 0)
06028       )
06029     )
06030        {
06031     delete kPDB;
06032        }
06033 
06034     }
06035     /*      END OF FUNCTION build_constraint_params    */
06036 
06037 
06038 /************************************************************************/
06039 /*                  */
06040 /*      FUNCTION build_movdrag_params  */
06041 /*                  */
06042 /*   INPUTS:        */
06043 /*  movDragFile - value of movDragFile from the config file */
06044 /*  movDragCol - value of movDragCol from the config file */
06045 /*  movDragVelFile - value of movDragVelFile from the config file */
06046 /*  initial_pdb - PDB object that contains initial positions  */
06047 /*  cwd - Current working directory          */
06048 /*                  */
06049 /*  This function builds all the parameters that are necessary  */
06050 /*  to do moving drag. This involves looking through one or more    */
06051 /*  PDB objects to determine which atoms are dragged,  and what the */
06052 /*  drag parameters for each atom are. This information is then stored */
06053 /*  in the arrays movDragIndexes and movDragParams. */
06054 /*                  */
06055 /************************************************************************/
06056 
06057 void Molecule::build_movdrag_params(StringList *movDragFile, 
06058                                     StringList *movDragCol, 
06059                                     StringList *movDragVelFile, 
06060                                     PDB *initial_pdb,
06061                                     char *cwd)
06062   
06063 {
06064   PDB *tPDB, *vPDB;        //  Pointers to other PDB file(s)
06065   register int i;          //  Loop counter
06066   int current_index=0;     //  Index into values used
06067   int dtcol = 4;           //  Column to look for drag tag in
06068   Real dtval = 0;          //  Drag tag value retreived
06069   char mainfilename[129];  //  main moving drag PDB filename
06070   char velfilename[129];   //  moving drag velocity PDB filename
06071   
06072   //  Get the PDB to read the moving drag tags from. Again, if the
06073   //  user gave us another file name, open that one.  Otherwise, just
06074   //  use the PDB with the initial coordinates
06075   if (movDragFile == NULL) {
06076     tPDB = initial_pdb;
06077     
06078   } else {
06079 
06080     if (movDragFile->next != NULL) {
06081       NAMD_die("Multiple definitions of moving drag tag file in configuration file");
06082     }
06083     
06084     if ( (cwd == NULL) || (movDragFile->data[0] == '/') ) {
06085       strcpy(mainfilename, movDragFile->data);
06086     } else {
06087       strcpy(mainfilename, cwd);
06088       strcat(mainfilename, movDragFile->data);
06089       }
06090     
06091     tPDB = new PDB(mainfilename);
06092     if ( tPDB == NULL ) {
06093       NAMD_die("Memory allocation failed in Molecule::build_movdrag_params");
06094     }
06095     
06096     if (tPDB->num_atoms() != numAtoms) {
06097       NAMD_die("Number of atoms in moving drag tag PDB doesn't match coordinate PDB");
06098     }
06099   }
06100   
06101   // Get the PDB to read atom velocities. If no name given, use
06102   // movDragFile if it is defined. Can NOT use the PDB coordinate
06103   // file!
06104   
06105   if (movDragVelFile == NULL) {
06106     if (movDragFile == NULL) {
06107       NAMD_die("Moving drag velocity file can not be same as coordinate PDB file");
06108     } else {
06109       if (movDragVelFile->next != NULL) {
06110         NAMD_die("Multiple definitions of moving drag velocity file in configuration file");
06111       };
06112       vPDB = tPDB;
06113     };
06114 
06115   } else {
06116 
06117     if ( (cwd == NULL) || (movDragVelFile->data[0] == '/') ) {
06118       strcpy(velfilename, movDragVelFile->data);
06119     } else {
06120       strcpy(velfilename, cwd);
06121       strcat(velfilename, movDragVelFile->data);
06122     }
06123     
06124     vPDB = new PDB(velfilename);
06125     if ( vPDB == NULL ) {
06126       NAMD_die("Memory allocation failed in Molecule::build_movdrag_params");
06127     }
06128     
06129     if (vPDB->num_atoms() != numAtoms) {
06130       NAMD_die("Number of atoms in moving drag velocity PDB doesn't match coordinate PDB");
06131     }
06132   };
06133   
06134   
06135   //  Get the column that the drag tag is going to be in. If
06136   //  movDragFile is defined, it can be in any of the 5 floating point
06137   //  fields in the PDB (X, Y, Z, O, or B) which correspond to the
06138   //  1st, 2nd, ... 5th floating point fields. If movDragFile is NOT
06139   //  defined, it can only be O or B fileds. The default is the O
06140   //  (4th) field, which is the occupancy.
06141 
06142   if (movDragCol == NULL) {
06143     dtcol = 4;
06144   } else {
06145     if (movDragCol->next != NULL) {
06146       NAMD_die("Multiple definitions of drag column in config file");
06147     };
06148     
06149     if (movDragFile == NULL
06150         && strcasecmp(movDragCol->data, "B")
06151         && strcasecmp(movDragCol->data, "O")) {
06152       NAMD_die("Can not read moving drag tags from X, Y, or Z column of the coordinate or velocity file");
06153     };
06154     if (!strcasecmp(movDragCol->data, "X")) {
06155       dtcol=1;
06156     } else if (!strcasecmp(movDragCol->data, "Y")) {
06157       dtcol=2;
06158     } else if (!strcasecmp(movDragCol->data, "Z")) {
06159       dtcol=3;
06160     } else if (!strcasecmp(movDragCol->data, "O")) {
06161       dtcol=4;
06162     } else if (!strcasecmp(movDragCol->data, "B")) {
06163       dtcol=5;
06164     }
06165     else {
06166       NAMD_die("movDragCol must have value of X, Y, Z, O, or B");
06167     };
06168   };
06169   
06170   //  Allocate an array that will store an index into the drag
06171   //  parameters for each atom.  If the atom is not dragged, its
06172   //  value will be set to -1 in this array.
06173   movDragIndexes = new int32[numAtoms];
06174     if (movDragIndexes == NULL) {
06175     NAMD_die("memory allocation failed in Molecule::build_movdrag_params()");
06176   };
06177   
06178   //  Loop through all the atoms and find out which ones are dragged
06179   for (i=0; i<numAtoms; i++) {
06180     switch (dtcol) {
06181     case 1:
06182       dtval = (tPDB->atom(i))->xcoor();
06183       break;
06184     case 2:
06185       dtval = (tPDB->atom(i))->ycoor();
06186       break;
06187     case 3:
06188       dtval = (tPDB->atom(i))->zcoor();
06189       break;
06190     case 4:
06191       dtval = (tPDB->atom(i))->occupancy();
06192       break;
06193     case 5:
06194       dtval = (tPDB->atom(i))->temperaturefactor();
06195       break;
06196     }
06197     
06198     if (dtval != 0.0) {
06199       //  This atom is dragged
06200       movDragIndexes[i] = current_index;
06201       current_index++;
06202     } else {
06203       //  This atom is not dragged
06204       movDragIndexes[i] = -1;
06205     }
06206   }
06207   
06208   if (current_index == 0) {
06209     //  Drag was turned on, but there weren't really any dragged
06210     iout << iWARN << "NO DRAGGED ATOMS WERE FOUND, BUT MOVING DRAG IS ON . . . " << endi;
06211   } else {
06212     //  Allocate an array to hold the drag parameters
06213     movDragParams = new MovDragParams[current_index];
06214     if (movDragParams == NULL) {
06215       NAMD_die("memory allocation failed in Molecule::build_movdrag_params");
06216     }
06217   };
06218   
06219   numMovDrag = current_index;
06220   
06221   //  Loop through all the atoms and assign the parameters for those
06222   //  that are dragged
06223   for (i=0; i<numAtoms; i++) {
06224     if (movDragIndexes[i] != -1) {
06225       movDragParams[movDragIndexes[i]].v[0] = (vPDB->atom(i))->xcoor();
06226       movDragParams[movDragIndexes[i]].v[1] = (vPDB->atom(i))->ycoor();
06227       movDragParams[movDragIndexes[i]].v[2] = (vPDB->atom(i))->zcoor();
06228     };
06229   };
06230       
06231   if (movDragFile != NULL) delete tPDB;
06232   if (movDragVelFile != NULL) delete vPDB;
06233 }
06234 /*      END OF FUNCTION build_movdrag_params    */
06235 
06236 
06237 /************************************************************************/
06238 /*                  */
06239 /*      FUNCTION build_rotdrag_params  */
06240 /*                  */
06241 /*   INPUTS:        */
06242 /*  rotDragFile - value of rotDragFile from the config file */
06243 /*  rotDragCol - value of rotDragCol from the config file */
06244 /*  rotDragAxisFile - value of rotDragAxisFile from the config file */
06245 /*  rotDragPivotFile - value of rotDragPivotFile from the config file */
06246 /*  rotDragVelFile - value of rotDragVelFile from the config file */
06247 /*  rotDragVelCol - value of rotDragVelCol from the config file */
06248 /*  initial_pdb - PDB object that contains initial positions  */
06249 /*  cwd - Current working directory          */
06250 /*                  */
06251 /*  This function builds all the parameters that are necessary  */
06252 /*  to do moving drag. This involves looking through one or more    */
06253 /*  PDB objects to determine which atoms are dragged,  and what the */
06254 /*  drag parameters for each atom are. This information is then stored */
06255 /*  in the arrays rotDragIndexes and rotDragParams. */
06256 /*                  */
06257 /************************************************************************/
06258 
06259 void Molecule::build_rotdrag_params(StringList *rotDragFile, 
06260                                     StringList *rotDragCol, 
06261                                     StringList *rotDragAxisFile, 
06262                                     StringList *rotDragPivotFile, 
06263                                     StringList *rotDragVelFile, 
06264                                     StringList *rotDragVelCol, 
06265                                     PDB *initial_pdb,
06266                                     char *cwd)
06267   
06268 {
06269   PDB *tPDB, *aPDB, *pPDB, *vPDB; //  Pointers to other PDB file(s)
06270   register int i;          //  Loop counter
06271   int current_index=0;     //  Index into values used
06272   int dtcol = 4;           //  Column to look for drag tag in
06273   Real dtval = 0;          //  Drag tag value retreived
06274   int dvcol = 4;           //  Column to look for angular velocity in
06275   Real dvval = 0;          //  Angular velocity value retreived
06276   char mainfilename[129];  //  main rotating drag PDB filename
06277   char axisfilename[129];  //  rotating drag axis PDB filename
06278   char pivotfilename[129]; //  rotating drag pivot point PDB filename
06279   char velfilename[129];   //  rotating drag angular velocity PDB filename
06280   
06281   //  Get the PDB to read the rotating drag tags from. Again, if the
06282   //  user gave us another file name, open that one.  Otherwise, just
06283   //  use the PDB with the initial coordinates
06284   if (rotDragFile == NULL) {
06285     tPDB = initial_pdb;
06286     
06287   } else {
06288 
06289     if (rotDragFile->next != NULL) {
06290       NAMD_die("Multiple definitions of rotating drag tag file in configuration file");
06291     }
06292     
06293     if ( (cwd == NULL) || (rotDragFile->data[0] == '/') ) {
06294       strcpy(mainfilename, rotDragFile->data);
06295     } else {
06296       strcpy(mainfilename, cwd);
06297       strcat(mainfilename, rotDragFile->data);
06298       }
06299     
06300     tPDB = new PDB(mainfilename);
06301     if ( tPDB == NULL ) {
06302       NAMD_die("Memory allocation failed in Molecule::build_rotdrag_params");
06303     }
06304     
06305     if (tPDB->num_atoms() != numAtoms) {
06306       NAMD_die("Number of atoms in rotating drag tag PDB doesn't match coordinate PDB");
06307     }
06308   }
06309   
06310   // Get the PDB to read atom rotation axes. If no name given, use
06311   // rotDragFile if both it AND rotDragPivotFile are defined. Can NOT
06312   // use the PDB coordinate file, nor rotDragPivotFile!
06313 
06314   if (rotDragAxisFile == NULL) {
06315     if (rotDragFile == NULL) {
06316       NAMD_die("Rotating drag axis file can not be same as coordinate PDB file");
06317     } else {
06318       if (rotDragAxisFile->next != NULL) {
06319         NAMD_die("Multiple definitions of rotating drag axis file in configuration file");
06320       };
06321       if (rotDragPivotFile == NULL) {
06322         NAMD_die("Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
06323       };
06324       aPDB = tPDB;
06325     };
06326 
06327   } else {
06328 
06329     if ( (cwd == NULL) || (rotDragAxisFile->data[0] == '/') ) {
06330       strcpy(axisfilename, rotDragAxisFile->data);
06331     } else {
06332       strcpy(axisfilename, cwd);
06333       strcat(axisfilename, rotDragAxisFile->data);
06334     }
06335     
06336     aPDB = new PDB(axisfilename);
06337     if ( aPDB == NULL ) {
06338       NAMD_die("Memory allocation failed in Molecule::build_rotdrag_params");
06339     }
06340     
06341     if (aPDB->num_atoms() != numAtoms) {
06342       NAMD_die("Number of atoms in rotating drag axis PDB doesn't match coordinate PDB");
06343     }
06344   };
06345   
06346   // Get the PDB to read atom rotation pivot points. If no name given,
06347   // use rotDragFile if both it AND rotDragAxisFile are defined. Can
06348   // NOT use the PDB coordinate file, nor rotDragAxisFile!
06349 
06350   if (rotDragPivotFile == NULL) {
06351     if (rotDragFile == NULL) {
06352       NAMD_die("Rotating drag pivot point file can not be same as coordinate PDB file");
06353     } else {
06354       if (rotDragPivotFile->next != NULL) {
06355         NAMD_die("Multiple definitions of rotating drag pivot point file in configuration file");
06356       };
06357       if (rotDragAxisFile == NULL) {
06358         NAMD_die("Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
06359       };
06360       pPDB = tPDB;
06361     };
06362 
06363   } else {
06364 
06365     if ( (cwd == NULL) || (rotDragPivotFile->data[0] == '/') ) {
06366       strcpy(pivotfilename, rotDragPivotFile->data);
06367     } else {
06368       strcpy(pivotfilename, cwd);
06369       strcat(pivotfilename, rotDragPivotFile->data);
06370     }
06371     
06372     pPDB = new PDB(pivotfilename);
06373     if ( pPDB == NULL ) {
06374       NAMD_die("Memory allocation failed in Molecule::build_rotdrag_params");
06375     }
06376     
06377     if (pPDB->num_atoms() != numAtoms) {
06378       NAMD_die("Number of atoms in rotating drag pivot point PDB doesn't match coordinate PDB");
06379     }
06380   };
06381   
06382   
06383   // Get the PDB to read atom angular velocities. If no name given,
06384   // use rotDragFile (or the coordinate PDB file if rotDragFile is not
06385   // defined).
06386 
06387   if (rotDragVelFile == NULL) {
06388     vPDB = tPDB;
06389   } else {
06390     if (rotDragVelFile->next != NULL) {
06391       NAMD_die("Multiple definitions of rotating drag velocity file in configuration file");
06392     };
06393     
06394     if ( (cwd == NULL) || (rotDragVelFile->data[0] == '/') ) {
06395       strcpy(velfilename, rotDragVelFile->data);
06396     } else {
06397       strcpy(velfilename, cwd);
06398       strcat(velfilename, rotDragVelFile->data);
06399     }
06400     
06401     vPDB = new PDB(velfilename);
06402     if ( vPDB == NULL ) {
06403       NAMD_die("Memory allocation failed in Molecule::build_rotdrag_params");
06404     }
06405     
06406     if (vPDB->num_atoms() != numAtoms) {
06407       NAMD_die("Number of atoms in rotating drag velocity PDB doesn't match coordinate PDB");
06408     }
06409   };
06410   
06411   //  Get the column that the drag tag is going to be in. If
06412   //  rotDragFile is defined, it can be in any of the 5 floating point
06413   //  fields in the PDB (X, Y, Z, O, or B) which correspond to the
06414   //  1st, 2nd, ... 5th floating point fields. If rotDragFile is NOT
06415   //  defined, it can only be O or B fileds. The default is the O
06416   //  (4th) field, which is the occupancy.
06417 
06418   if (rotDragCol == NULL) {
06419     dtcol = 4;
06420   } else {
06421     if (rotDragCol->next != NULL) {
06422       NAMD_die("Multiple definitions of drag tag column in config file");
06423     };
06424     
06425     if ( rotDragFile == NULL
06426          && (!strcasecmp(rotDragCol->data, "X")
06427              || !strcasecmp(rotDragCol->data, "Y")
06428              || !strcasecmp(rotDragCol->data, "Z"))) {
06429       NAMD_die("Can not read rotating drag tags from X, Y, or Z column of the PDB coordinate file");
06430     };
06431     if (!strcasecmp(rotDragCol->data, "X")) {
06432       dtcol=1;
06433     } else if (!strcasecmp(rotDragCol->data, "Y")) {
06434       dtcol=2;
06435     } else if (!strcasecmp(rotDragCol->data, "Z")) {
06436       dtcol=3;
06437     } else if (!strcasecmp(rotDragCol->data, "O")) {
06438       dtcol=4;
06439     } else if (!strcasecmp(rotDragCol->data, "B")) {
06440       dtcol=5;
06441     }
06442     else {
06443       NAMD_die("rotDragCol must have value of X, Y, Z, O, or B");
06444     };
06445   };
06446   
06447   //  Get the column that the drag angular velocity is going to be
06448   //  in. If rotDragVelFile is defined, it can be in any of the 5
06449   //  floating point fields in the PDB (X, Y, Z, O, or B) which
06450   //  correspond to the 1st, 2nd, ... 5th floating point fields. If
06451   //  NEITHER of rotDragVelFile OR rotDragFile is defined, it can
06452   //  only be O or B fileds. The default is the O (4th) field, which
06453   //  is the occupancy.
06454 
06455   if (rotDragVelCol == NULL) {
06456     dvcol = 4;
06457   } else {
06458     if (rotDragVelCol->next != NULL) {
06459       NAMD_die("Multiple definitions of drag angular velocity column in config file");
06460     };
06461     
06462     if (rotDragVelFile == NULL
06463         && rotDragFile == NULL
06464         && strcasecmp(rotDragCol->data, "B")
06465         && strcasecmp(rotDragCol->data, "O")) {
06466       NAMD_die("Can not read rotating drag angular velocities from X, Y, or Z column of the PDB coordinate file");
06467     };
06468     if (!strcasecmp(rotDragVelCol->data, "X")) {
06469       dvcol=1;
06470     } else if (!strcasecmp(rotDragVelCol->data, "Y")) {
06471       dvcol=2;
06472     } else if (!strcasecmp(rotDragVelCol->data, "Z")) {
06473       dvcol=3;
06474     } else if (!strcasecmp(rotDragVelCol->data, "O")) {
06475       dvcol=4;
06476     } else if (!strcasecmp(rotDragVelCol->data, "B")) {
06477       dvcol=5;
06478     }
06479     else {
06480       NAMD_die("rotDragVelCol must have value of X, Y, Z, O, or B");
06481     };
06482   };
06483   
06484   //  Allocate an array that will store an index into the drag
06485   //  parameters for each atom.  If the atom is not dragged, its
06486   //  value will be set to -1 in this array.
06487   rotDragIndexes = new int32[numAtoms];
06488   if (rotDragIndexes == NULL) {
06489       NAMD_die("memory allocation failed in Molecule::build_rotdrag_params()");
06490   };
06491   
06492   //  Loop through all the atoms and find out which ones are dragged
06493   for (i=0; i<numAtoms; i++) {
06494     switch (dtcol) {
06495     case 1:
06496       dtval = (tPDB->atom(i))->xcoor();
06497       break;
06498     case 2:
06499       dtval = (tPDB->atom(i))->ycoor();
06500       break;
06501     case 3:
06502       dtval = (tPDB->atom(i))->zcoor();
06503       break;
06504     case 4:
06505       dtval = (tPDB->atom(i))->occupancy();
06506       break;
06507     case 5:
06508       dtval = (tPDB->atom(i))->temperaturefactor();
06509       break;
06510     }
06511     
06512     if (dtval != 0.0) {
06513       //  This atom is dragged
06514       rotDragIndexes[i] = current_index;
06515       current_index++;
06516     } else {
06517       //  This atom is not dragged
06518       rotDragIndexes[i] = -1;
06519     }
06520   }
06521   
06522   if (current_index == 0) {
06523     iout << iWARN << "NO DRAGGED ATOMS WERE FOUND, BUT ROTATING DRAG IS ON . . . " << endi;
06524   } else {
06525     rotDragParams = new RotDragParams[current_index];
06526     if (rotDragParams == NULL) {
06527       NAMD_die("memory allocation failed in Molecule::build_rotdrag_params");
06528     }
06529   };
06530   
06531   numRotDrag = current_index;
06532   
06533   //  Loop through all the atoms and assign the parameters for those
06534   //  that are dragged
06535   for (i=0; i<numAtoms; i++) {
06536     if (rotDragIndexes[i] != -1) {
06537       rotDragParams[rotDragIndexes[i]].a[0] = (aPDB->atom(i))->xcoor();
06538       rotDragParams[rotDragIndexes[i]].a[1] = (aPDB->atom(i))->ycoor();
06539       rotDragParams[rotDragIndexes[i]].a[2] = (aPDB->atom(i))->zcoor();
06540       rotDragParams[rotDragIndexes[i]].p[0] = (pPDB->atom(i))->xcoor();
06541       rotDragParams[rotDragIndexes[i]].p[1] = (pPDB->atom(i))->ycoor();
06542       rotDragParams[rotDragIndexes[i]].p[2] = (pPDB->atom(i))->zcoor();
06543       switch (dvcol) {
06544       case 1:
06545         rotDragParams[rotDragIndexes[i]].v = (vPDB->atom(i))->xcoor();
06546         break;
06547       case 2:
06548         rotDragParams[rotDragIndexes[i]].v = (vPDB->atom(i))->ycoor();
06549         break;
06550       case 3:
06551         rotDragParams[rotDragIndexes[i]].v = (vPDB->atom(i))->zcoor();
06552         break;
06553       case 4:
06554         rotDragParams[rotDragIndexes[i]].v = (vPDB->atom(i))->occupancy();
06555         break;
06556       case 5:
06557         rotDragParams[rotDragIndexes[i]].v = (vPDB->atom(i))->temperaturefactor();
06558         break;
06559       };
06560     };
06561   };
06562       
06563   if (rotDragFile != NULL) delete tPDB;
06564   if (rotDragAxisFile != NULL) delete aPDB;
06565   if (rotDragPivotFile != NULL) delete pPDB;
06566   if (rotDragVelFile != NULL) delete vPDB;
06567 }
06568 /*      END OF FUNCTION build_rotdrag_params    */
06569 
06570 
06571 /************************************************************************/
06572 /*                  */
06573 /*      FUNCTION build_constorque_params  */
06574 /*                  */
06575 /*   INPUTS:        */
06576 /*  consTorqueFile - value of consTorqueFile from the config file */
06577 /*  consTorqueCol - value of consTorqueCol from the config file */
06578 /*  consTorqueAxisFile - value of consTorqueAxisFile from the config file */
06579 /*  consTorquePivotFile - value of consTorquePivotFile from the config file */
06580 /*  consTorqueValFile - value of consTorqueValFile from the config file */
06581 /*  consTorqueValCol - value of consTorqueValCol from the config file */
06582 /*  initial_pdb - PDB object that contains initial positions  */
06583 /*  cwd - Current working directory          */
06584 /*                  */
06585 /*  This function builds all the parameters that are necessary  */
06586 /*  to do "constant" torque. This involves looking through one or more    */
06587 /*  PDB objects to determine which atoms are torqued,  and what the */
06588 /*  torque parameters for each atom are. This information is then stored */
06589 /*  in the arrays consTorqueIndexes and consTorqueParams. */
06590 /*                  */
06591 /************************************************************************/
06592 
06593 void Molecule::build_constorque_params(StringList *consTorqueFile, 
06594                                     StringList *consTorqueCol, 
06595                                     StringList *consTorqueAxisFile, 
06596                                     StringList *consTorquePivotFile, 
06597                                     StringList *consTorqueValFile, 
06598                                     StringList *consTorqueValCol, 
06599                                     PDB *initial_pdb,
06600                                     char *cwd)
06601   
06602 {
06603   PDB *tPDB, *aPDB, *pPDB, *vPDB; //  Pointers to other PDB file(s)
06604   register int i;          //  Loop counter
06605   int current_index=0;     //  Index into values used
06606   int dtcol = 4;           //  Column to look for torque tag in
06607   Real dtval = 0;          //  Torque tag value retreived
06608   int dvcol = 4;           //  Column to look for angular velocity in
06609   Real dvval = 0;          //  Angular velocity value retreived
06610   char mainfilename[129];  //  main "constant" torque PDB filename
06611   char axisfilename[129];  //  "constant" torque axis PDB filename
06612   char pivotfilename[129]; //  "constant" torque pivot point PDB filename
06613   char velfilename[129];   //  "constant" torque angular velocity PDB filename
06614   
06615   //  Get the PDB to read the "constant" torque tags from. Again, if the
06616   //  user gave us another file name, open that one.  Otherwise, just
06617   //  use the PDB with the initial coordinates
06618   if (consTorqueFile == NULL) {
06619     tPDB = initial_pdb;
06620     
06621   } else {
06622 
06623     if (consTorqueFile->next != NULL) {
06624       NAMD_die("Multiple definitions of \"constant\" torque tag file in configuration file");
06625     }
06626     
06627     if ( (cwd == NULL) || (consTorqueFile->data[0] == '/') ) {
06628       strcpy(mainfilename, consTorqueFile->data);
06629     } else {
06630       strcpy(mainfilename, cwd);
06631       strcat(mainfilename, consTorqueFile->data);
06632       }
06633     
06634     tPDB = new PDB(mainfilename);
06635     if ( tPDB == NULL ) {
06636       NAMD_die("Memory allocation failed in Molecule::build_constorque_params");
06637     }
06638     
06639     if (tPDB->num_atoms() != numAtoms) {
06640       NAMD_die("Number of atoms in \"constant\" torque tag PDB doesn't match coordinate PDB");
06641     }
06642   }
06643   
06644   // Get the PDB to read atom rotation axes. If no name given, use
06645   // consTorqueFile if both it AND consTorquePivotFile are defined. Can NOT
06646   // use the PDB coordinate file, nor consTorquePivotFile!
06647 
06648   if (consTorqueAxisFile == NULL) {
06649     if (consTorqueFile == NULL) {
06650       NAMD_die("\"Constant\" torque axis file can not be same as coordinate PDB file");
06651     } else {
06652       if (consTorqueAxisFile->next != NULL) {
06653         NAMD_die("Multiple definitions of \"constant\" torque axis file in configuration file");
06654       };
06655       if (consTorquePivotFile == NULL) {
06656         NAMD_die("Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
06657       };
06658       aPDB = tPDB;
06659     };
06660 
06661   } else {
06662 
06663     if ( (cwd == NULL) || (consTorqueAxisFile->data[0] == '/') ) {
06664       strcpy(axisfilename, consTorqueAxisFile->data);
06665     } else {
06666       strcpy(axisfilename, cwd);
06667       strcat(axisfilename, consTorqueAxisFile->data);
06668     }
06669     
06670     aPDB = new PDB(axisfilename);
06671     if ( aPDB == NULL ) {
06672       NAMD_die("Memory allocation failed in Molecule::build_constorque_params");
06673     }
06674     
06675     if (aPDB->num_atoms() != numAtoms) {
06676       NAMD_die("Number of atoms in \"constant\" torque axis PDB doesn't match coordinate PDB");
06677     }
06678   };
06679   
06680   // Get the PDB to read atom rotation pivot points. If no name given,
06681   // use consTorqueFile if both it AND consTorqueAxisFile are defined. Can
06682   // NOT use the PDB coordinate file, nor consTorqueAxisFile!
06683 
06684   if (consTorquePivotFile == NULL) {
06685     if (consTorqueFile == NULL) {
06686       NAMD_die("\"Constant\" torque pivot point file can not be same as coordinate PDB file");
06687     } else {
06688