ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
(Generate patch)

Comparing branches/development/src/io/DumpWriter.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1752 by gezelter, Sun Jun 10 14:05:02 2012 UTC

# Line 36 | Line 36
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include "io/DumpWriter.hpp"
# Line 58 | Line 59 | namespace OpenMD {
59      : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
60  
61      Globals* simParams = info->getSimParams();
62 <    needCompression_ = simParams->getCompressDumpFile();
63 <    needForceVector_ = simParams->getOutputForceVector();
62 >    needCompression_   = simParams->getCompressDumpFile();
63 >    needForceVector_   = simParams->getOutputForceVector();
64 >    needParticlePot_   = simParams->getOutputParticlePotential();
65 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
66 >    needElectricField_ = simParams->getOutputElectricField();
67 >
68 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
69 >      doSiteData_ = true;
70 >    } else {
71 >      doSiteData_ = false;
72 >    }
73 >
74      createDumpFile_ = true;
75   #ifdef HAVE_LIBZ
76      if (needCompression_) {
# Line 97 | Line 108 | namespace OpenMD {
108      Globals* simParams = info->getSimParams();
109      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
110  
111 <    needCompression_ = simParams->getCompressDumpFile();
112 <    needForceVector_ = simParams->getOutputForceVector();
111 >    needCompression_   = simParams->getCompressDumpFile();
112 >    needForceVector_   = simParams->getOutputForceVector();
113 >    needParticlePot_   = simParams->getOutputParticlePotential();
114 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
115 >    needElectricField_ = simParams->getOutputElectricField();
116 >
117 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
118 >      doSiteData_ = true;
119 >    } else {
120 >      doSiteData_ = false;
121 >    }
122 >
123      createDumpFile_ = true;
124   #ifdef HAVE_LIBZ
125      if (needCompression_) {
# Line 136 | Line 157 | namespace OpenMD {
157      Globals* simParams = info->getSimParams();
158      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
159      
160 <    needCompression_ = simParams->getCompressDumpFile();
161 <    needForceVector_ = simParams->getOutputForceVector();
162 <    
160 >    needCompression_   = simParams->getCompressDumpFile();
161 >    needForceVector_   = simParams->getOutputForceVector();
162 >    needParticlePot_   = simParams->getOutputParticlePotential();
163 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
164 >    needElectricField_ = simParams->getOutputElectricField();
165 >
166 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
167 >      doSiteData_ = true;
168 >    } else {
169 >      doSiteData_ = false;
170 >    }
171 >
172   #ifdef HAVE_LIBZ
173      if (needCompression_) {
174        filename_ += ".gz";
# Line 272 | Line 302 | namespace OpenMD {
302      StuntDouble* integrableObject;
303      SimInfo::MoleculeIterator mi;
304      Molecule::IntegrableObjectIterator ii;
305 +    RigidBody::AtomIterator ai;
306 +    Atom* atom;
307  
308   #ifndef IS_MPI
309      os << "  <Snapshot>\n";
# Line 289 | Line 321 | namespace OpenMD {
321        }
322      }    
323      os << "    </StuntDoubles>\n";
324 <    
324 >
325 >    if (doSiteData_) {
326 >      os << "    <SiteData>\n";
327 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
328 >              
329 >        for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
330 >           integrableObject = mol->nextIntegrableObject(ii)) {  
331 >
332 >          int ioIndex = integrableObject->getGlobalIntegrableObjectIndex();
333 >          // do one for the IO itself
334 >          os << prepareSiteLine(integrableObject, ioIndex, 0);
335 >
336 >          if (integrableObject->isRigidBody()) {
337 >            
338 >            RigidBody* rb = static_cast<RigidBody*>(integrableObject);
339 >            int siteIndex = 0;
340 >            for (atom = rb->beginAtom(ai); atom != NULL;  
341 >                 atom = rb->nextAtom(ai)) {                                            
342 >              os << prepareSiteLine(atom, ioIndex, siteIndex);
343 >              siteIndex++;
344 >            }
345 >          }
346 >        }
347 >      }    
348 >      os << "    </SiteData>\n";
349 >    }
350      os << "  </Snapshot>\n";
351  
352      os.flush();
# Line 306 | Line 363 | namespace OpenMD {
363      }
364      
365      const int masterNode = 0;
366 <
366 >    int nProc;
367 >    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
368      if (worldRank == masterNode) {      
369        os << "  <Snapshot>\n";  
370        writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
# Line 314 | Line 372 | namespace OpenMD {
372          
373        os << buffer;
374  
317 –      int nProc;
318 –      MPI_Comm_size(MPI_COMM_WORLD, &nProc);
375        for (int i = 1; i < nProc; ++i) {
376  
377          // receive the length of the string buffer that was
378          // prepared by processor i
379  
380 +        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
381          int recvLength;
382          MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
383          char* recvBuffer = new char[recvLength];
# Line 337 | Line 394 | namespace OpenMD {
394        os.flush();
395      } else {
396        int sendBufferLength = buffer.size() + 1;
397 <      MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
398 <      MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
397 >      int myturn = 0;
398 >      for (int i = 1; i < nProc; ++i){
399 >        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
400 >        if (myturn == worldRank){
401 >          MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
402 >          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
403 >        }
404 >      }
405      }
406  
407   #endif // is_mpi
# Line 420 | Line 483 | namespace OpenMD {
483  
484      if (needForceVector_) {
485        type += "f";
486 <      Vector3d frc;
424 <
425 <      frc = integrableObject->getFrc();
426 <
486 >      Vector3d frc = integrableObject->getFrc();
487        if (isinf(frc[0]) || isnan(frc[0]) ||
488            isinf(frc[1]) || isnan(frc[1]) ||
489            isinf(frc[2]) || isnan(frc[2]) ) {      
# Line 439 | Line 499 | namespace OpenMD {
499        
500        if (integrableObject->isDirectional()) {
501          type += "t";
502 <        Vector3d trq;
443 <        
444 <        trq = integrableObject->getTrq();
445 <        
502 >        Vector3d trq = integrableObject->getTrq();        
503          if (isinf(trq[0]) || isnan(trq[0]) ||
504              isinf(trq[1]) || isnan(trq[1]) ||
505              isinf(trq[2]) || isnan(trq[2]) ) {      
# Line 451 | Line 508 | namespace OpenMD {
508                     " for object %d", index);      
509            painCave.isFatal = 1;
510            simError();
511 <        }
455 <        
511 >        }        
512          sprintf(tempBuffer, " %13e %13e %13e",
513                  trq[0], trq[1], trq[2]);
514          line += tempBuffer;
515        }      
516      }
517 <    
517 >
518      sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
519      return std::string(tempBuffer);
520    }
521  
522 +  std::string DumpWriter::prepareSiteLine(StuntDouble* integrableObject, int ioIndex, int siteIndex) {
523 +        
524 +
525 +    std::string id;
526 +    std::string type;
527 +    std::string line;
528 +    char tempBuffer[4096];
529 +
530 +    if (integrableObject->isRigidBody()) {
531 +      sprintf(tempBuffer, "%10d           ", ioIndex);
532 +      id = std::string(tempBuffer);
533 +    } else {
534 +      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
535 +      id = std::string(tempBuffer);
536 +    }
537 +              
538 +    if (needFlucQ_) {
539 +      type += "cw";
540 +      RealType fqPos = integrableObject->getFlucQPos();
541 +      if (isinf(fqPos) || isnan(fqPos) ) {      
542 +        sprintf( painCave.errMsg,
543 +                 "DumpWriter detected a numerical error writing the"
544 +                 " fluctuating charge for object %s", id.c_str());      
545 +        painCave.isFatal = 1;
546 +        simError();
547 +      }
548 +      sprintf(tempBuffer, " %13e ", fqPos);
549 +      line += tempBuffer;
550 +    
551 +      RealType fqVel = integrableObject->getFlucQVel();
552 +      if (isinf(fqVel) || isnan(fqVel) ) {      
553 +        sprintf( painCave.errMsg,
554 +                 "DumpWriter detected a numerical error writing the"
555 +                 " fluctuating charge velocity for object %s", id.c_str());      
556 +        painCave.isFatal = 1;
557 +        simError();
558 +      }
559 +      sprintf(tempBuffer, " %13e ", fqVel);
560 +      line += tempBuffer;
561 +
562 +      if (needForceVector_) {
563 +        type += "g";
564 +        RealType fqFrc = integrableObject->getFlucQFrc();        
565 +        if (isinf(fqFrc) || isnan(fqFrc) ) {      
566 +          sprintf( painCave.errMsg,
567 +                   "DumpWriter detected a numerical error writing the"
568 +                   " fluctuating charge force for object %s", id.c_str());      
569 +          painCave.isFatal = 1;
570 +          simError();
571 +        }
572 +        sprintf(tempBuffer, " %13e ", fqFrc);        
573 +        line += tempBuffer;
574 +      }
575 +    }
576 +
577 +    if (needElectricField_) {
578 +      type += "e";
579 +      Vector3d eField= integrableObject->getElectricField();
580 +      if (isinf(eField[0]) || isnan(eField[0]) ||
581 +          isinf(eField[1]) || isnan(eField[1]) ||
582 +          isinf(eField[2]) || isnan(eField[2]) ) {      
583 +        sprintf( painCave.errMsg,
584 +                 "DumpWriter detected a numerical error writing the electric"
585 +                 " field for object %s", id.c_str());      
586 +        painCave.isFatal = 1;
587 +        simError();
588 +      }
589 +      sprintf(tempBuffer, " %13e %13e %13e",
590 +              eField[0], eField[1], eField[2]);
591 +      line += tempBuffer;
592 +    }
593 +
594 +
595 +    if (needParticlePot_) {
596 +      type += "u";
597 +      RealType particlePot = integrableObject->getParticlePot();
598 +      if (isinf(particlePot) || isnan(particlePot)) {      
599 +        sprintf( painCave.errMsg,
600 +                 "DumpWriter detected a numerical error writing the particle "
601 +                 " potential for object %s", id.c_str());      
602 +        painCave.isFatal = 1;
603 +        simError();
604 +      }
605 +      sprintf(tempBuffer, " %13e", particlePot);
606 +      line += tempBuffer;
607 +    }
608 +    
609 +
610 +    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
611 +    return std::string(tempBuffer);
612 +  }
613 +
614    void DumpWriter::writeDump() {
615      writeFrame(*dumpFile_);
616    }
# Line 540 | Line 688 | namespace OpenMD {
688      newOStream = new std::ofstream(filename.c_str());
689   #endif
690      //write out MetaData first
691 <    (*newOStream) << "<OpenMD version=1>" << std::endl;
691 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
692      (*newOStream) << "  <MetaData>" << std::endl;
693      (*newOStream) << info_->getRawMetaData();
694      (*newOStream) << "  </MetaData>" << std::endl;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines