NAMD
CompressPsf.C
Go to the documentation of this file.
1 #include <algorithm>
2 #include "CompressPsf.h"
3 #include "strlib.h"
4 #include "Molecule.h"
5 #include "Parameters.h"
6 #include "SimParameters.h"
7 #include "InfoStream.h"
8 #include "UniqueSet.h"
9 #include "UniqueSetIter.h"
10 
11 using namespace std;
12 
31 //global variables recording compressed psf file information
32 Molecule *g_mol = NULL;
36 
38 {
47  int resID;
48 };
49 
51  flipNum((char *)&sSet, sizeof(short), sizeof(sSet)/sizeof(short));
52  flipNum((char *)&iSet, sizeof(int), sizeof(iSet)/sizeof(int));
53  flipNum((char *)&fSet, sizeof(float), sizeof(fSet)/sizeof(float));
54 }
55 
57 {
58  vector<SigIndex> bondSigIndices;
59  vector<SigIndex> angleSigIndices;
60  vector<SigIndex> dihedralSigIndices;
61  vector<SigIndex> improperSigIndices;
62  vector<SigIndex> crosstermSigIndices;
63 
65  {}
67  {
68  bondSigIndices.clear();
69  for(int i=0; i<sig.bondSigIndices.size(); i++)
70  bondSigIndices.push_back(sig.bondSigIndices[i]);
71 
72  angleSigIndices.clear();
73  for(int i=0; i<sig.angleSigIndices.size(); i++)
74  angleSigIndices.push_back(sig.angleSigIndices[i]);
75 
76  dihedralSigIndices.clear();
77  for(int i=0; i<sig.dihedralSigIndices.size(); i++)
78  dihedralSigIndices.push_back(sig.dihedralSigIndices[i]);
79 
80  improperSigIndices.clear();
81  for(int i=0; i<sig.improperSigIndices.size(); i++)
82  improperSigIndices.push_back(sig.improperSigIndices[i]);
83 
84  crosstermSigIndices.clear();
85  for(int i=0; i<sig.crosstermSigIndices.size(); i++)
86  crosstermSigIndices.push_back(sig.crosstermSigIndices[i]);
87  }
88 
90  {
91  bondSigIndices.clear();
92  angleSigIndices.clear();
93  dihedralSigIndices.clear();
94  improperSigIndices.clear();
95  crosstermSigIndices.clear();
96  }
97 
99  {
100  sort(bondSigIndices.begin(), bondSigIndices.end());
101  sort(angleSigIndices.begin(), angleSigIndices.end());
102  sort(dihedralSigIndices.begin(), dihedralSigIndices.end());
103  sort(improperSigIndices.begin(), improperSigIndices.end());
104  sort(crosstermSigIndices.begin(), crosstermSigIndices.end());
105  }
106 
107  inline CkHashCode hash() const {
108  // What's a good hash function for this? Lets make something up
109  // Concatenate all the index lists into a list of chars, then hash that
110  // string using Charm's string hash function
111 
112  // To keep memory allocation cheap, we'll just use a 32-byte buffer
113  // and wrap around if we have more sigs
114  const int maxlen = 32;
115  unsigned char keydata[maxlen+1];
116  const int maxchar = 256;
117  int i,j;
118  for(j=0;j<=maxlen;j++) keydata[j] = 0;
119  j=0;
120  for(i=0; i<bondSigIndices.size(); i++,j++) {
121  keydata[j % maxlen] ^= (bondSigIndices[i] % maxchar);
122  }
123  for(i=0; i<angleSigIndices.size(); i++,j++) {
124  keydata[j % maxlen] ^= (angleSigIndices[i] % maxchar);
125  }
126  for(i=0; i<dihedralSigIndices.size(); i++,j++) {
127  keydata[j % maxlen] ^= (dihedralSigIndices[i] % maxchar);
128  }
129  for(i=0; i<improperSigIndices.size(); i++,j++) {
130  keydata[j % maxlen] ^= (improperSigIndices[i] % maxchar);
131  }
132  for(i=0; i<crosstermSigIndices.size(); i++,j++) {
133  keydata[j % maxlen] ^= (crosstermSigIndices[i] % maxchar);
134  }
135 // CmiPrintf("Computed hash string len %d,%d\n",j,maxlen);
136  if (j > maxlen) j = maxlen;
137 // for(i=0; i < j; i++) {
138 // if (keydata[i] == 0)
139 // keydata[i] = 255;
140 // CmiPrintf("key[%d]=%d %p\n",i,keydata[i],keydata);
141 // }
142  return CkHashFunction_default((const void*)keydata,(size_t)j);
143  }
144 };
145 
146 int operator==(const AtomSigInfo &s1, const AtomSigInfo& s2)
147 {
148  if(s1.bondSigIndices.size() != s2.bondSigIndices.size())
149  return 0;
150  if(s1.angleSigIndices.size() != s2.angleSigIndices.size())
151  return 0;
152  if(s1.dihedralSigIndices.size() != s2.dihedralSigIndices.size())
153  return 0;
154  if(s1.improperSigIndices.size() != s2.improperSigIndices.size())
155  return 0;
156  if(s1.crosstermSigIndices.size() != s2.crosstermSigIndices.size())
157  return 0;
158 
159  int equalCnt;
160  equalCnt=0;
161  int bondSigCnt = s1.bondSigIndices.size();
162  for(int i=0; i<bondSigCnt; i++)
163  equalCnt += (s1.bondSigIndices[i]==s2.bondSigIndices[i]);
164  if(equalCnt!=bondSigCnt)
165  return 0;
166 
167  equalCnt=0;
168  int angleSigCnt = s1.angleSigIndices.size();
169  for(int i=0; i<angleSigCnt; i++)
170  equalCnt += (s1.angleSigIndices[i]==s2.angleSigIndices[i]);
171  if(equalCnt!=angleSigCnt)
172  return 0;
173 
174  equalCnt=0;
175  int dihedralSigCnt = s1.dihedralSigIndices.size();
176  for(int i=0; i<dihedralSigCnt; i++)
177  equalCnt += (s1.dihedralSigIndices[i]==s2.dihedralSigIndices[i]);
178  if(equalCnt!=dihedralSigCnt)
179  return 0;
180 
181  equalCnt=0;
182  int improperSigCnt = s1.improperSigIndices.size();
183  for(int i=0; i<improperSigCnt; i++)
184  equalCnt += (s1.improperSigIndices[i]==s2.improperSigIndices[i]);
185  if(equalCnt!=improperSigCnt)
186  return 0;
187 
188  equalCnt=0;
189  int crosstermSigCnt = s1.crosstermSigIndices.size();
190  for(int i=0; i<crosstermSigCnt; i++)
191  equalCnt += (s1.crosstermSigIndices[i]==s2.crosstermSigIndices[i]);
192  if(equalCnt!=crosstermSigCnt)
193  return 0;
194 
195  return 1;
196 }
197 
199 {
200  vector<int> fullExclOffset;
201  vector<int> modExclOffset;
202 
204  {}
206  {
207  fullExclOffset.clear();
208  for(int i=0; i<sig.fullExclOffset.size(); i++)
209  fullExclOffset.push_back(sig.fullExclOffset[i]);
210 
211  modExclOffset.clear();
212  for(int i=0; i<sig.modExclOffset.size(); i++)
213  modExclOffset.push_back(sig.modExclOffset[i]);
214  }
215 
217  {
218  fullExclOffset.clear();
219  modExclOffset.clear();
220  }
221 
223  {
224  sort(fullExclOffset.begin(), fullExclOffset.end());
225  sort(modExclOffset.begin(), modExclOffset.end());
226  }
227 
228  int hash() const {
229  unsigned int code = 0x1234;
230  unsigned int codesz = 8 * sizeof(int);
231  const unsigned int numFoffset = fullExclOffset.size();
232  const unsigned int numMoffset = modExclOffset.size();
233  const unsigned int numOffsets = numFoffset + numMoffset;
234 
235  // No excluded atoms? Just hash to 0
236  if (numOffsets == 0)
237  return 0;
238 
239  unsigned int shift = codesz / numOffsets;
240  if (shift == 0) shift=1;
241  unsigned int i;
242  for(i=0; i < numFoffset; i++) {
243  code = circShift(code,shift);
244  code ^= fullExclOffset[i];
245  }
246  for(i=0; i < numMoffset; i++) {
247  code = circShift(code,shift);
248  code ^= modExclOffset[i];
249  }
250  return code;
251  }
252 };
253 int operator==(const ExclSigInfo &s1, const ExclSigInfo &s2)
254 {
255  if(s1.fullExclOffset.size()!=s2.fullExclOffset.size())
256  return 0;
257  if(s1.modExclOffset.size()!=s2.modExclOffset.size())
258  return 0;
259 
260  for(int i=0; i<s1.fullExclOffset.size(); i++)
261  {
262  if(s1.fullExclOffset[i] != s2.fullExclOffset[i])
263  return 0;
264  }
265 
266  for(int i=0; i<s1.modExclOffset.size(); i++)
267  {
268  if(s1.modExclOffset[i] != s2.modExclOffset[i])
269  return 0;
270  }
271  return 1;
272 }
273 
274 class HashString : public string {
275 public:
276  int hash() const {
277  const char* d = this->c_str();
278  int ret=0;
279  for (int i=0;d[i]!=0;i++) {
280  int shift1=((5*i)%16)+0;
281  int shift2=((6*i)%16)+8;
282  ret+=((0xa5^d[i])<<shift2)+(d[i]<<shift1);
283  }
284  return ret;
285  }
286 };
287 
288 class HashReal {
289  Real val;
290 public:
291  HashReal(Real v) : val(v) {}
292  operator Real & () { return val; }
293  operator const Real & () const { return val; }
294 
295  int hash() const {
296  const char* d = (const char *)&val;
297  int ret=0;
298  for (int i=0;i < sizeof(Real);i++) {
299  int shift1=((5*i)%16)+0;
300  int shift2=((6*i)%16)+8;
301  ret+=((0xa5^d[i])<<shift2)+(d[i]<<shift1);
302  }
303  return ret;
304  };
305 };
306 
315 
316 //Recording cluster information after reading all bonds info
317 int *eachAtomClusterID = NULL;
318 vector<int> eachClusterSize;
319 vector<int> eachClusterID;
321 
328 
331 
332 //Structures for extraBond options
333 vector<Bond> extraBonds;
334 vector<Angle> extraAngles;
335 vector<Dihedral> extraDihedrals;
336 vector<Improper> extraImpropers;
337 
338 vector<BondValue> extraBondParams;
339 vector<AngleValue> extraAngleParams;
340 vector<DihedralValue> extraDihedralParams;
341 vector<ImproperValue> extraImproperParams;
342 
343 int operator==(const BondValue &b1, const BondValue &b2)
344 {
345  return (b1.k==b2.k) && (b1.x0==b2.x0);
346 }
347 
348 int operator==(const AngleValue &a1, const AngleValue &a2)
349 {
350  return (a1.k==a2.k) && (a1.k_ub==a2.k_ub) && (a1.r_ub==a2.r_ub) && (a1.theta0==a2.theta0);
351 }
352 
353 int operator!=(const FourBodyConsts& f1, const FourBodyConsts& f2)
354 {
355  return (f1.delta!=f2.delta) || (f1.k!=f2.k) || (f1.n!=f2.n);
356 }
357 
358 int operator==(const DihedralValue &d1, const DihedralValue &d2)
359 {
360  if(d1.multiplicity != d2.multiplicity)
361  return 0;
362  for(int i=0; i<MAX_MULTIPLICITY; i++)
363  {
364  if(d1.values[i] != d2.values[i])
365  return 0;
366  }
367  return 1;
368 }
369 
370 int operator==(const ImproperValue &d1, const ImproperValue &d2)
371 {
372  if(d1.multiplicity != d2.multiplicity)
373  return 0;
374  for(int i=0; i<MAX_MULTIPLICITY; i++)
375  {
376  if(d1.values[i] != d2.values[i])
377  return 0;
378  }
379  return 1;
380 }
381 
382 void loadMolInfo();
383 
384 void integrateAllAtomSigs();
385 void outputCompressedFile(FILE *txtOfp, FILE *binOfp);
386 
387 //reading extraBond's information
388 void getExtraBonds(StringList *file);
389 
390 void buildAtomData();
391 void buildBondData();
392 void buildAngleData();
393 void buildDihedralData();
394 void buildImproperData();
395 void buildCrosstermData();
396 
397 void buildExclusionData();
398 
399 //Functions related with building exclusions
400 void buildExclusions();
401 void build12Excls(UniqueSet<Exclusion>&, vector<int> *);
402 void build13Excls(UniqueSet<Exclusion>&, vector<int> *);
403 void build14Excls(UniqueSet<Exclusion>&, vector<int> *, int);
404 
405 //reverse the byte-order of every element starting at "elem"
406 void flipNum(char *elem, int elemSize, int numElems){
407  int mid = elemSize/2;
408  char *ptr = elem;
409  for(int i=0; i<numElems; i++) {
410  for(int j=0; j<mid; j++) {
411  char tmp = ptr[j];
412  ptr[j] = ptr[elemSize-1-j];
413  ptr[elemSize-1-j] = tmp;
414  }
415  ptr += elemSize;
416  }
417 }
418 
420 {
421  segNamePool.clear();
422  resNamePool.clear();
423  atomNamePool.clear();
424  atomTypePool.clear();
425  chargePool.clear();
426  massPool.clear();
427  sigsOfBonds.clear();
428  sigsOfAngles.clear();
429  sigsOfDihedrals.clear();
430  sigsOfImpropers.clear();
431  sigsOfExclusions.clear();
432 
433  eachClusterSize.clear();
434 }
435 
436 void compress_molecule_info(Molecule *mol, char *psfFileName, Parameters *param, SimParameters *simParam, ConfigList* cfgList)
437 {
438  g_mol = mol;
439  g_param = param;
440  g_simParam = simParam; //used for building exclusions
441  g_cfgList = cfgList; //used for integrating extra bonds
442 
443  //read psf files
444  //readPsfFile(psfFileName);
445  loadMolInfo();
446 
448 
449  buildExclusions();
450 
451  //buildParamData();
452 
453  char *outFileName = new char[strlen(psfFileName)+20];
454  sprintf(outFileName, "%s.inter", psfFileName);
455  //the text file for signatures and other non-per-atom info
456  FILE *txtOfp = fopen(outFileName, "w");
457  sprintf(outFileName, "%s.inter.bin", psfFileName);
458  //the binary file for per-atom info
459  FILE *binOfp = fopen(outFileName, "wb");
460  delete [] outFileName;
461 
462  //output compressed psf file
463  outputCompressedFile(txtOfp, binOfp);
464 
465  fclose(txtOfp);
466  fclose(binOfp);
467 }
468 
474 {
475  char err_msg[512]; // Error message for NAMD_die
476  char buffer[512]; // Buffer for file reading
477  int i; // Loop counter
478  int NumTitle; // Number of Title lines in .psf file
479  int ret_code; // ret_code from NAMD_read_line calls
480  FILE *psf_file;
481  Parameters *params = g_param;
482 
483  buildAtomData();
484 
485  //read extra bonds/angles/dihedrals/impropers information first
486  //and then integrate them into the following reading procedures
487  /*if(g_simParam->extraBondsOn)
488  {
489  getExtraBonds(g_cfgList->find("extraBondsFile"));
490  }*/
491 
492  //initialize eachAtomSigs
494 
495  buildBondData();
496  buildAngleData();
500 
501  // analyze the data and find the status of each atom
502  //build_atom_status();
503 }
504 
506 {
507  printf("Bond sigs: %d\n", (int)sigsOfBonds.size());
508  printf("Angle sigs: %d\n", (int)sigsOfAngles.size());
509  printf("Dihedral sigs: %d\n", (int)sigsOfDihedrals.size());
510  printf("Improper sigs: %d\n", (int)sigsOfImpropers.size());
511  printf("Crossterm sigs: %d\n", (int)sigsOfCrossterms.size());
512 
513 
514  for(int i=0; i<g_mol->numAtoms; i++)
515  {
517  int poolIndex = atomSigPool.lookupCstPool(eachAtomSigs[i]);
518  if(poolIndex==-1)
519  {
520  atomSigPool.push_back(eachAtomSigs[i]);
521  poolIndex = atomSigPool.size()-1;
522  }
523  atomData[i].atomSigIdx = poolIndex;
524  }
525 
526  printf("Atom's sigs: %d\n", (int)atomSigPool.size());
527 
528  delete[] eachAtomSigs;
529 }
530 
541 void outputCompressedFile(FILE *txtOfp, FILE *binOfp)
542 {
543 #ifndef MEM_OPT_VERSION
544  fprintf(txtOfp, "FORMAT VERSION: %f\n", COMPRESSED_PSF_VER);
545 
546  fprintf(txtOfp, "%d !NSEGMENTNAMES\n", segNamePool.size());
547  for(int i=0; i<segNamePool.size(); i++)
548  {
549  fprintf(txtOfp, "%s\n", segNamePool[i].c_str());
550  }
551 
552  fprintf(txtOfp, "%d !NRESIDUENAMES\n", resNamePool.size());
553  for(int i=0; i<resNamePool.size(); i++)
554  {
555  fprintf(txtOfp, "%s\n", resNamePool[i].c_str());
556  }
557 
558  fprintf(txtOfp, "%d !NATOMNAMES\n", atomNamePool.size());
559  for(int i=0; i<atomNamePool.size(); i++)
560  {
561  fprintf(txtOfp, "%s\n", atomNamePool[i].c_str());
562  }
563 
564  fprintf(txtOfp, "%d !NATOMTYPES\n", atomTypePool.size());
565  for(int i=0; i<atomTypePool.size(); i++)
566  {
567  fprintf(txtOfp, "%s\n", atomTypePool[i].c_str());
568  }
569 
570  fprintf(txtOfp, "%d !NCHARGES\n", chargePool.size());
571  for(int i=0; i<chargePool.size(); i++)
572  {
573  const Real charge = chargePool[i];
574  fprintf(txtOfp, "%f\n", charge);
575  }
576 
577  fprintf(txtOfp, "%d !NMASSES\n", massPool.size());
578  for(int i=0; i<massPool.size(); i++)
579  {
580  const Real mass = massPool[i];
581  fprintf(txtOfp, "%f\n", mass);
582  }
583 
584 
585  fprintf(txtOfp, "%d !NATOMSIGS\n", atomSigPool.size());
586  for(int i=0; i<atomSigPool.size(); i++)
587  {
588  AtomSigInfo& oneAtomSig = atomSigPool[i];
589  int oneTypeCnt = oneAtomSig.bondSigIndices.size();
590  fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NBOND");
591  for(int j=0; j<oneTypeCnt; j++)
592  {
593  SigIndex idx = oneAtomSig.bondSigIndices[j];
594  TupleSignature& tSig = sigsOfBonds[idx];
595  tSig.output(txtOfp);
596  }
597 
598  oneTypeCnt = oneAtomSig.angleSigIndices.size();
599  fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NTHETA");
600  for(int j=0; j<oneTypeCnt; j++)
601  {
602  SigIndex idx = oneAtomSig.angleSigIndices[j];
603  TupleSignature& tSig = sigsOfAngles[idx];
604  tSig.output(txtOfp);
605  }
606 
607  oneTypeCnt = oneAtomSig.dihedralSigIndices.size();
608  fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NPHI");
609  for(int j=0; j<oneTypeCnt; j++)
610  {
611  SigIndex idx = oneAtomSig.dihedralSigIndices[j];
612  TupleSignature& tSig = sigsOfDihedrals[idx];
613  tSig.output(txtOfp);
614  }
615 
616  oneTypeCnt = oneAtomSig.improperSigIndices.size();
617  fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NIMPHI");
618  for(int j=0; j<oneTypeCnt; j++)
619  {
620  SigIndex idx = oneAtomSig.improperSigIndices[j];
621  TupleSignature& tSig = sigsOfImpropers[idx];
622  tSig.output(txtOfp);
623  }
624 
625  oneTypeCnt = oneAtomSig.crosstermSigIndices.size();
626  fprintf(txtOfp, "%d !%sSIGS\n", oneTypeCnt, "NCRTERM");
627  for(int j=0; j<oneTypeCnt; j++)
628  {
629  SigIndex idx = oneAtomSig.crosstermSigIndices[j];
630  TupleSignature& tSig = sigsOfCrossterms[idx];
631  tSig.output(txtOfp);
632  }
633  }
634 
635  //2. Output exclusion signatures
636  int exclSigCnt = sigsOfExclusions.size();
637  fprintf(txtOfp, "%d !NEXCLSIGS\n", exclSigCnt);
638  for(int i=0; i<exclSigCnt; i++)
639  {
640  ExclSigInfo *sig = &sigsOfExclusions[i];
641  //first line is for full exclusions (1-2, 1-3) in the format of count offset1 offset2 offset3 ...
642  fprintf(txtOfp, "%lu", (unsigned long) sig->fullExclOffset.size());
643  for(int j=0; j<sig->fullExclOffset.size(); j++)
644  fprintf(txtOfp, " %d", sig->fullExclOffset[j]);
645  fprintf(txtOfp, "\n");
646 
647  //second line is for modified exclusions (1-4)
648  fprintf(txtOfp, "%lu", (unsigned long) sig->modExclOffset.size());
649  for(int j=0; j<sig->modExclOffset.size(); j++)
650  fprintf(txtOfp, " %d", sig->modExclOffset[j]);
651  fprintf(txtOfp, "\n");
652  }
653 
654  //3. Output the cluster information
655  fprintf(txtOfp, "%d !NCLUSTERS\n", g_numClusters);
656 
657  //4. Output atom info
658  fprintf(txtOfp, "%d !NATOM\n", g_mol->numAtoms);
659  fprintf(txtOfp, "%d !NHYDROGENGROUP\n", g_mol->numHydrogenGroups);
660  fprintf(txtOfp, "%d !MAXHYDROGENGROUPSIZE\n", g_mol->maxHydrogenGroupSize);
661  fprintf(txtOfp, "%d !NMIGRATIONGROUP\n", g_mol->numMigrationGroups);
662  fprintf(txtOfp, "%d !MAXMIGRATIONGROUPSIZE\n", g_mol->maxMigrationGroupSize);
663 
664  //5. Output rigid bond type
665  fprintf(txtOfp, "%d !RIGIDBONDTYPE\n", g_simParam->rigidBonds);
666 #if 0
667  const float *atomOccupancy = g_mol->getOccupancyData();
668  const float *atomBFactor = g_mol->getBFactorData();
669  fprintf(txtOfp, "%d !OCCUPANCYVALID\n", (atomOccupancy==NULL)?0:1);
670  fprintf(txtOfp, "%d !TEMPFACTORVALID\n", (atomBFactor==NULL)?0:1);
671 
672  float *zeroFloats = NULL;
673  if(atomOccupancy==NULL || atomBFactor==NULL) {
674  zeroFloats = new float[g_mol->numAtoms];
675  memset(zeroFloats, 0, sizeof(float)*g_mol->numAtoms);
676  if(atomOccupancy==NULL) atomOccupancy = (const float *)zeroFloats;
677  if(atomBFactor==NULL) atomBFactor = (const float *)zeroFloats;
678  }
679 #endif
680 
681  Atom *atoms = g_mol->getAtoms(); //need to output its partner and hydrogenList
683 
684  //First, output magic number
685  int magicNum = COMPRESSED_PSF_MAGICNUM;
686  fwrite(&magicNum, sizeof(int), 1, binOfp);
687  //Second, version number
688  float verNum = (float)COMPRESSED_PSF_VER;
689  fwrite(&verNum, sizeof(float), 1, binOfp);
690  //Third, the per-atom record size
691  int recSize = sizeof(OutputAtomRecord);
692  fwrite(&recSize, sizeof(int), 1, binOfp);
693  //Fourth, each atom info
694  OutputAtomRecord oneRec;
695  for(int i=0; i<g_mol->numAtoms; i++)
696  {
697  oneRec.sSet.segNameIdx = atomData[i].segNameIdx;
698  oneRec.sSet.resNameIdx = atomData[i].resNameIdx;
699  oneRec.sSet.atomNameIdx = atomData[i].atomNameIdx;
700  oneRec.sSet.atomTypeIdx = atomData[i].atomTypeIdx;
701  oneRec.sSet.chargeIdx = atomData[i].chargeIdx;
702  oneRec.sSet.massIdx = atomData[i].massIdx;
703  oneRec.iSet.atomSigIdx = atomData[i].atomSigIdx;
704 
705  oneRec.iSet.exclSigIdx = atomData[i].exclSigIdx;
706  oneRec.sSet.vdw_type = atoms[i].vdw_type;
707  oneRec.iSet.resID = atomData[i].resID;
708  int hydIdx = atoms[i].hydrogenList;
709  oneRec.iSet.hydrogenList = hydIdx;
710  oneRec.iSet.atomsInGroup = hg[hydIdx].atomsInGroup;
711  oneRec.iSet.GPID = hg[hydIdx].GPID;
712  //oneRec.waterVal = hg[hydIdx].waterVal;
714  oneRec.iSet.MPID = hg[hydIdx].MPID;
715  //oneRec.fSet.occupancy = atomOccupancy[i];
716  //oneRec.fSet.bfactor = atomBFactor[i];
718  fwrite(&oneRec, sizeof(OutputAtomRecord), 1, binOfp);
719  }
720  //if(zeroFloats) delete[] zeroFloats;
721  delete[] atomData;
722 
723  //Output the info required for parallel output:
724  //1. Cluster ids of each atom
725  //Since the cluster info is only when doing output under
726  //wrapAll or wrapWater conditions, for the purpose of parallel
727  //output, it is reasonable to separate the cluster info from
728  //the above per-atom info. So whole the binary per-atom file
729  //contains two parts, one for parallel input to create patches
730  //and computes; the other for parallel output -Chao Mei
731  //Fifth: output the cluster id of each atoms
732  fwrite(eachAtomClusterID, sizeof(int), g_mol->numAtoms, binOfp);
733  //2. Whether each atom is water or not
734  char *isWater = new char[g_mol->numAtoms];
735  for(int i=0; i<g_mol->numAtoms; i++){
736  isWater[i] = g_mol->is_water(i);
737  }
738  fwrite(isWater, sizeof(char), g_mol->numAtoms, binOfp);
739  delete [] isWater;
740  delete[] atoms;
742  delete[] eachAtomClusterID;
743 
744  //Output the parameter new values if extraBonds are present.
745  //The parameters are not needed since now extraBonds' parameters will be
746  //read again during running the simulation
747 
748  //6. Output the "multiplicity" field TUPLE_array of the Parameter object
749  fprintf(txtOfp, "!DIHEDRALPARAMARRAY\n");
750  for(int i=0; i<g_param->NumDihedralParams; i++)
751  {
752  fprintf(txtOfp, "%d ", g_param->dihedral_array[i].multiplicity);
753  }
754  fprintf(txtOfp, "\n");
755  fprintf(txtOfp, "!IMPROPERPARAMARRAY\n");
756  for(int i=0; i<g_param->NumImproperParams; i++)
757  {
758  fprintf(txtOfp, "%d ", g_param->improper_array[i].multiplicity);
759  }
760  fprintf(txtOfp, "\n");
761 #endif
762 }
763 
765 {
766 #ifndef MEM_OPT_VERSION
767  int numAtoms = g_mol->numAtoms;
768 
769  //1. parse atom data to build constant pool (atom name, mass, charge etc.)
770  atomData = new BasicAtomInfo[numAtoms];
771  Atom *atoms = g_mol->getAtoms();
772  AtomNameInfo *atomNames = g_mol->getAtomNames();
773  AtomSegResInfo *atomSegResids = g_mol->getAtomSegResInfo();
774 
775  for(int atomID=0; atomID < numAtoms; atomID++)
776  {
777  //building constant pool
778  int poolIndex;
779  HashString fieldName;
780  fieldName.assign(atomSegResids[atomID].segname);
781  poolIndex = segNamePool.lookupCstPool(fieldName);
782  if(poolIndex==-1)
783  {
784  segNamePool.push_back(fieldName);
785  poolIndex = segNamePool.size()-1;
786  }
787  atomData[atomID].segNameIdx = poolIndex;
788 
789  atomData[atomID].resID = atomSegResids[atomID].resid;
790 
791  fieldName.assign(atomNames[atomID].resname);
792  poolIndex = resNamePool.lookupCstPool(fieldName);
793  if(poolIndex==-1)
794  {
795  resNamePool.push_back(fieldName);
796  poolIndex = resNamePool.size()-1;
797  }
798  atomData[atomID].resNameIdx = poolIndex;
799 
800  fieldName.assign(atomNames[atomID].atomname);
801  poolIndex = atomNamePool.lookupCstPool(fieldName);
802  if(poolIndex==-1)
803  {
804  atomNamePool.push_back(fieldName);
805  poolIndex = atomNamePool.size()-1;
806  }
807  atomData[atomID].atomNameIdx = poolIndex;
808 
809  fieldName.assign(atomNames[atomID].atomtype);
810  poolIndex = atomTypePool.lookupCstPool(fieldName);
811  if(poolIndex==-1)
812  {
813  atomTypePool.push_back(fieldName);
814  poolIndex = atomTypePool.size()-1;
815  }
816  atomData[atomID].atomTypeIdx = poolIndex;
817 
818  poolIndex = chargePool.lookupCstPool(atoms[atomID].charge);
819  if(poolIndex==-1)
820  {
821  chargePool.push_back(atoms[atomID].charge);
822  poolIndex = chargePool.size()-1;
823  }
824  atomData[atomID].chargeIdx = poolIndex;
825 
826  poolIndex = massPool.lookupCstPool(atoms[atomID].mass);
827  if(poolIndex==-1)
828  {
829  massPool.push_back(atoms[atomID].mass);
830  poolIndex = massPool.size()-1;
831  }
832  atomData[atomID].massIdx = poolIndex;
833  }
834 
835  //Free those space to reduce transient memory usage
836  //delete [] atoms; (deleted until per-atom info is output)
837  delete [] atomNames;
838  delete [] atomSegResids;
839 #endif
840 }
841 
842 
844 {
845 #ifndef MEM_OPT_VERSION
846  Bond *bonds = g_mol->getAllBonds();
847 
848  //then creating bond's tupleSignature
849  for(int i=0; i<g_mol->numBonds; i++)
850  {
851  Bond *b = bonds+i;
852  TupleSignature oneSig(1,BOND,b->bond_type);
853  oneSig.offset[0] = b->atom2 - b->atom1;
854  oneSig.isReal = (i<g_mol->numRealBonds);
855 
856  int poolIndex = sigsOfBonds.lookupCstPool(oneSig);
857  int newSig=0;
858  if(poolIndex == -1)
859  {
860  sigsOfBonds.push_back(oneSig);
861  poolIndex = (SigIndex)sigsOfBonds.size()-1;
862  newSig=1;
863  }
864 
865  if(!newSig)
866  {//check duplicate bonds in the form of (a, b) && (a, b);
867  int dupIdx = lookupCstPool(eachAtomSigs[b->atom1].bondSigIndices, (SigIndex)poolIndex);
868  if(dupIdx!=-1)
869  {
870  char err_msg[128];
871  sprintf(err_msg, "Duplicate bond %d-%d!", b->atom1+1, b->atom2+1);
872  NAMD_die(err_msg);
873  }
874  }
875  eachAtomSigs[b->atom1].bondSigIndices.push_back(poolIndex);
876  }
877 
878  //check duplicate bonds in the form of (a, b) && (b, a)
879  for(int i=0; i<g_mol->numBonds; i++)
880  {
881  Bond *b=bonds+i;
882  int atom2 = b->atom2;
883  int thisOffset = atom2 - b->atom1;
884  for(int j=0; j<eachAtomSigs[atom2].bondSigIndices.size(); j++)
885  {
886  SigIndex atom2BondId = eachAtomSigs[atom2].bondSigIndices[j];
887  TupleSignature *secSig = &(sigsOfBonds[atom2BondId]);
888  if(thisOffset== -(secSig->offset[0]))
889  {
890  char err_msg[128];
891  sprintf(err_msg, "Duplicate bond %d-%d because two atoms are just reversed!", b->atom1+1, atom2+1);
892  NAMD_die(err_msg);
893  }
894  }
895  }
896 
897  //building clusters for this simulation system in two steps
898  //1. create a list for each atom where each atom in the list is bonded with that atom
899  vector<int> *atomListOfBonded = new vector<int>[g_mol->numAtoms];
900 
901  for(int i=0; i<g_mol->numRealBonds; i++)
902  {
903  Bond *b=bonds+i;
904  int atom1 = b->atom1;
905  int atom2 = b->atom2;
906  atomListOfBonded[atom1].push_back(atom2);
907  atomListOfBonded[atom2].push_back(atom1);
908  }
909 
910  delete [] bonds;
911 
912  //2. using breadth-first-search to build the clusters. Here, we avoid recursive call
913  // because the depth of calls may be of thousands which will blow up the stack, and
914  //recursive call is slower than the stack-based BFS.
915  //Considering such structure
916  //1->1245; 7->1243; 1243->1245
917  eachAtomClusterID = new int[g_mol->numAtoms];
918  for(int i=0; i<g_mol->numAtoms; i++)
919  eachAtomClusterID[i] = -1;
920 
921  //It is guaranteed that the clusters found in this way use the
922  //smallest atom id of this cluster because each atom is at least
923  //connected to one of the atoms in its cluster (by atomListOfBonded
924  //constructed above).
925  //--Chao Mei
926  for(int i=0; i<g_mol->numAtoms; i++)
927  {
928  int curClusterID=eachAtomClusterID[i];
929  //if the atom's cluster id is not -1, the atom has been visited
930  if(curClusterID!=-1) continue;
931 
932  curClusterID=i;
933  deque<int> toVisitAtoms;
934  eachAtomClusterID[i] = curClusterID;
935  toVisitAtoms.push_back(i);
936  while(!toVisitAtoms.empty())
937  {
938  int visAtomID = toVisitAtoms.front();
939  toVisitAtoms.pop_front();
940  for(int j=0; j<atomListOfBonded[visAtomID].size(); j++)
941  {
942  int otherAtom = atomListOfBonded[visAtomID][j];
943  if(eachAtomClusterID[otherAtom]!=curClusterID){
944  eachAtomClusterID[otherAtom]=curClusterID;
945  toVisitAtoms.push_back(otherAtom);
946  }
947  }
948  }
949  }
950 
951 #if 0
952  //Now the clusterID of each atom should be usually in the non-decreasing
953  //order. In other words, the atom ids of a cluster are generally contiguous.
954  //If this is the case, the temporary memory usage of output IO during
955  //the simulation can be dramatically reduced. So g_isClusterContiguous
956  //is used to differentiate the two cases:
957 
958  //1. if the cluster id of atoms is monotonically increasing
959  //(g_isClusterContiguous=1), the size of the cluster can be used as this
960  //this cluster's signature (represented by "eachAtomClusterID").
961  //The atom who is the first atom (in terms of atom id) in this cluster will
962  //store the cluster size as its signature. The remaining atoms in this
963  //cluster store -1.
964 
965  //2. if the cluster id of atoms is not monotonically increasing, that is,
966  //the atom ids of a cluster are not contiguous (g_isClusterContiguous=0).
967  //Then we have to still record each atom's cluster id in variable
968  //"eachAtomClusterID", and we have to change each unique cluster id (valued
969  //at the scale of numAtoms) into indexes at the scale of numClusters. For
970  //example, the cluster ids of atoms up to this point may be (0...0, 24...24,
971  //789...789,...), after re-scaling, cluster ids should be (0...0,1...1,
972  //2...2,....) for atoms.
973  //We made an ASSUMPTION here: the cluster ids are in the non-decreasing
974  //order, and when a cluster discontiguity happens, the cluster id has
975  //appeared before! Such ASSUMPTION is to make implementation easy.
976 
977  int curClusterID;
978  int prevClusterID = eachAtomClusterID[0];
979  int curClusterSize = 1;
980  //step1: segment all atoms according to each cluster
981  for(int i=1; i<g_mol->numAtoms; i++){
982  curClusterID = eachAtomClusterID[i];
983  if(curClusterID == prevClusterID){
984  curClusterSize++;
985  }else{
986  eachClusterSize.push_back(curClusterSize);
987  eachClusterID.push_back(prevClusterID);
988  curClusterSize=1;
989  }
990  prevClusterID = curClusterID;
991  }
992  //record the info of the last cluster segment
993  eachClusterSize.push_back(curClusterSize);
994  eachClusterID.push_back(prevClusterID);
995 
996  //step2: detect contiguity of atoms in each cluster and re-scale
997  //cluster id.
998  g_isClusterContiguous = 1;
999  int *newClusterIDs = new int[eachClusterID.size()];
1000  memset(newClusterIDs, 0, sizeof(int)*eachClusterID.size());
1001  prevClusterID = eachClusterID[0];
1002  int newCId = 0;
1003  newClusterIDs[0] = newCId;
1004  for(int seg=1; seg<eachClusterID.size(); seg++){
1005  curClusterID = eachClusterID[seg];
1006  if(curClusterID > prevClusterID){
1007  newClusterIDs[seg] = ++newCId;
1008  prevClusterID = curClusterID;
1009  }else{
1010  //non-contiguity happens
1011  g_isClusterContiguous = 0;
1012 
1013  //we ASSUME this id appears before
1014  //binary search in eachAtomClusterID[0...seg-1).
1015  int jl=0, jh=seg-2;
1016  int isFound = 0;
1017  while(jh>=jl){
1018  int mid = (jl+jh)/2;
1019  if(curClusterID > eachClusterID[mid])
1020  jl = mid+1;
1021  else if(curClusterID < eachClusterID[mid])
1022  jh = mid-1;
1023  else{
1024  newClusterIDs[seg] = newClusterIDs[mid];
1025  isFound = 1;
1026  break;
1027  }
1028  }
1029  if(!isFound){
1030  //ASSUMPTION is wrong and abort
1031  char errmsg[300];
1032  sprintf(errmsg, "Assumption about building cluster is broken in file %s at line %d\n", __FILE__, __LINE__);
1033  NAMD_die(errmsg);
1034  }
1035  }
1036  }
1037 
1038  //step 3: modify eachAtomClusterID according to g_isClusterContiguous
1039  //newCId is the id of the last cluster, as id starts from 0, so the
1040  //total number clusters should be newCId+1
1041  g_numClusters = newCId+1;
1042  if(g_isClusterContiguous){
1043  int aid=0;
1044  for(int seg=0; seg<eachClusterSize.size(); seg++)
1045  {
1046  int curSize = eachClusterSize[seg];
1047  eachAtomClusterID[aid] = curSize;
1048  for(int i=aid+1; i<aid+curSize; i++)
1049  eachAtomClusterID[i] = -1;
1050  aid += curSize;
1051  }
1052  }else{
1053  int aid=0;
1054  for(int seg=0; seg<eachClusterSize.size(); seg++)
1055  {
1056  int curSize = eachClusterSize[seg];
1057  for(int i=aid; i<aid+curSize; i++)
1058  eachAtomClusterID[i] = newClusterIDs[seg];
1059  aid += curSize;
1060  }
1061  }
1062  free(newClusterIDs);
1063  eachClusterSize.clear();
1064  eachClusterID.clear();
1065 #endif
1066 
1067 /*
1068  //check whether cluster is built correctly
1069  printf("num clusters: %d\n", g_numClusters);
1070  FILE *checkFile = fopen("cluster.opt", "w");
1071  for(int i=0; i<g_mol->numAtoms; i++) fprintf(checkFile, "%d\n", eachAtomClusterID[i]);
1072  fclose(checkFile);
1073 */
1074 
1075  for(int i=0; i<g_mol->numAtoms; i++)
1076  atomListOfBonded[i].clear();
1077  delete [] atomListOfBonded;
1078 #endif
1079 }
1080 
1082 {
1083 #ifndef MEM_OPT_VERSION
1084  Angle *angles = g_mol->getAllAngles();
1085  //create angles' tupleSignature
1086  for(int i=0; i<g_mol->numAngles; i++)
1087  {
1088  Angle *tuple = angles+i;
1089  TupleSignature oneSig(2,ANGLE,tuple->angle_type);
1090  int offset[2];
1091  offset[0] = tuple->atom2 - tuple->atom1;
1092  offset[1] = tuple->atom3 - tuple->atom1;
1093  oneSig.setOffsets(offset);
1094 
1095  int poolIndex = sigsOfAngles.lookupCstPool(oneSig);
1096  if(poolIndex == -1)
1097  {
1098  sigsOfAngles.push_back(oneSig);
1099  poolIndex = (SigIndex)sigsOfAngles.size()-1;
1100  }
1101  eachAtomSigs[tuple->atom1].angleSigIndices.push_back(poolIndex);
1102  }
1103  delete [] angles;
1104 #endif
1105 }
1106 
1108 {
1109 #ifndef MEM_OPT_VERSION
1110  Dihedral *dihedrals = g_mol->getAllDihedrals();
1111 
1112  //create dihedrals' tupleSignature
1113  for(int i=0; i<g_mol->numDihedrals; i++)
1114  {
1115  Dihedral *tuple = dihedrals+i;
1116  TupleSignature oneSig(3,DIHEDRAL,tuple->dihedral_type);
1117  int offset[3];
1118  offset[0] = tuple->atom2 - tuple->atom1;
1119  offset[1] = tuple->atom3 - tuple->atom1;
1120  offset[2] = tuple->atom4 - tuple->atom1;
1121  oneSig.setOffsets(offset);
1122 
1123  int poolIndex = sigsOfDihedrals.lookupCstPool(oneSig);
1124  if(poolIndex == -1)
1125  {
1126  sigsOfDihedrals.push_back(oneSig);
1127  poolIndex = (SigIndex)sigsOfDihedrals.size()-1;
1128  }
1129  eachAtomSigs[tuple->atom1].dihedralSigIndices.push_back(poolIndex);
1130  }
1131 
1132  delete[] dihedrals;
1133 #endif
1134 }
1135 
1137 {
1138 #ifndef MEM_OPT_VERSION
1139  Improper *impropers=g_mol->getAllImpropers();
1140 
1141  //create improper's tupleSignature
1142  for(int i=0; i<g_mol->numImpropers; i++)
1143  {
1144  Improper *tuple = impropers+i;
1145  TupleSignature oneSig(3,IMPROPER,tuple->improper_type);
1146  int offset[3];
1147  offset[0] = tuple->atom2 - tuple->atom1;
1148  offset[1] = tuple->atom3 - tuple->atom1;
1149  offset[2] = tuple->atom4 - tuple->atom1;
1150  oneSig.setOffsets(offset);
1151 
1152  int poolIndex = sigsOfImpropers.lookupCstPool(oneSig);
1153  if(poolIndex == -1)
1154  {
1155  sigsOfImpropers.push_back(oneSig);
1156  poolIndex = (SigIndex)sigsOfImpropers.size()-1;
1157  }
1158  eachAtomSigs[tuple->atom1].improperSigIndices.push_back(poolIndex);
1159  }
1160 
1161  delete[] impropers;
1162 #endif
1163 }
1164 
1166 {
1167 #ifndef MEM_OPT_VERSION
1168  Crossterm *crossterms = g_mol->getAllCrossterms();
1169  //create crossterm's tupleSignature
1170  for(int i=0; i<g_mol->numCrossterms; i++)
1171  {
1172  Crossterm *tuple = crossterms+i;
1173  TupleSignature oneSig(7, CROSSTERM, tuple->crossterm_type);
1174  int offset[7];
1175  offset[0] = tuple->atom2 - tuple->atom1;
1176  offset[1] = tuple->atom3 - tuple->atom1;
1177  offset[2] = tuple->atom4 - tuple->atom1;
1178  offset[3] = tuple->atom5 - tuple->atom1;
1179  offset[4] = tuple->atom6 - tuple->atom1;
1180  offset[5] = tuple->atom7 - tuple->atom1;
1181  offset[6] = tuple->atom8 - tuple->atom1;
1182  oneSig.setOffsets(offset);
1183 
1184  int poolIndex = sigsOfCrossterms.lookupCstPool(oneSig);
1185  if(poolIndex == -1)
1186  {
1187  sigsOfCrossterms.push_back(oneSig);
1188  poolIndex = (SigIndex)sigsOfCrossterms.size()-1;
1189  }
1190  eachAtomSigs[tuple->atom1].crosstermSigIndices.push_back(poolIndex);
1191  }
1192 
1193  delete[] crossterms;
1194 #endif
1195 }
1196 
1198 {
1199  //1. Build exclusions: mainly accomplish the function of
1200  //Molecule::build_exclusions (based on the bonds)
1201  UniqueSet<Exclusion> allExclusions;
1202 
1203  int exclude_flag; //Exclusion policy
1204  exclude_flag = g_simParam->exclude;
1205  //int stripHGroupExclFlag = (simParams->splitPatch == SPLIT_PATCH_HYDROGEN);
1206 
1207  //Commented now since no explicit exclusions are read
1208  // Go through the explicit exclusions and add them to the arrays
1209  //for(i=0; i<numExclusions; i++){
1210  // exclusionSet.add(exclusions[i]);
1211  //}
1212 
1213  // If this is AMBER force field, and readExclusions is TRUE,
1214  // then all the exclusions were read from parm file, and we
1215  // shouldn't generate any of them.
1216  // Comment on stripHGroupExcl:
1217  // 1. Inside this function, hydrogenGroup is initialized in
1218  // build_atom_status, therefore, not available when reading psf files
1219  // 2. this function's main purpose is to reduce memory usage. Since exclusion
1220  // signatures are used, this function could be overlooked --Chao Mei
1221 
1222  vector<int> *eachAtomNeighbors = new vector<int>[g_mol->numAtoms];
1223  for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
1224  {
1225  AtomSigInfo *aSig = &atomSigPool[atomData[atom1].atomSigIdx];
1226  for(int j=0; j<aSig->bondSigIndices.size(); j++)
1227  {
1228  TupleSignature *tSig = &sigsOfBonds[aSig->bondSigIndices[j]];
1229  if(!tSig->isReal) continue;
1230  int atom2 = atom1+tSig->offset[0];
1231  eachAtomNeighbors[atom1].push_back(atom2);
1232  eachAtomNeighbors[atom2].push_back(atom1);
1233  }
1234  }
1235 
1237  { // Now calculate the bonded exlcusions based on the exclusion policy
1238  switch (exclude_flag)
1239  {
1240  case NONE:
1241  break;
1242  case ONETWO:
1243  build12Excls(allExclusions, eachAtomNeighbors);
1244  break;
1245  case ONETHREE:
1246  build12Excls(allExclusions, eachAtomNeighbors);
1247  build13Excls(allExclusions, eachAtomNeighbors);
1248  //if ( stripHGroupExclFlag ) stripHGroupExcl();
1249  break;
1250  case ONEFOUR:
1251  build12Excls(allExclusions, eachAtomNeighbors);
1252  build13Excls(allExclusions, eachAtomNeighbors);
1253  build14Excls(allExclusions, eachAtomNeighbors, 0);
1254  //if ( stripHGroupExclFlag ) stripHGroupExcl();
1255  break;
1256  case SCALED14:
1257  build12Excls(allExclusions, eachAtomNeighbors);
1258  build13Excls(allExclusions, eachAtomNeighbors);
1259  build14Excls(allExclusions, eachAtomNeighbors, 1);
1260  //if ( stripHGroupExclFlag ) stripHGroupExcl();
1261  break;
1262  }
1263  }
1264  //else if (stripHGroupExclFlag && exclude_flag!=NONE && exclude_flag!=ONETWO)
1265  // stripHGroupExcl();
1266 
1267  //Commented since atomFepFlags information is not available when reading psf file
1268  //stripFepExcl();
1269 
1270  for(int i=0; i<g_mol->numAtoms; i++)
1271  eachAtomNeighbors[i].clear();
1272  delete [] eachAtomNeighbors;
1273 
1274  //2. Build each atom's list of exclusions
1275  iout << iINFO << "ADDED " << allExclusions.size() << " IMPLICIT EXCLUSIONS\n" << endi;
1276  UniqueSetIter<Exclusion> exclIter(allExclusions);
1278  for(exclIter=exclIter.begin(); exclIter!=exclIter.end(); exclIter++)
1279  {
1280  int atom1 = exclIter->atom1;
1281  int atom2 = exclIter->atom2;
1282  int offset21 = atom2-atom1;
1283  if(exclIter->modified)
1284  {
1285  eachAtomExclSigs[atom1].modExclOffset.push_back(offset21);
1286  eachAtomExclSigs[atom2].modExclOffset.push_back(-offset21);
1287  }
1288  else
1289  {
1290  eachAtomExclSigs[atom1].fullExclOffset.push_back(offset21);
1291  eachAtomExclSigs[atom2].fullExclOffset.push_back(-offset21);
1292  }
1293  }
1294  allExclusions.clear();
1295 
1296  //3. Build up exclusion signatures and determine each atom's
1297  //exclusion signature index
1298  for(int i=0; i<g_mol->numAtoms; i++)
1299  {
1301  int poolIndex = sigsOfExclusions.lookupCstPool(eachAtomExclSigs[i]);
1302  if(poolIndex==-1)
1303  {
1304  poolIndex = sigsOfExclusions.size();
1305  sigsOfExclusions.push_back(eachAtomExclSigs[i]);
1306  }
1307  atomData[i].exclSigIdx = poolIndex;
1308  }
1309  delete [] eachAtomExclSigs;
1310  eachAtomExclSigs = NULL;
1311  printf("Exclusion signatures: %d\n", (int)sigsOfExclusions.size());
1312 }
1313 
1314 void build12Excls(UniqueSet<Exclusion>& allExcls, vector<int> *eachAtomNeighbors)
1315 {
1316  for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
1317  {
1318  vector<int> *atom1List = &eachAtomNeighbors[atom1];
1319  for(int j=0; j<atom1List->size(); j++)
1320  {
1321  int atom2 = atom1List->at(j);
1322  if(atom1<atom2)
1323  allExcls.add(Exclusion(atom1, atom2));
1324  else
1325  allExcls.add(Exclusion(atom2, atom1));
1326  }
1327  }
1328 }
1329 
1330 void build13Excls(UniqueSet<Exclusion>& allExcls, vector<int> *eachAtomNeighbors)
1331 {
1332  for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
1333  {
1334  vector<int> *atom1List = &eachAtomNeighbors[atom1];
1335  for(int j=0; j<atom1List->size(); j++)
1336  {
1337  int atom2 = atom1List->at(j);
1338  vector<int> *atom2List = &eachAtomNeighbors[atom2];
1339  for(int k=0; k<atom2List->size(); k++)
1340  {
1341  int atom3 = atom2List->at(k);
1342  //atom1-atom2, so atom2List contains atom1 which should not be considered
1343  if(atom3 == atom1)
1344  continue;
1345  if(atom1<atom3)
1346  allExcls.add(Exclusion(atom1, atom3));
1347  else
1348  allExcls.add(Exclusion(atom3, atom1));
1349  }
1350  }
1351  }
1352 }
1353 
1354 void build14Excls(UniqueSet<Exclusion>& allExcls, vector<int> *eachAtomNeighbors, int modified)
1355 {
1356  for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
1357  {
1358  vector<int> *atom1List = &eachAtomNeighbors[atom1];
1359  for(int j=0; j<atom1List->size(); j++)
1360  {
1361  int atom2 = atom1List->at(j);
1362  vector<int> *atom2List = &eachAtomNeighbors[atom2];
1363  for(int k=0; k<atom2List->size(); k++)
1364  {
1365  int atom3 = atom2List->at(k);
1366  //atom1-atom2, so atom2List contains atom1 which should not be considered
1367  if(atom3 == atom1)
1368  continue;
1369  vector<int> *atom3List = &eachAtomNeighbors[atom3];
1370  for(int l=0; l<atom3List->size(); l++)
1371  {
1372  int atom4 = atom3List->at(l);
1373  //atom1-atom2, so atom2List contains atom1 which should not be considered
1374  if(atom4 == atom2 || atom4 == atom1)
1375  continue;
1376  if(atom1<atom4)
1377  allExcls.add(Exclusion(atom1, atom4, modified));
1378  else
1379  allExcls.add(Exclusion(atom4, atom1, modified));
1380  }
1381  }
1382  }
1383  }
1384 }
1385 
1386 template <class T> void HashPool<T>::dump_tables()
1387 {
1388  for(int j=0; j < pool.size(); j++) {
1389  HashPoolAdaptorT<T>* pval = pool[j];
1390  CmiPrintf("Pool[%d]=%p %p hash = %d\n",j,pool[j],pval,pval->hash());
1391  }
1392  CkHashtableIterator *iter = index_table.iterator();
1393  void *key,*indx;
1394  while (iter->hasNext()) {
1395  indx = iter->next(&key);
1397  CmiPrintf("key %p indx %p %d hash=%d\n",key,indx,*((int *)indx),pkey->hash());
1398  }
1399 }
AtomSigInfo(const AtomSigInfo &sig)
Definition: CompressPsf.C:66
#define SCALED14
Definition: SimParameters.h:43
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:81
Index crossterm_type
Definition: structures.h:89
CkHashCode hash(void) const
Definition: CompressPsf.h:77
vector< Bond > extraBonds
Definition: CompressPsf.C:333
int32 atom4
Definition: structures.h:75
Index improper_type
Definition: structures.h:76
void flipNum(char *elem, int elemSize, int numElems)
Definition: CompressPsf.C:406
int numBonds
Definition: Molecule.h:560
#define COMPRESSED_PSF_MAGICNUM
Definition: CompressPsf.h:13
vector< DihedralValue > extraDihedralParams
Definition: CompressPsf.C:340
vector< Angle > extraAngles
Definition: CompressPsf.C:334
void clear(void)
Definition: UniqueSet.h:62
int numHydrogenGroups
Definition: Molecule.h:600
int32 atom3
Definition: structures.h:65
void buildExclusions()
Definition: CompressPsf.C:1197
void buildDihedralData()
Definition: CompressPsf.C:1107
int hash() const
Definition: CompressPsf.C:276
void buildAtomData()
Definition: CompressPsf.C:764
Real r_ub
Definition: Parameters.h:96
Index resNameIdx
Definition: CompressPsf.C:40
static __thread atom * atoms
Atom * getAtoms() const
Definition: Molecule.h:491
float Real
Definition: common.h:109
vector< Dihedral > extraDihedrals
Definition: CompressPsf.C:335
int32 atom5
Definition: structures.h:85
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
void build14Excls(UniqueSet< Exclusion > &, vector< int > *, int)
Definition: CompressPsf.C:1354
Real rigid_bond_length(int atomnum) const
Definition: Molecule.h:1457
BasicAtomInfo * atomData
Definition: CompressPsf.C:314
int32 atom8
Definition: structures.h:88
void outputCompressedFile(FILE *txtOfp, FILE *binOfp)
Definition: CompressPsf.C:541
int numRealBonds
Definition: Molecule.h:559
Index atomTypeIdx
Definition: CompressPsf.C:42
Crossterm * getAllCrossterms() const
Definition: Molecule.h:1098
#define ONETHREE
Definition: SimParameters.h:41
DihedralValue * dihedral_array
Definition: Parameters.h:245
int32 atom1
Definition: structures.h:81
int add(const Elem &elem)
Definition: UniqueSet.h:52
int size(void) const
Definition: UniqueSet.h:58
#define iout
Definition: InfoStream.h:51
vector< SigIndex > crosstermSigIndices
Definition: CompressPsf.C:62
vector< SigIndex > improperSigIndices
Definition: CompressPsf.C:61
int32 atom4
Definition: structures.h:84
HashPool< TupleSignature > sigsOfAngles
Definition: CompressPsf.C:323
SimParameters * g_simParam
Definition: CompressPsf.C:34
vector< AngleValue > extraAngleParams
Definition: CompressPsf.C:339
int NumDihedralParams
Definition: Parameters.h:264
ExclSigInfo(const ExclSigInfo &sig)
Definition: CompressPsf.C:205
int32 atom3
Definition: structures.h:83
Index atomNameIdx
Definition: CompressPsf.C:41
int atomsInGroup
Definition: Hydrogen.h:19
void loadMolInfo()
Definition: CompressPsf.C:473
#define COMPRESSED_PSF_VER
Definition: CompressPsf.h:9
HashPool< TupleSignature > sigsOfBonds
Definition: CompressPsf.C:322
vector< ImproperValue > extraImproperParams
Definition: CompressPsf.C:341
UniqueSetIter< T > begin(void) const
Definition: UniqueSetIter.h:55
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.h:110
int32 atom4
Definition: structures.h:66
HydrogenGroup hydrogenGroup
Definition: Molecule.h:640
int32 atom1
Definition: structures.h:48
int Index
Definition: structures.h:26
void integrateAllAtomSigs()
Definition: CompressPsf.C:505
void setOffsets(int *offs)
Definition: structures.h:262
void buildExclusionData()
#define MAX_MULTIPLICITY
Definition: Parameters.h:70
vector< int > fullExclOffset
Definition: CompressPsf.C:200
int hash() const
Definition: CompressPsf.C:228
HashPool< ExclSigInfo > sigsOfExclusions
Definition: CompressPsf.C:329
vector< BondValue > extraBondParams
Definition: CompressPsf.C:338
const float * getBFactorData()
Definition: Molecule.h:1036
int lookupCstPool(const std::vector< T > &pool, const T &val)
Definition: CompressPsf.h:56
#define ONETWO
Definition: SimParameters.h:40
int32 atom2
Definition: structures.h:82
Index dihedral_type
Definition: structures.h:67
int32 atom2
Definition: structures.h:56
int32 atom7
Definition: structures.h:87
Index atomSigIdx
Definition: CompressPsf.C:45
int NumImproperParams
Definition: Parameters.h:265
HashPool< HashReal > chargePool
Definition: CompressPsf.C:311
vector< SigIndex > bondSigIndices
Definition: CompressPsf.C:58
int32 atom1
Definition: structures.h:55
Index segNameIdx
Definition: CompressPsf.C:39
unsigned int circShift(unsigned int h, unsigned int by)
Definition: structures.h:186
int32 atom3
Definition: structures.h:57
int32 hydrogenList
Definition: structures.h:43
void buildImproperData()
Definition: CompressPsf.C:1136
AtomSigInfo * eachAtomSigs
Definition: CompressPsf.C:327
Index chargeIdx
Definition: CompressPsf.C:43
void clearGlobalVectors()
Definition: CompressPsf.C:419
vector< int > eachClusterID
Definition: CompressPsf.C:319
int numAngles
Definition: Molecule.h:561
const float * getOccupancyData()
Definition: Molecule.h:1032
int numAtoms
Definition: Molecule.h:557
ImproperValue * improper_array
Definition: Parameters.h:246
#define NONE
Definition: SimParameters.h:39
ExclSigInfo * eachAtomExclSigs
Definition: CompressPsf.C:330
int numCrossterms
Definition: Molecule.h:568
void NAMD_die(const char *err_msg)
Definition: common.C:85
Parameters * g_param
Definition: CompressPsf.C:33
HashPool< HashString > atomNamePool
Definition: CompressPsf.C:309
HashPool< TupleSignature > sigsOfImpropers
Definition: CompressPsf.C:325
HashPool< HashString > resNamePool
Definition: CompressPsf.C:308
int hash() const
Definition: CompressPsf.C:295
int operator!=(const FourBodyConsts &f1, const FourBodyConsts &f2)
Definition: CompressPsf.C:353
Index angle_type
Definition: structures.h:58
AtomNameInfo * getAtomNames() const
Definition: Molecule.h:492
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.h:116
int maxMigrationGroupSize
Definition: Molecule.h:603
int numDihedrals
Definition: Molecule.h:562
int32 atom2
Definition: structures.h:64
int numImpropers
Definition: Molecule.h:567
#define ONEFOUR
Definition: SimParameters.h:42
int32 atom1
Definition: structures.h:72
int g_numClusters
Definition: CompressPsf.C:320
void compress_molecule_info(Molecule *mol, char *psfFileName, Parameters *param, SimParameters *simParam, ConfigList *cfgList)
Definition: CompressPsf.C:436
struct OutputAtomRecord::floatVals fSet
UniqueSetIter< T > end(void) const
Definition: UniqueSetIter.h:64
BlockRadixSort::TempStorage sort
void build13Excls(UniqueSet< Exclusion > &, vector< int > *)
Definition: CompressPsf.C:1330
vector< SigIndex > angleSigIndices
Definition: CompressPsf.C:59
vector< int > eachClusterSize
Definition: CompressPsf.C:318
int SigIndex
Definition: NamdTypes.h:81
Dihedral * getAllDihedrals() const
Definition: Molecule.h:1097
int atomsInMigrationGroup
Definition: Hydrogen.h:29
AtomSegResInfo * getAtomSegResInfo() const
Definition: Molecule.h:495
HashPool< TupleSignature > sigsOfCrossterms
Definition: CompressPsf.C:326
Index vdw_type
Definition: structures.h:40
void buildAngleData()
Definition: CompressPsf.C:1081
void resize(int i)
Definition: ResizeArray.h:84
struct OutputAtomRecord::integerVals iSet
int maxHydrogenGroupSize
Definition: Molecule.h:601
void buildCrosstermData()
Definition: CompressPsf.C:1165
void dump_tables()
Definition: CompressPsf.C:1386
int numMigrationGroups
Definition: Molecule.h:602
Real x0
Definition: Parameters.h:87
void output(FILE *ofp)
Definition: structures.h:298
int32 atom2
Definition: structures.h:73
Index exclSigIdx
Definition: CompressPsf.C:46
vector< int > modExclOffset
Definition: CompressPsf.C:201
HashPool< HashString > segNamePool
Definition: CompressPsf.C:307
HashPool< HashString > atomTypePool
Definition: CompressPsf.C:310
int32 atom6
Definition: structures.h:86
int operator==(const AtomSigInfo &s1, const AtomSigInfo &s2)
Definition: CompressPsf.C:146
vector< Improper > extraImpropers
Definition: CompressPsf.C:336
HashPool< TupleSignature > sigsOfDihedrals
Definition: CompressPsf.C:324
void sortTupleSigIndices()
Definition: CompressPsf.C:98
int * eachAtomClusterID
Definition: CompressPsf.C:317
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
Angle * getAllAngles() const
Definition: Molecule.h:1095
Real theta0
Definition: Parameters.h:94
void buildBondData()
Definition: CompressPsf.C:843
vector< SigIndex > dihedralSigIndices
Definition: CompressPsf.C:60
void sortExclOffset()
Definition: CompressPsf.C:222
CkHashCode hash() const
Definition: CompressPsf.C:107
void getExtraBonds(StringList *file)
Index bond_type
Definition: structures.h:50
int resid
Definition: Molecule.h:147
ExclusionSettings exclude
void build12Excls(UniqueSet< Exclusion > &, vector< int > *)
Definition: CompressPsf.C:1314
HashPool< HashReal > massPool
Definition: CompressPsf.C:312
struct OutputAtomRecord::shortVals sSet
Bool is_water(int)
Bond * getAllBonds() const
Definition: Molecule.h:1094
int32 atom3
Definition: structures.h:74
Molecule * g_mol
Definition: CompressPsf.C:32
Real k_ub
Definition: Parameters.h:95
Real k
Definition: Parameters.h:86
int32 atom1
Definition: structures.h:63
Improper * getAllImpropers() const
Definition: Molecule.h:1096
int32 atom2
Definition: structures.h:49
ConfigList * g_cfgList
Definition: CompressPsf.C:35
HashReal(Real v)
Definition: CompressPsf.C:291
HashPool< AtomSigInfo > atomSigPool
Definition: CompressPsf.C:313
iterator begin(void)
Definition: ResizeArray.h:36