CompressPsf.C

Go to the documentation of this file.
00001 #include <algorithm>
00002 #include "CompressPsf.h"
00003 #include "strlib.h"
00004 #include "Molecule.h"
00005 #include "Parameters.h"
00006 #include "SimParameters.h"
00007 #include "InfoStream.h"
00008 #include "UniqueSet.h"
00009 #include "UniqueSetIter.h"
00010 
00011 using namespace std;
00012 
00031 //global variables recording compressed psf file information
00032 Molecule *g_mol = NULL;
00033 Parameters *g_param = NULL;
00034 SimParameters *g_simParam = NULL;
00035 ConfigList *g_cfgList = NULL;
00036 
00037 struct BasicAtomInfo
00038 {
00039     Index segNameIdx;
00040     Index resNameIdx;
00041     Index atomNameIdx;
00042     Index atomTypeIdx;
00043     Index chargeIdx;
00044     Index massIdx;    
00045     Index atomSigIdx;
00046     Index exclSigIdx;
00047     int resID;
00048 };
00049 
00050 void OutputAtomRecord::flip(){
00051     flipNum((char *)&sSet, sizeof(short), sizeof(sSet)/sizeof(short));
00052     flipNum((char *)&iSet, sizeof(int), sizeof(iSet)/sizeof(int));
00053     flipNum((char *)&fSet, sizeof(float), sizeof(fSet)/sizeof(float));
00054 }
00055 
00056 struct AtomSigInfo
00057 {
00058     vector<SigIndex> bondSigIndices;
00059     vector<SigIndex> angleSigIndices;
00060     vector<SigIndex> dihedralSigIndices;
00061     vector<SigIndex> improperSigIndices;
00062     vector<SigIndex> crosstermSigIndices;
00063 
00064     AtomSigInfo()
00065     {}
00066     AtomSigInfo(const AtomSigInfo& sig)
00067     {
00068         bondSigIndices.clear();
00069         for(int i=0; i<sig.bondSigIndices.size(); i++)
00070             bondSigIndices.push_back(sig.bondSigIndices[i]);
00071 
00072         angleSigIndices.clear();
00073         for(int i=0; i<sig.angleSigIndices.size(); i++)
00074             angleSigIndices.push_back(sig.angleSigIndices[i]);
00075 
00076         dihedralSigIndices.clear();
00077         for(int i=0; i<sig.dihedralSigIndices.size(); i++)
00078             dihedralSigIndices.push_back(sig.dihedralSigIndices[i]);
00079 
00080         improperSigIndices.clear();
00081         for(int i=0; i<sig.improperSigIndices.size(); i++)
00082             improperSigIndices.push_back(sig.improperSigIndices[i]);
00083 
00084         crosstermSigIndices.clear();
00085         for(int i=0; i<sig.crosstermSigIndices.size(); i++)
00086             crosstermSigIndices.push_back(sig.crosstermSigIndices[i]);
00087     }
00088 
00089     ~AtomSigInfo()
00090     {
00091         bondSigIndices.clear();
00092         angleSigIndices.clear();
00093         dihedralSigIndices.clear();
00094         improperSigIndices.clear();
00095         crosstermSigIndices.clear();
00096     }
00097 
00098     void sortTupleSigIndices()
00099     {
00100         sort(bondSigIndices.begin(), bondSigIndices.end());
00101         sort(angleSigIndices.begin(), angleSigIndices.end());
00102         sort(dihedralSigIndices.begin(), dihedralSigIndices.end());
00103         sort(improperSigIndices.begin(), improperSigIndices.end());
00104         sort(crosstermSigIndices.begin(), crosstermSigIndices.end());
00105     }
00106   
00107     inline CkHashCode hash() const {
00108       // What's a good hash function for this? Lets make something up
00109       // Concatenate all the index lists into a list of chars, then hash that
00110       // string using Charm's string hash function
00111 
00112       // To keep memory allocation cheap, we'll just use a 32-byte buffer
00113       // and wrap around if we have more sigs
00114       const int maxlen = 32;
00115       unsigned char keydata[maxlen+1];
00116       const int maxchar = 256;
00117       int i,j;
00118       for(j=0;j<=maxlen;j++) keydata[j] = 0;
00119       j=0;
00120       for(i=0; i<bondSigIndices.size(); i++,j++) {
00121         keydata[j % maxlen] ^= (bondSigIndices[i] % maxchar);
00122       }
00123       for(i=0; i<angleSigIndices.size(); i++,j++) {
00124         keydata[j % maxlen] ^= (angleSigIndices[i] % maxchar);
00125       }
00126       for(i=0; i<dihedralSigIndices.size(); i++,j++) {
00127         keydata[j % maxlen] ^= (dihedralSigIndices[i] % maxchar);
00128       }
00129       for(i=0; i<improperSigIndices.size(); i++,j++) {
00130         keydata[j % maxlen] ^= (improperSigIndices[i] % maxchar);
00131       }
00132       for(i=0; i<crosstermSigIndices.size(); i++,j++) {
00133         keydata[j % maxlen] ^= (crosstermSigIndices[i] % maxchar);
00134       }
00135 //      CmiPrintf("Computed hash string len %d,%d\n",j,maxlen);
00136       if (j > maxlen) j = maxlen;
00137 //      for(i=0; i < j; i++) {
00138 //        if (keydata[i] == 0)
00139 //          keydata[i] = 255;
00140 //        CmiPrintf("key[%d]=%d %p\n",i,keydata[i],keydata);
00141 //      }
00142       return CkHashFunction_default((const void*)keydata,(size_t)j);
00143     }
00144 };
00145 
00146 int operator==(const AtomSigInfo &s1, const AtomSigInfo& s2)
00147 {
00148     if(s1.bondSigIndices.size() != s2.bondSigIndices.size())
00149         return 0;
00150     if(s1.angleSigIndices.size() != s2.angleSigIndices.size())
00151         return 0;
00152     if(s1.dihedralSigIndices.size() != s2.dihedralSigIndices.size())
00153         return 0;
00154     if(s1.improperSigIndices.size() != s2.improperSigIndices.size())
00155         return 0;
00156     if(s1.crosstermSigIndices.size() != s2.crosstermSigIndices.size())
00157         return 0;
00158 
00159     int equalCnt;
00160     equalCnt=0;
00161     int bondSigCnt = s1.bondSigIndices.size();
00162     for(int i=0; i<bondSigCnt; i++)
00163         equalCnt += (s1.bondSigIndices[i]==s2.bondSigIndices[i]);
00164     if(equalCnt!=bondSigCnt)
00165         return 0;
00166 
00167     equalCnt=0;
00168     int angleSigCnt = s1.angleSigIndices.size();
00169     for(int i=0; i<angleSigCnt; i++)
00170         equalCnt += (s1.angleSigIndices[i]==s2.angleSigIndices[i]);
00171     if(equalCnt!=angleSigCnt)
00172         return 0;
00173 
00174     equalCnt=0;
00175     int dihedralSigCnt = s1.dihedralSigIndices.size();
00176     for(int i=0; i<dihedralSigCnt; i++)
00177         equalCnt += (s1.dihedralSigIndices[i]==s2.dihedralSigIndices[i]);
00178     if(equalCnt!=dihedralSigCnt)
00179         return 0;
00180 
00181     equalCnt=0;
00182     int improperSigCnt = s1.improperSigIndices.size();
00183     for(int i=0; i<improperSigCnt; i++)
00184         equalCnt += (s1.improperSigIndices[i]==s2.improperSigIndices[i]);
00185     if(equalCnt!=improperSigCnt)
00186         return 0;
00187 
00188     equalCnt=0;
00189     int crosstermSigCnt = s1.crosstermSigIndices.size();
00190     for(int i=0; i<crosstermSigCnt; i++)
00191         equalCnt += (s1.crosstermSigIndices[i]==s2.crosstermSigIndices[i]);
00192     if(equalCnt!=crosstermSigCnt)
00193         return 0;
00194 
00195     return 1;
00196 }
00197 
00198 struct ExclSigInfo
00199 {
00200     vector<int> fullExclOffset;
00201     vector<int> modExclOffset;
00202 
00203     ExclSigInfo()
00204     {}
00205     ExclSigInfo(const ExclSigInfo& sig)
00206     {
00207         fullExclOffset.clear();
00208         for(int i=0; i<sig.fullExclOffset.size(); i++)
00209             fullExclOffset.push_back(sig.fullExclOffset[i]);
00210 
00211         modExclOffset.clear();
00212         for(int i=0; i<sig.modExclOffset.size(); i++)
00213             modExclOffset.push_back(sig.modExclOffset[i]);
00214     }
00215 
00216     ~ExclSigInfo()
00217     {
00218         fullExclOffset.clear();
00219         modExclOffset.clear();
00220     }
00221 
00222     void sortExclOffset()
00223     {
00224         sort(fullExclOffset.begin(), fullExclOffset.end());
00225         sort(modExclOffset.begin(), modExclOffset.end());
00226     }
00227 
00228     int hash() const {
00229       unsigned int code = 0x1234;
00230       unsigned int codesz = 8 * sizeof(int);
00231       const unsigned int numFoffset = fullExclOffset.size();
00232       const unsigned int numMoffset = modExclOffset.size();
00233       const unsigned int numOffsets = numFoffset + numMoffset;
00234       
00235       // No excluded atoms? Just hash to 0
00236       if (numOffsets == 0)
00237         return 0;
00238       
00239       unsigned int shift = codesz / numOffsets;
00240       if (shift == 0) shift=1;
00241       unsigned int i;
00242       for(i=0; i < numFoffset; i++) {
00243         code = circShift(code,shift);
00244         code ^= fullExclOffset[i];
00245       }
00246       for(i=0; i < numMoffset; i++) {
00247         code = circShift(code,shift);
00248         code ^= modExclOffset[i];
00249       }
00250       return code;
00251     }
00252 };
00253 int operator==(const ExclSigInfo &s1, const ExclSigInfo &s2)
00254 {
00255     if(s1.fullExclOffset.size()!=s2.fullExclOffset.size())
00256         return 0;
00257     if(s1.modExclOffset.size()!=s2.modExclOffset.size())
00258         return 0;
00259 
00260     for(int i=0; i<s1.fullExclOffset.size(); i++)
00261     {
00262         if(s1.fullExclOffset[i] != s2.fullExclOffset[i])
00263             return 0;
00264     }
00265 
00266     for(int i=0; i<s1.modExclOffset.size(); i++)
00267     {
00268         if(s1.modExclOffset[i] != s2.modExclOffset[i])
00269             return 0;
00270     }
00271     return 1;
00272 }
00273 
00274 class HashString : public string {
00275 public:
00276   int hash() const {
00277     const char* d = this->c_str();
00278     int ret=0;
00279     for (int i=0;d[i]!=0;i++) {
00280       int shift1=((5*i)%16)+0;
00281       int shift2=((6*i)%16)+8;
00282       ret+=((0xa5^d[i])<<shift2)+(d[i]<<shift1);
00283     }
00284     return ret;
00285   }
00286 };
00287 
00288 class HashReal {
00289   Real val;
00290 public:
00291   HashReal(Real v) : val(v) {}
00292   operator Real & () { return val; }
00293   operator const Real & () const { return val; }
00294   
00295   int hash() const {
00296     const char* d = (const char *)&val;
00297     int ret=0;
00298     for (int i=0;i < sizeof(Real);i++) {
00299       int shift1=((5*i)%16)+0;
00300       int shift2=((6*i)%16)+8;
00301       ret+=((0xa5^d[i])<<shift2)+(d[i]<<shift1);
00302     }
00303     return ret;
00304   };
00305 };
00306 
00307 HashPool<HashString> segNamePool;
00308 HashPool<HashString> resNamePool;
00309 HashPool<HashString> atomNamePool;
00310 HashPool<HashString> atomTypePool;
00311 HashPool<HashReal> chargePool;
00312 HashPool<HashReal> massPool;
00313 HashPool<AtomSigInfo> atomSigPool;
00314 BasicAtomInfo *atomData;
00315 
00316 //Recording cluster information after reading all bonds info
00317 int *eachAtomClusterID = NULL;
00318 vector<int> eachClusterSize;
00319 vector<int> eachClusterID;
00320 int g_numClusters = 0;
00321 
00322 HashPool<TupleSignature> sigsOfBonds;
00323 HashPool<TupleSignature> sigsOfAngles;
00324 HashPool<TupleSignature> sigsOfDihedrals;
00325 HashPool<TupleSignature> sigsOfImpropers;
00326 HashPool<TupleSignature> sigsOfCrossterms;
00327 AtomSigInfo *eachAtomSigs;
00328 
00329 HashPool<ExclSigInfo> sigsOfExclusions;
00330 ExclSigInfo *eachAtomExclSigs;
00331 
00332 //Structures for extraBond options
00333 vector<Bond> extraBonds;
00334 vector<Angle> extraAngles;
00335 vector<Dihedral> extraDihedrals;
00336 vector<Improper> extraImpropers;
00337 
00338 vector<BondValue> extraBondParams;
00339 vector<AngleValue> extraAngleParams;
00340 vector<DihedralValue> extraDihedralParams;
00341 vector<ImproperValue> extraImproperParams;
00342 
00343 int operator==(const BondValue &b1, const BondValue &b2)
00344 {
00345     return (b1.k==b2.k) && (b1.x0==b2.x0);
00346 }
00347 
00348 int operator==(const AngleValue &a1, const AngleValue &a2)
00349 {
00350     return (a1.k==a2.k) && (a1.k_ub==a2.k_ub) && (a1.r_ub==a2.r_ub) && (a1.theta0==a2.theta0);
00351 }
00352 
00353 int operator!=(const FourBodyConsts& f1, const FourBodyConsts& f2)
00354 {
00355     return (f1.delta!=f2.delta) || (f1.k!=f2.k) || (f1.n!=f2.n);
00356 }
00357 
00358 int operator==(const DihedralValue &d1, const DihedralValue &d2)
00359 {
00360     if(d1.multiplicity != d2.multiplicity)
00361         return 0;
00362     for(int i=0; i<MAX_MULTIPLICITY; i++)
00363     {
00364         if(d1.values[i] != d2.values[i])
00365             return 0;
00366     }
00367     return 1;
00368 }
00369 
00370 int operator==(const ImproperValue &d1, const ImproperValue &d2)
00371 {
00372     if(d1.multiplicity != d2.multiplicity)
00373         return 0;
00374     for(int i=0; i<MAX_MULTIPLICITY; i++)
00375     {
00376         if(d1.values[i] != d2.values[i])
00377             return 0;
00378     }
00379     return 1;
00380 }
00381 
00382 void loadMolInfo();
00383 
00384 void integrateAllAtomSigs();
00385 void outputCompressedFile(FILE *txtOfp, FILE *binOfp);
00386 
00387 //reading extraBond's information
00388 void getExtraBonds(StringList *file);
00389 
00390 void buildAtomData();
00391 void buildBondData();
00392 void buildAngleData();
00393 void buildDihedralData();
00394 void buildImproperData();
00395 void buildCrosstermData();
00396 
00397 void buildExclusionData();
00398 
00399 //Functions related with building exclusions
00400 void buildExclusions();
00401 void build12Excls(UniqueSet<Exclusion>&, vector<int> *);
00402 void build13Excls(UniqueSet<Exclusion>&, vector<int> *);
00403 void build14Excls(UniqueSet<Exclusion>&, vector<int> *, int);
00404 
00405 //reverse the byte-order of every element starting at "elem"
00406 void flipNum(char *elem, int elemSize, int numElems){
00407     int mid = elemSize/2;
00408     char *ptr = elem;
00409     for(int i=0; i<numElems; i++) {
00410         for(int j=0; j<mid; j++) {
00411             char tmp = ptr[j];
00412             ptr[j] = ptr[elemSize-1-j];
00413             ptr[elemSize-1-j] = tmp;
00414         }
00415         ptr += elemSize;
00416     }
00417 }
00418 
00419 void clearGlobalVectors()
00420 {
00421     segNamePool.clear();
00422     resNamePool.clear();
00423     atomNamePool.clear();
00424     atomTypePool.clear();
00425     chargePool.clear();
00426     massPool.clear();
00427     sigsOfBonds.clear();
00428     sigsOfAngles.clear();
00429     sigsOfDihedrals.clear();
00430     sigsOfImpropers.clear();
00431     sigsOfExclusions.clear();
00432 
00433     eachClusterSize.clear();
00434 }
00435 
00436 void compress_molecule_info(Molecule *mol, char *psfFileName, Parameters *param, SimParameters *simParam, ConfigList* cfgList)
00437 {
00438     g_mol = mol;
00439     g_param = param;
00440     g_simParam = simParam; //used for building exclusions
00441     g_cfgList = cfgList; //used for integrating extra bonds
00442 
00443     //read psf files
00444     //readPsfFile(psfFileName);
00445     loadMolInfo();
00446 
00447     integrateAllAtomSigs();
00448 
00449     buildExclusions();
00450 
00451     //buildParamData();
00452 
00453     char *outFileName = new char[strlen(psfFileName)+20];
00454     sprintf(outFileName, "%s.inter", psfFileName);
00455     //the text file for signatures and other non-per-atom info
00456     FILE *txtOfp = fopen(outFileName, "w");
00457     sprintf(outFileName, "%s.inter.bin", psfFileName);
00458     //the binary file for per-atom info
00459     FILE *binOfp = fopen(outFileName, "wb");
00460     delete [] outFileName;
00461 
00462     //output compressed psf file
00463     outputCompressedFile(txtOfp, binOfp);
00464 
00465     fclose(txtOfp);
00466     fclose(binOfp);
00467 }
00468 
00473 void loadMolInfo()
00474 {
00475     char err_msg[512];  //  Error message for NAMD_die
00476     char buffer[512];  //  Buffer for file reading
00477     int i;      //  Loop counter
00478     int NumTitle;    //  Number of Title lines in .psf file
00479     int ret_code;    //  ret_code from NAMD_read_line calls
00480     FILE *psf_file;
00481     Parameters *params = g_param;
00482 
00483     buildAtomData();
00484 
00485     //read extra bonds/angles/dihedrals/impropers information first
00486     //and then integrate them into the following reading procedures
00487     /*if(g_simParam->extraBondsOn)
00488     {
00489         getExtraBonds(g_cfgList->find("extraBondsFile"));
00490     }*/
00491 
00492     //initialize eachAtomSigs
00493     eachAtomSigs = new AtomSigInfo[g_mol->numAtoms];    
00494 
00495     buildBondData();
00496     buildAngleData();
00497     buildDihedralData();
00498     buildImproperData();
00499     buildCrosstermData();
00500 
00501     //  analyze the data and find the status of each atom
00502     //build_atom_status();
00503 }
00504 
00505 void integrateAllAtomSigs()
00506 {
00507     printf("Bond sigs:  %d\n", (int)sigsOfBonds.size());
00508     printf("Angle sigs:  %d\n", (int)sigsOfAngles.size());
00509     printf("Dihedral sigs:  %d\n", (int)sigsOfDihedrals.size());
00510     printf("Improper sigs:  %d\n", (int)sigsOfImpropers.size());
00511     printf("Crossterm sigs:  %d\n", (int)sigsOfCrossterms.size());
00512 
00513 
00514     for(int i=0; i<g_mol->numAtoms; i++)
00515     {
00516         eachAtomSigs[i].sortTupleSigIndices();
00517         int poolIndex = atomSigPool.lookupCstPool(eachAtomSigs[i]);
00518         if(poolIndex==-1)
00519         {
00520             atomSigPool.push_back(eachAtomSigs[i]);
00521             poolIndex = atomSigPool.size()-1;
00522         }
00523         atomData[i].atomSigIdx = poolIndex;
00524     }
00525 
00526     printf("Atom's sigs: %d\n", (int)atomSigPool.size());
00527 
00528     delete[] eachAtomSigs;
00529 }
00530 
00541 void outputCompressedFile(FILE *txtOfp, FILE *binOfp)
00542 {
00543 #ifndef MEM_OPT_VERSION
00544     fprintf(txtOfp, "FORMAT VERSION: %f\n", COMPRESSED_PSF_VER);
00545 
00546     fprintf(txtOfp, "%d !NSEGMENTNAMES\n", segNamePool.size());
00547     for(int i=0; i<segNamePool.size(); i++)
00548     {
00549         fprintf(txtOfp, "%s\n", segNamePool[i].c_str());
00550     }
00551 
00552     fprintf(txtOfp, "%d !NRESIDUENAMES\n", resNamePool.size());
00553     for(int i=0; i<resNamePool.size(); i++)
00554     {
00555         fprintf(txtOfp, "%s\n", resNamePool[i].c_str());
00556     }
00557 
00558     fprintf(txtOfp, "%d !NATOMNAMES\n", atomNamePool.size());
00559     for(int i=0; i<atomNamePool.size(); i++)
00560     {
00561         fprintf(txtOfp, "%s\n", atomNamePool[i].c_str());
00562     }
00563 
00564     fprintf(txtOfp, "%d !NATOMTYPES\n", atomTypePool.size());
00565     for(int i=0; i<atomTypePool.size(); i++)
00566     {
00567         fprintf(txtOfp, "%s\n", atomTypePool[i].c_str());
00568     }
00569 
00570     fprintf(txtOfp, "%d !NCHARGES\n", chargePool.size());
00571     for(int i=0; i<chargePool.size(); i++)
00572     {
00573         const Real charge = chargePool[i];
00574         fprintf(txtOfp, "%f\n", charge);
00575     }
00576 
00577     fprintf(txtOfp, "%d !NMASSES\n", massPool.size());
00578     for(int i=0; i<massPool.size(); i++)
00579     {
00580         const Real mass = massPool[i];
00581         fprintf(txtOfp, "%f\n", mass);
00582     }
00583 
00584 
00585     fprintf(txtOfp, "%d !NATOMSIGS\n", atomSigPool.size());
00586     for(int i=0; i<atomSigPool.size(); i++)
00587     {
00588         AtomSigInfo& oneAtomSig = atomSigPool[i];
00589         int oneTypeCnt = oneAtomSig.bondSigIndices.size();
00590         fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NBOND");
00591         for(int j=0; j<oneTypeCnt; j++)
00592         {
00593             SigIndex idx = oneAtomSig.bondSigIndices[j];
00594             TupleSignature& tSig = sigsOfBonds[idx];
00595             tSig.output(txtOfp);
00596         }
00597 
00598         oneTypeCnt = oneAtomSig.angleSigIndices.size();
00599         fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NTHETA");
00600         for(int j=0; j<oneTypeCnt; j++)
00601         {
00602             SigIndex idx = oneAtomSig.angleSigIndices[j];
00603             TupleSignature& tSig = sigsOfAngles[idx];
00604             tSig.output(txtOfp);
00605         }
00606 
00607         oneTypeCnt = oneAtomSig.dihedralSigIndices.size();
00608         fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NPHI");
00609         for(int j=0; j<oneTypeCnt; j++)
00610         {
00611             SigIndex idx = oneAtomSig.dihedralSigIndices[j];
00612             TupleSignature& tSig = sigsOfDihedrals[idx];
00613             tSig.output(txtOfp);
00614         }
00615 
00616         oneTypeCnt = oneAtomSig.improperSigIndices.size();
00617         fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NIMPHI");
00618         for(int j=0; j<oneTypeCnt; j++)
00619         {
00620             SigIndex idx = oneAtomSig.improperSigIndices[j];
00621             TupleSignature& tSig = sigsOfImpropers[idx];
00622             tSig.output(txtOfp);
00623         }
00624 
00625         oneTypeCnt = oneAtomSig.crosstermSigIndices.size();
00626         fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NCRTERM");
00627         for(int j=0; j<oneTypeCnt; j++)
00628         {
00629             SigIndex idx = oneAtomSig.crosstermSigIndices[j];
00630             TupleSignature& tSig = sigsOfCrossterms[idx];
00631             tSig.output(txtOfp);
00632         }
00633     }
00634 
00635     //2. Output exclusion signatures
00636     int exclSigCnt = sigsOfExclusions.size();
00637     fprintf(txtOfp, "%d !NEXCLSIGS\n", exclSigCnt);
00638     for(int i=0; i<exclSigCnt; i++)
00639     {
00640         ExclSigInfo *sig = &sigsOfExclusions[i];
00641         //first line is for full exclusions (1-2, 1-3) in the format of count offset1 offset2 offset3 ...
00642         fprintf(txtOfp, "%d", sig->fullExclOffset.size());
00643         for(int j=0; j<sig->fullExclOffset.size(); j++)
00644             fprintf(txtOfp, " %d", sig->fullExclOffset[j]);
00645         fprintf(txtOfp, "\n");
00646 
00647         //second line is for modified exclusions (1-4)
00648         fprintf(txtOfp, "%d", sig->modExclOffset.size());
00649         for(int j=0; j<sig->modExclOffset.size(); j++)
00650             fprintf(txtOfp, " %d", sig->modExclOffset[j]);
00651         fprintf(txtOfp, "\n");
00652     }
00653 
00654     //3. Output the cluster information
00655     fprintf(txtOfp, "%d !NCLUSTERS\n", g_numClusters);
00656 
00657     //4. Output atom info
00658     fprintf(txtOfp, "%d !NATOM\n", g_mol->numAtoms);
00659     fprintf(txtOfp, "%d !NHYDROGENGROUP\n", g_mol->numHydrogenGroups);
00660     fprintf(txtOfp, "%d !MAXHYDROGENGROUPSIZE\n", g_mol->maxHydrogenGroupSize);
00661     fprintf(txtOfp, "%d !NMIGRATIONGROUP\n", g_mol->numMigrationGroups);
00662     fprintf(txtOfp, "%d !MAXMIGRATIONGROUPSIZE\n", g_mol->maxMigrationGroupSize);
00663 
00664     //5. Output rigid bond type
00665     fprintf(txtOfp, "%d !RIGIDBONDTYPE\n", g_simParam->rigidBonds);
00666 #if 0
00667     const float *atomOccupancy = g_mol->getOccupancyData();
00668     const float *atomBFactor = g_mol->getBFactorData();
00669     fprintf(txtOfp, "%d !OCCUPANCYVALID\n", (atomOccupancy==NULL)?0:1);
00670     fprintf(txtOfp, "%d !TEMPFACTORVALID\n", (atomBFactor==NULL)?0:1);
00671 
00672     float *zeroFloats = NULL;
00673     if(atomOccupancy==NULL || atomBFactor==NULL) {
00674         zeroFloats = new float[g_mol->numAtoms];
00675         memset(zeroFloats, 0, sizeof(float)*g_mol->numAtoms);
00676         if(atomOccupancy==NULL) atomOccupancy = (const float *)zeroFloats;
00677         if(atomBFactor==NULL) atomBFactor = (const float *)zeroFloats;
00678     }
00679 #endif
00680 
00681     Atom *atoms = g_mol->getAtoms(); //need to output its partner and hydrogenList
00682     HydrogenGroupID *hg = g_mol->hydrogenGroup.begin();
00683 
00684     //First, output magic number
00685     int magicNum = COMPRESSED_PSF_MAGICNUM;
00686     fwrite(&magicNum, sizeof(int), 1, binOfp);
00687     //Second, version number
00688     float verNum = (float)COMPRESSED_PSF_VER;
00689     fwrite(&verNum, sizeof(float), 1, binOfp);
00690     //Third, the per-atom record size
00691     int recSize = sizeof(OutputAtomRecord);
00692     fwrite(&recSize, sizeof(int), 1, binOfp);
00693     //Fourth, each atom info
00694     OutputAtomRecord oneRec;
00695     for(int i=0; i<g_mol->numAtoms; i++)
00696     {                
00697         oneRec.sSet.segNameIdx = atomData[i].segNameIdx;
00698         oneRec.sSet.resNameIdx = atomData[i].resNameIdx;
00699         oneRec.sSet.atomNameIdx = atomData[i].atomNameIdx;
00700         oneRec.sSet.atomTypeIdx = atomData[i].atomTypeIdx;
00701         oneRec.sSet.chargeIdx = atomData[i].chargeIdx;
00702         oneRec.sSet.massIdx = atomData[i].massIdx;
00703         oneRec.iSet.atomSigIdx = atomData[i].atomSigIdx;
00704 
00705         oneRec.iSet.exclSigIdx = atomData[i].exclSigIdx;        
00706         oneRec.sSet.vdw_type = atoms[i].vdw_type;
00707         oneRec.iSet.resID = atomData[i].resID;        
00708         int hydIdx = atoms[i].hydrogenList;       
00709         oneRec.iSet.hydrogenList = hydIdx;
00710         oneRec.iSet.atomsInGroup = hg[hydIdx].atomsInGroup;
00711         oneRec.iSet.GPID = hg[hydIdx].GPID;
00712         //oneRec.waterVal = hg[hydIdx].waterVal;
00713         oneRec.iSet.atomsInMigrationGroup = hg[hydIdx].atomsInMigrationGroup;
00714         oneRec.iSet.MPID = hg[hydIdx].MPID;
00715         //oneRec.fSet.occupancy = atomOccupancy[i];
00716         //oneRec.fSet.bfactor = atomBFactor[i];
00717         oneRec.fSet.rigidBondLength = g_mol->rigid_bond_length(i);
00718         fwrite(&oneRec, sizeof(OutputAtomRecord), 1, binOfp);
00719     }
00720     //if(zeroFloats) delete[] zeroFloats;
00721     delete[] atomData;
00722 
00723     //Output the info required for parallel output:
00724     //1. Cluster ids of each atom
00725     //Since the cluster info is only when doing output under
00726     //wrapAll or wrapWater conditions, for the purpose of parallel
00727     //output, it is reasonable to separate the cluster info from
00728     //the above per-atom info. So whole the binary per-atom file
00729     //contains two parts, one for parallel input to create patches
00730     //and computes; the other for parallel output -Chao Mei
00731     //Fifth: output the cluster id of each atoms
00732     fwrite(eachAtomClusterID, sizeof(int), g_mol->numAtoms, binOfp);    
00733     //2. Whether each atom is water or not
00734     char *isWater = new char[g_mol->numAtoms];
00735     for(int i=0; i<g_mol->numAtoms; i++){
00736         isWater[i] = g_mol->is_water(i);
00737     }
00738     fwrite(isWater, sizeof(char), g_mol->numAtoms, binOfp);
00739     delete [] isWater;
00740     delete[] atoms;
00741     g_mol->hydrogenGroup.resize(0);
00742     delete[] eachAtomClusterID;    
00743 
00744     //Output the parameter new values if extraBonds are present.
00745     //The parameters are not needed since now extraBonds' parameters will be
00746     //read again during running the simulation
00747 
00748     //6. Output the "multiplicity" field TUPLE_array of the Parameter object
00749     fprintf(txtOfp, "!DIHEDRALPARAMARRAY\n");
00750     for(int i=0; i<g_param->NumDihedralParams; i++)
00751     {
00752         fprintf(txtOfp, "%d ", g_param->dihedral_array[i].multiplicity);
00753     }
00754     fprintf(txtOfp, "\n");
00755     fprintf(txtOfp, "!IMPROPERPARAMARRAY\n");
00756     for(int i=0; i<g_param->NumImproperParams; i++)
00757     {
00758         fprintf(txtOfp, "%d ", g_param->improper_array[i].multiplicity);
00759     }
00760     fprintf(txtOfp, "\n");
00761 #endif
00762 }
00763 
00764 void buildAtomData()
00765 {
00766 #ifndef MEM_OPT_VERSION
00767     int numAtoms = g_mol->numAtoms;
00768 
00769     //1. parse atom data to build constant pool (atom name, mass, charge etc.)
00770     atomData = new BasicAtomInfo[numAtoms];
00771     Atom *atoms = g_mol->getAtoms();
00772     AtomNameInfo *atomNames = g_mol->getAtomNames();
00773     AtomSegResInfo *atomSegResids = g_mol->getAtomSegResInfo();
00774 
00775     for(int atomID=0; atomID < numAtoms; atomID++)
00776     {
00777         //building constant pool
00778         int poolIndex;
00779         HashString fieldName;
00780         fieldName.assign(atomSegResids[atomID].segname);
00781         poolIndex = segNamePool.lookupCstPool(fieldName);
00782         if(poolIndex==-1)
00783         {
00784             segNamePool.push_back(fieldName);
00785             poolIndex = segNamePool.size()-1;
00786         }
00787         atomData[atomID].segNameIdx = poolIndex;
00788         
00789         atomData[atomID].resID = atomSegResids[atomID].resid;
00790 
00791         fieldName.assign(atomNames[atomID].resname);
00792         poolIndex = resNamePool.lookupCstPool(fieldName);
00793         if(poolIndex==-1)
00794         {
00795             resNamePool.push_back(fieldName);
00796             poolIndex = resNamePool.size()-1;
00797         }
00798         atomData[atomID].resNameIdx = poolIndex;
00799 
00800         fieldName.assign(atomNames[atomID].atomname);
00801         poolIndex = atomNamePool.lookupCstPool(fieldName);
00802         if(poolIndex==-1)
00803         {
00804             atomNamePool.push_back(fieldName);
00805             poolIndex = atomNamePool.size()-1;
00806         }
00807         atomData[atomID].atomNameIdx = poolIndex;
00808 
00809         fieldName.assign(atomNames[atomID].atomtype);
00810         poolIndex = atomTypePool.lookupCstPool(fieldName);
00811         if(poolIndex==-1)
00812         {
00813             atomTypePool.push_back(fieldName);
00814             poolIndex = atomTypePool.size()-1;
00815         }
00816         atomData[atomID].atomTypeIdx = poolIndex;
00817 
00818         poolIndex = chargePool.lookupCstPool(atoms[atomID].charge);
00819         if(poolIndex==-1)
00820         {
00821             chargePool.push_back(atoms[atomID].charge);
00822             poolIndex = chargePool.size()-1;
00823         }
00824         atomData[atomID].chargeIdx = poolIndex;
00825 
00826         poolIndex = massPool.lookupCstPool(atoms[atomID].mass);
00827         if(poolIndex==-1)
00828         {
00829             massPool.push_back(atoms[atomID].mass);
00830             poolIndex = massPool.size()-1;
00831         }
00832         atomData[atomID].massIdx = poolIndex;
00833     }
00834 
00835     //Free those space to reduce transient memory usage
00836     //delete [] atoms; (deleted until per-atom info is output)
00837     delete [] atomNames;
00838     delete [] atomSegResids;
00839 #endif
00840 }
00841 
00842 
00843 void buildBondData()
00844 {
00845 #ifndef MEM_OPT_VERSION
00846     Bond *bonds = g_mol->getAllBonds();    
00847 
00848     //then creating bond's tupleSignature
00849     for(int i=0; i<g_mol->numBonds; i++)
00850     {
00851         Bond *b = bonds+i;
00852         TupleSignature oneSig(1,BOND,b->bond_type);
00853         oneSig.offset[0] = b->atom2 - b->atom1;
00854         oneSig.isReal = (i<g_mol->numRealBonds);
00855 
00856         int poolIndex = sigsOfBonds.lookupCstPool(oneSig);
00857         int newSig=0;
00858         if(poolIndex == -1)
00859         {
00860             sigsOfBonds.push_back(oneSig);
00861             poolIndex = (SigIndex)sigsOfBonds.size()-1;
00862             newSig=1;
00863         }
00864 
00865       if(!newSig)
00866         {//check duplicate bonds in the form of (a, b) && (a, b);
00867           int dupIdx = lookupCstPool(eachAtomSigs[b->atom1].bondSigIndices, (SigIndex)poolIndex);
00868             if(dupIdx!=-1)
00869             {
00870                 char err_msg[128];
00871                 sprintf(err_msg, "Duplicate bond %d-%d!", b->atom1+1, b->atom2+1);
00872                 NAMD_die(err_msg);
00873             }
00874         }
00875         eachAtomSigs[b->atom1].bondSigIndices.push_back(poolIndex);
00876     }
00877 
00878     //check duplicate bonds in the form of (a, b) && (b, a)
00879     for(int i=0; i<g_mol->numBonds; i++)
00880     {
00881         Bond *b=bonds+i;
00882         int atom2 = b->atom2;
00883         int thisOffset = atom2 - b->atom1;
00884         for(int j=0; j<eachAtomSigs[atom2].bondSigIndices.size(); j++)
00885         {
00886             SigIndex atom2BondId = eachAtomSigs[atom2].bondSigIndices[j];
00887             TupleSignature *secSig = &(sigsOfBonds[atom2BondId]);
00888             if(thisOffset== -(secSig->offset[0]))
00889             {
00890                 char err_msg[128];
00891                 sprintf(err_msg, "Duplicate bond %d-%d because two atoms are just reversed!", b->atom1+1, atom2+1);
00892                 NAMD_die(err_msg);
00893             }
00894         }
00895     }      
00896 
00897     //building clusters for this simulation system in two steps
00898     //1. create a list for each atom where each atom in the list is bonded with that atom
00899     vector<int> *atomListOfBonded = new vector<int>[g_mol->numAtoms];
00900 
00901     for(int i=0; i<g_mol->numRealBonds; i++)
00902     {
00903         Bond *b=bonds+i;
00904         int atom1 = b->atom1;
00905         int atom2 = b->atom2;
00906         atomListOfBonded[atom1].push_back(atom2);
00907         atomListOfBonded[atom2].push_back(atom1);
00908     }
00909 
00910     delete [] bonds;
00911 
00912     //2. using breadth-first-search to build the clusters. Here, we avoid recursive call
00913     // because the depth of calls may be of thousands which will blow up the stack, and
00914     //recursive call is slower than the stack-based BFS.
00915     //Considering such structure
00916     //1->1245; 7->1243; 1243->1245
00917     eachAtomClusterID = new int[g_mol->numAtoms];
00918     for(int i=0; i<g_mol->numAtoms; i++)
00919         eachAtomClusterID[i] = -1;
00920 
00921     //It is guaranteed that the clusters found in this way use the
00922     //smallest atom id of this cluster because each atom is at least
00923     //connected to one of the atoms in its cluster (by atomListOfBonded
00924     //constructed above).
00925     //--Chao Mei
00926     for(int i=0; i<g_mol->numAtoms; i++)
00927     {
00928         int curClusterID=eachAtomClusterID[i];
00929         //if the atom's cluster id is not -1, the atom has been visited
00930         if(curClusterID!=-1) continue;
00931 
00932         curClusterID=i;
00933         deque<int> toVisitAtoms;
00934         eachAtomClusterID[i] = curClusterID;
00935         toVisitAtoms.push_back(i);
00936         while(!toVisitAtoms.empty())
00937         {
00938             int visAtomID = toVisitAtoms.front();
00939             toVisitAtoms.pop_front();
00940             for(int j=0; j<atomListOfBonded[visAtomID].size(); j++)
00941             {
00942                 int otherAtom = atomListOfBonded[visAtomID][j];
00943                 if(eachAtomClusterID[otherAtom]!=curClusterID){
00944                     eachAtomClusterID[otherAtom]=curClusterID;
00945                     toVisitAtoms.push_back(otherAtom);
00946                 }
00947             }
00948         }
00949     }
00950 
00951 #if 0
00952     //Now the clusterID of each atom should be usually in the non-decreasing
00953     //order. In other words, the atom ids of a cluster are generally contiguous.
00954     //If this is the case, the temporary memory usage of output IO during
00955     //the simulation can be dramatically reduced. So g_isClusterContiguous
00956     //is used to differentiate the two cases: 
00957 
00958     //1. if the cluster id of atoms is monotonically increasing 
00959     //(g_isClusterContiguous=1), the size of the cluster can be used as this
00960     //this cluster's signature (represented by "eachAtomClusterID").
00961     //The atom who is the first atom (in terms of atom id) in this cluster will 
00962     //store the cluster size as its signature. The remaining atoms in this 
00963     //cluster store -1.
00964     
00965     //2. if the cluster id of atoms is not monotonically increasing, that is,
00966     //the atom ids of a cluster are not contiguous (g_isClusterContiguous=0).
00967     //Then we have to still record each atom's cluster id in variable 
00968     //"eachAtomClusterID", and we have to change each unique cluster id (valued
00969     //at the scale of numAtoms) into indexes at the scale of numClusters. For
00970     //example, the cluster ids of atoms up to this point may be (0...0, 24...24,
00971     //789...789,...), after re-scaling, cluster ids should be (0...0,1...1,
00972     //2...2,....) for atoms.
00973     //We made an ASSUMPTION here: the cluster ids are in the non-decreasing
00974     //order, and when a cluster discontiguity happens, the cluster id has
00975     //appeared before! Such ASSUMPTION is to make implementation easy.
00976     
00977     int curClusterID;
00978     int prevClusterID = eachAtomClusterID[0];
00979     int curClusterSize = 1;
00980     //step1: segment all atoms according to each cluster
00981     for(int i=1; i<g_mol->numAtoms; i++){
00982         curClusterID = eachAtomClusterID[i];
00983         if(curClusterID == prevClusterID){
00984                 curClusterSize++;
00985         }else{
00986                 eachClusterSize.push_back(curClusterSize);
00987                 eachClusterID.push_back(prevClusterID);
00988                 curClusterSize=1;
00989         }
00990         prevClusterID = curClusterID;   
00991     }
00992     //record the info of the last cluster segment
00993     eachClusterSize.push_back(curClusterSize);
00994     eachClusterID.push_back(prevClusterID);
00995     
00996     //step2: detect contiguity of atoms in each cluster and re-scale
00997     //cluster id.
00998     g_isClusterContiguous = 1;    
00999     int *newClusterIDs = new int[eachClusterID.size()];
01000     memset(newClusterIDs, 0, sizeof(int)*eachClusterID.size());
01001     prevClusterID = eachClusterID[0];
01002     int newCId = 0;    
01003     newClusterIDs[0] = newCId;
01004     for(int seg=1; seg<eachClusterID.size(); seg++){
01005         curClusterID = eachClusterID[seg];
01006         if(curClusterID > prevClusterID){
01007                 newClusterIDs[seg] = ++newCId;
01008                 prevClusterID = curClusterID;
01009         }else{
01010                 //non-contiguity happens                
01011                 g_isClusterContiguous = 0;
01012                 
01013                 //we ASSUME this id appears before
01014                 //binary search in eachAtomClusterID[0...seg-1).
01015                 int jl=0, jh=seg-2;
01016                 int isFound = 0;
01017                 while(jh>=jl){
01018                         int mid = (jl+jh)/2;
01019                         if(curClusterID > eachClusterID[mid])
01020                                 jl = mid+1;
01021                         else if(curClusterID < eachClusterID[mid])
01022                                 jh = mid-1;
01023                         else{
01024                                 newClusterIDs[seg] = newClusterIDs[mid];                                
01025                                 isFound = 1;
01026                                 break;
01027                         }                       
01028                 }
01029                 if(!isFound){
01030                         //ASSUMPTION is wrong and abort
01031                         char errmsg[300];
01032                         sprintf(errmsg, "Assumption about building cluster is broken in file %s at line %d\n", __FILE__, __LINE__);
01033                         NAMD_die(errmsg);
01034                 }                               
01035         }       
01036     }
01037     
01038     //step 3: modify eachAtomClusterID according to g_isClusterContiguous
01039     //newCId is the id of the last cluster, as id starts from 0, so the
01040     //total number clusters should be newCId+1
01041     g_numClusters = newCId+1;
01042     if(g_isClusterContiguous){
01043         int aid=0;
01044         for(int seg=0; seg<eachClusterSize.size(); seg++)
01045         {
01046             int curSize = eachClusterSize[seg];
01047             eachAtomClusterID[aid] = curSize;
01048             for(int i=aid+1; i<aid+curSize; i++)
01049                 eachAtomClusterID[i] = -1;
01050             aid += curSize;
01051         }       
01052     }else{
01053         int aid=0;
01054         for(int seg=0; seg<eachClusterSize.size(); seg++)
01055         {
01056             int curSize = eachClusterSize[seg];            
01057             for(int i=aid; i<aid+curSize; i++)
01058                 eachAtomClusterID[i] = newClusterIDs[seg];
01059             aid += curSize;
01060         }
01061     }
01062     free(newClusterIDs);
01063     eachClusterSize.clear();
01064     eachClusterID.clear();
01065 #endif 
01066                 
01067 /*
01068     //check whether cluster is built correctly
01069     printf("num clusters: %d\n", g_numClusters);
01070     FILE *checkFile = fopen("cluster.opt", "w");
01071     for(int i=0; i<g_mol->numAtoms; i++)  fprintf(checkFile, "%d\n", eachAtomClusterID[i]);
01072     fclose(checkFile);
01073 */
01074     
01075     for(int i=0; i<g_mol->numAtoms; i++)
01076         atomListOfBonded[i].clear();
01077     delete [] atomListOfBonded;
01078 #endif
01079 }
01080 
01081 void buildAngleData()
01082 {
01083 #ifndef MEM_OPT_VERSION
01084     Angle *angles = g_mol->getAllAngles();
01085     //create angles' tupleSignature
01086     for(int i=0; i<g_mol->numAngles; i++)
01087     {
01088         Angle *tuple = angles+i;
01089         TupleSignature oneSig(2,ANGLE,tuple->angle_type);
01090         int offset[2];
01091         offset[0] = tuple->atom2 - tuple->atom1;
01092         offset[1] = tuple->atom3 - tuple->atom1;
01093         oneSig.setOffsets(offset);
01094 
01095         int poolIndex = sigsOfAngles.lookupCstPool(oneSig);
01096         if(poolIndex == -1)
01097         {
01098             sigsOfAngles.push_back(oneSig);
01099             poolIndex = (SigIndex)sigsOfAngles.size()-1;
01100         }
01101         eachAtomSigs[tuple->atom1].angleSigIndices.push_back(poolIndex);
01102     }
01103     delete [] angles;
01104 #endif
01105 }
01106 
01107 void buildDihedralData()
01108 {
01109 #ifndef MEM_OPT_VERSION
01110     Dihedral *dihedrals = g_mol->getAllDihedrals();    
01111 
01112     //create dihedrals' tupleSignature
01113     for(int i=0; i<g_mol->numDihedrals; i++)
01114     {
01115         Dihedral *tuple = dihedrals+i;
01116         TupleSignature oneSig(3,DIHEDRAL,tuple->dihedral_type);
01117         int offset[3];
01118         offset[0] = tuple->atom2 - tuple->atom1;
01119         offset[1] = tuple->atom3 - tuple->atom1;
01120         offset[2] = tuple->atom4 - tuple->atom1;
01121         oneSig.setOffsets(offset);
01122 
01123         int poolIndex = sigsOfDihedrals.lookupCstPool(oneSig);
01124         if(poolIndex == -1)
01125         {
01126             sigsOfDihedrals.push_back(oneSig);
01127             poolIndex = (SigIndex)sigsOfDihedrals.size()-1;
01128         }
01129         eachAtomSigs[tuple->atom1].dihedralSigIndices.push_back(poolIndex);
01130     }
01131 
01132     delete[] dihedrals;
01133 #endif
01134 }
01135 
01136 void buildImproperData()
01137 { 
01138 #ifndef MEM_OPT_VERSION
01139     Improper *impropers=g_mol->getAllImpropers();
01140 
01141     //create improper's tupleSignature
01142     for(int i=0; i<g_mol->numImpropers; i++)
01143     {
01144         Improper *tuple = impropers+i;
01145         TupleSignature oneSig(3,IMPROPER,tuple->improper_type);
01146         int offset[3];
01147         offset[0] = tuple->atom2 - tuple->atom1;
01148         offset[1] = tuple->atom3 - tuple->atom1;
01149         offset[2] = tuple->atom4 - tuple->atom1;
01150         oneSig.setOffsets(offset);
01151 
01152         int poolIndex = sigsOfImpropers.lookupCstPool(oneSig);
01153         if(poolIndex == -1)
01154         {
01155             sigsOfImpropers.push_back(oneSig);
01156             poolIndex = (SigIndex)sigsOfImpropers.size()-1;
01157         }
01158         eachAtomSigs[tuple->atom1].improperSigIndices.push_back(poolIndex);
01159     }
01160 
01161     delete[] impropers;
01162 #endif
01163 }
01164 
01165 void buildCrosstermData()
01166 {
01167 #ifndef MEM_OPT_VERSION
01168   Crossterm *crossterms = g_mol->getAllCrossterms();
01169   //create crossterm's tupleSignature
01170   for(int i=0; i<g_mol->numCrossterms; i++)
01171   {
01172     Crossterm *tuple = crossterms+i;
01173     TupleSignature oneSig(7, CROSSTERM, tuple->crossterm_type);
01174     int offset[7];
01175     offset[0] = tuple->atom2 - tuple->atom1;
01176     offset[1] = tuple->atom3 - tuple->atom1;
01177     offset[2] = tuple->atom4 - tuple->atom1;
01178     offset[3] = tuple->atom5 - tuple->atom1;
01179     offset[4] = tuple->atom6 - tuple->atom1;
01180     offset[5] = tuple->atom7 - tuple->atom1;
01181     offset[6] = tuple->atom8 - tuple->atom1;
01182     oneSig.setOffsets(offset);
01183    
01184     int poolIndex = sigsOfCrossterms.lookupCstPool(oneSig);
01185     if(poolIndex == -1)
01186     {
01187       sigsOfCrossterms.push_back(oneSig);
01188       poolIndex = (SigIndex)sigsOfCrossterms.size()-1;
01189     }
01190     eachAtomSigs[tuple->atom1].crosstermSigIndices.push_back(poolIndex);
01191   }
01192   
01193   delete[] crossterms;
01194 #endif
01195 }
01196 
01197 void buildExclusions()
01198 {
01199     //1. Build exclusions: mainly accomplish the function of
01200     //Molecule::build_exclusions (based on the bonds)
01201     UniqueSet<Exclusion> allExclusions;
01202 
01203     int exclude_flag; //Exclusion policy
01204     exclude_flag = g_simParam->exclude;
01205     //int stripHGroupExclFlag = (simParams->splitPatch == SPLIT_PATCH_HYDROGEN);
01206 
01207     //Commented now since no explicit exclusions are read
01208     //  Go through the explicit exclusions and add them to the arrays
01209     //for(i=0; i<numExclusions; i++){
01210     //  exclusionSet.add(exclusions[i]);
01211     //}
01212 
01213     // If this is AMBER force field, and readExclusions is TRUE,
01214     // then all the exclusions were read from parm file, and we
01215     // shouldn't generate any of them.
01216     // Comment on stripHGroupExcl:
01217     // 1. Inside this function, hydrogenGroup is initialized in
01218     // build_atom_status, therefore, not available when reading psf files
01219     // 2. this function's main purpose is to reduce memory usage. Since exclusion
01220     // signatures are used, this function could be overlooked  --Chao Mei
01221 
01222     vector<int> *eachAtomNeighbors = new vector<int>[g_mol->numAtoms];   
01223     for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
01224     {
01225         AtomSigInfo *aSig = &atomSigPool[atomData[atom1].atomSigIdx];
01226         for(int j=0; j<aSig->bondSigIndices.size(); j++)
01227         {
01228             TupleSignature *tSig = &sigsOfBonds[aSig->bondSigIndices[j]];
01229             if(!tSig->isReal) continue;
01230             int atom2 = atom1+tSig->offset[0];
01231             eachAtomNeighbors[atom1].push_back(atom2);
01232             eachAtomNeighbors[atom2].push_back(atom1);
01233         }
01234     }
01235 
01236     if (!g_simParam->amberOn || !g_simParam->readExclusions)
01237     { //  Now calculate the bonded exlcusions based on the exclusion policy
01238         switch (exclude_flag)
01239         {
01240         case NONE:
01241             break;
01242         case ONETWO:
01243             build12Excls(allExclusions, eachAtomNeighbors);
01244             break;
01245         case ONETHREE:
01246             build12Excls(allExclusions, eachAtomNeighbors);
01247             build13Excls(allExclusions, eachAtomNeighbors);
01248             //if ( stripHGroupExclFlag ) stripHGroupExcl();
01249             break;
01250         case ONEFOUR:
01251             build12Excls(allExclusions, eachAtomNeighbors);
01252             build13Excls(allExclusions, eachAtomNeighbors);
01253             build14Excls(allExclusions, eachAtomNeighbors, 0);
01254             //if ( stripHGroupExclFlag ) stripHGroupExcl();
01255             break;
01256         case SCALED14:
01257             build12Excls(allExclusions, eachAtomNeighbors);
01258             build13Excls(allExclusions, eachAtomNeighbors);
01259             build14Excls(allExclusions, eachAtomNeighbors, 1);
01260             //if ( stripHGroupExclFlag ) stripHGroupExcl();
01261             break;
01262         }
01263     }
01264     //else if (stripHGroupExclFlag && exclude_flag!=NONE && exclude_flag!=ONETWO)
01265     //  stripHGroupExcl();
01266 
01267     //Commented since atomFepFlags information is not available when reading psf file
01268     //stripFepExcl();
01269 
01270     for(int i=0; i<g_mol->numAtoms; i++)
01271         eachAtomNeighbors[i].clear();
01272     delete [] eachAtomNeighbors;
01273 
01274     //2. Build each atom's list of exclusions
01275     iout << iINFO << "ADDED " << allExclusions.size() << " IMPLICIT EXCLUSIONS\n" << endi;
01276     UniqueSetIter<Exclusion> exclIter(allExclusions);
01277     eachAtomExclSigs = new ExclSigInfo[g_mol->numAtoms];
01278     for(exclIter=exclIter.begin(); exclIter!=exclIter.end(); exclIter++)
01279     {
01280         int atom1 = exclIter->atom1;
01281         int atom2 = exclIter->atom2;
01282         int offset21 = atom2-atom1;
01283         if(exclIter->modified)
01284         {
01285             eachAtomExclSigs[atom1].modExclOffset.push_back(offset21);
01286             eachAtomExclSigs[atom2].modExclOffset.push_back(-offset21);
01287         }
01288         else
01289         {
01290             eachAtomExclSigs[atom1].fullExclOffset.push_back(offset21);
01291             eachAtomExclSigs[atom2].fullExclOffset.push_back(-offset21);
01292         }
01293     }
01294     allExclusions.clear();
01295 
01296     //3. Build up exclusion signatures and determine each atom's
01297     //exclusion signature index
01298     for(int i=0; i<g_mol->numAtoms; i++)
01299     {
01300         eachAtomExclSigs[i].sortExclOffset();
01301         int poolIndex = sigsOfExclusions.lookupCstPool(eachAtomExclSigs[i]);
01302         if(poolIndex==-1)
01303         {
01304             poolIndex = sigsOfExclusions.size();
01305             sigsOfExclusions.push_back(eachAtomExclSigs[i]);
01306         }
01307         atomData[i].exclSigIdx = poolIndex;
01308     }
01309     delete [] eachAtomExclSigs;
01310     eachAtomExclSigs = NULL;
01311     printf("Exclusion signatures: %d\n", (int)sigsOfExclusions.size());
01312 }
01313 
01314 void build12Excls(UniqueSet<Exclusion>& allExcls, vector<int> *eachAtomNeighbors)
01315 {
01316     for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
01317     {
01318         vector<int> *atom1List = &eachAtomNeighbors[atom1];
01319         for(int j=0; j<atom1List->size(); j++)
01320         {
01321             int atom2 = atom1List->at(j);
01322             if(atom1<atom2)
01323                 allExcls.add(Exclusion(atom1, atom2));
01324             else
01325                 allExcls.add(Exclusion(atom2, atom1));
01326         }
01327     }
01328 }
01329 
01330 void build13Excls(UniqueSet<Exclusion>& allExcls, vector<int> *eachAtomNeighbors)
01331 {
01332     for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
01333     {
01334         vector<int> *atom1List = &eachAtomNeighbors[atom1];
01335         for(int j=0; j<atom1List->size(); j++)
01336         {
01337             int atom2 = atom1List->at(j);
01338             vector<int> *atom2List = &eachAtomNeighbors[atom2];
01339             for(int k=0; k<atom2List->size(); k++)
01340             {
01341                 int atom3 = atom2List->at(k);
01342                 //atom1-atom2, so atom2List contains atom1 which should not be considered
01343                 if(atom3 == atom1)
01344                     continue;
01345                 if(atom1<atom3)
01346                     allExcls.add(Exclusion(atom1, atom3));
01347                 else
01348                     allExcls.add(Exclusion(atom3, atom1));
01349             }
01350         }
01351     }
01352 }
01353 
01354 void build14Excls(UniqueSet<Exclusion>& allExcls, vector<int> *eachAtomNeighbors, int modified)
01355 {
01356     for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
01357     {
01358         vector<int> *atom1List = &eachAtomNeighbors[atom1];
01359         for(int j=0; j<atom1List->size(); j++)
01360         {
01361             int atom2 = atom1List->at(j);
01362             vector<int> *atom2List = &eachAtomNeighbors[atom2];
01363             for(int k=0; k<atom2List->size(); k++)
01364             {
01365                 int atom3 = atom2List->at(k);
01366                 //atom1-atom2, so atom2List contains atom1 which should not be considered
01367                 if(atom3 == atom1)
01368                     continue;
01369                 vector<int> *atom3List = &eachAtomNeighbors[atom3];
01370                 for(int l=0; l<atom3List->size(); l++)
01371                 {
01372                     int atom4 = atom3List->at(l);
01373                     //atom1-atom2, so atom2List contains atom1 which should not be considered
01374                     if(atom4 == atom2 || atom4 == atom1)
01375                         continue;
01376                     if(atom1<atom4)
01377                         allExcls.add(Exclusion(atom1, atom4, modified));
01378                     else
01379                         allExcls.add(Exclusion(atom4, atom1, modified));
01380                 }
01381             }
01382         }
01383     }
01384 }
01385 
01386 template <class T> void HashPool<T>::dump_tables()
01387 {
01388   for(int j=0; j < pool.size(); j++) {
01389     HashPoolAdaptorT<T>* pval = pool[j];
01390     CmiPrintf("Pool[%d]=%p %p  hash = %d\n",j,pool[j],pval,pval->hash());
01391   }
01392   CkHashtableIterator *iter = index_table.iterator();
01393   void *key,*indx;
01394   while (iter->hasNext()) {
01395     indx = iter->next(&key);
01396     HashPoolAdaptorT<T> *pkey = (HashPoolAdaptorT<T>*)key;
01397     CmiPrintf("key %p indx %p %d hash=%d\n",key,indx,*((int *)indx),pkey->hash());
01398   }
01399 }

Generated on Mon Nov 20 01:17:10 2017 for NAMD by  doxygen 1.4.7