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

Comparing trunk/src/io/DumpWriter.cpp (file contents):
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 1993 by gezelter, Tue Apr 29 17:32:31 2014 UTC

# Line 35 | Line 35
35   *                                                                      
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
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 "config.h"
44 +
45 + #ifdef IS_MPI
46 + #include <mpi.h>
47 + #endif
48  
49   #include "io/DumpWriter.hpp"
50   #include "primitives/Molecule.hpp"
# Line 54 | Line 60
60   #define isinf(x) (!_finite(x) && !_isnan(x))
61   #endif
62  
57 – #ifdef IS_MPI
58 – #include <mpi.h>
59 – #endif
60 –
63   using namespace std;
64   namespace OpenMD {
65  
66    DumpWriter::DumpWriter(SimInfo* info)
67 <    : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
67 >    : info_(info), filename_(info->getDumpFileName()),
68 >      eorFilename_(info->getFinalConfigFileName()){
69  
70      Globals* simParams = info->getSimParams();
71      needCompression_   = simParams->getCompressDumpFile();
# Line 70 | Line 73 | namespace OpenMD {
73      needParticlePot_   = simParams->getOutputParticlePotential();
74      needFlucQ_         = simParams->getOutputFluctuatingCharges();
75      needElectricField_ = simParams->getOutputElectricField();
76 +    needSitePotential_ = simParams->getOutputSitePotential();
77  
78 <    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
78 >    if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
79 >        needSitePotential_) {
80        doSiteData_ = true;
81      } else {
82        doSiteData_ = false;
# Line 119 | Line 124 | namespace OpenMD {
124      needParticlePot_   = simParams->getOutputParticlePotential();
125      needFlucQ_         = simParams->getOutputFluctuatingCharges();
126      needElectricField_ = simParams->getOutputElectricField();
127 +    needSitePotential_ = simParams->getOutputSitePotential();
128  
129 <    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
129 >    if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
130 >        needSitePotential_) {
131        doSiteData_ = true;
132      } else {
133        doSiteData_ = false;
# Line 168 | Line 175 | namespace OpenMD {
175      needParticlePot_   = simParams->getOutputParticlePotential();
176      needFlucQ_         = simParams->getOutputFluctuatingCharges();
177      needElectricField_ = simParams->getOutputElectricField();
178 +    needSitePotential_ = simParams->getOutputSitePotential();
179  
180 <    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
180 >    if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
181 >        needSitePotential_) {
182        doSiteData_ = true;
183      } else {
184        doSiteData_ = false;
# Line 310 | Line 319 | namespace OpenMD {
319      SimInfo::MoleculeIterator mi;
320      Molecule::IntegrableObjectIterator ii;
321      RigidBody::AtomIterator ai;
313 –    Atom* atom;
322  
323   #ifndef IS_MPI
324      os << "  <Snapshot>\n";
# Line 318 | Line 326 | namespace OpenMD {
326      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
327  
328      os << "    <StuntDoubles>\n";
329 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
330 <
329 >    for (mol = info_->beginMolecule(mi); mol != NULL;
330 >         mol = info_->nextMolecule(mi)) {
331        
332        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
333             sd = mol->nextIntegrableObject(ii)) {        
# Line 331 | Line 339 | namespace OpenMD {
339  
340      if (doSiteData_) {
341        os << "    <SiteData>\n";
342 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
342 >      for (mol = info_->beginMolecule(mi); mol != NULL;
343 >           mol = info_->nextMolecule(mi)) {
344                
345          for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
346             sd = mol->nextIntegrableObject(ii)) {        
# Line 344 | Line 353 | namespace OpenMD {
353              
354              RigidBody* rb = static_cast<RigidBody*>(sd);
355              int siteIndex = 0;
356 <            for (atom = rb->beginAtom(ai); atom != NULL;  
356 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
357                   atom = rb->nextAtom(ai)) {                                            
358                os << prepareSiteLine(atom, ioIndex, siteIndex);
359                siteIndex++;
# Line 358 | Line 367 | namespace OpenMD {
367  
368      os.flush();
369   #else
361 –    //every node prepares the dump lines for integrable objects belong to itself
362 –    std::string buffer;
363 –    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
370  
371 +    const int masterNode = 0;
372 +    int worldRank;
373 +    int nProc;
374 +
375 +    MPI_Comm_size( MPI_COMM_WORLD, &nProc);
376 +    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
377 +
378 +
379 +    if (worldRank == masterNode) {      
380 +      os << "  <Snapshot>\n";  
381 +      writeFrameProperties(os,
382 +                           info_->getSnapshotManager()->getCurrentSnapshot());
383 +      os << "    <StuntDoubles>\n";
384 +    }
385  
386 +    //every node prepares the dump lines for integrable objects belong to itself
387 +    std::string buffer;
388 +    for (mol = info_->beginMolecule(mi); mol != NULL;
389 +         mol = info_->nextMolecule(mi)) {
390        for (sd = mol->beginIntegrableObject(ii); sd != NULL;
391             sd = mol->nextIntegrableObject(ii)) {        
392 <          buffer += prepareDumpLine(sd);
392 >        buffer += prepareDumpLine(sd);
393        }
394      }
395      
372 –    const int masterNode = 0;
373 –    int nProc;
374 –    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
396      if (worldRank == masterNode) {      
376 –      os << "  <Snapshot>\n";  
377 –      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
378 –      os << "    <StuntDoubles>\n";
379 –        
397        os << buffer;
398 <
398 >      
399        for (int i = 1; i < nProc; ++i) {
400 +        // tell processor i to start sending us data:
401 +        MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
402  
403          // receive the length of the string buffer that was
404 <        // prepared by processor i
386 <
387 <        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
404 >        // prepared by processor i:        
405          int recvLength;
406 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
406 >        MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
407 >                 &istatus);
408 >
409 >        // create a buffer to receive the data
410          char* recvBuffer = new char[recvLength];
411          if (recvBuffer == NULL) {
412          } else {
413 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
413 >          // receive the data:
414 >          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
415 >                               MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
416 >          // send it to the file:
417            os << recvBuffer;
418 +          // get rid of the receive buffer:
419            delete [] recvBuffer;
420          }
421        }
398 –      os << "    </StuntDoubles>\n";
399 –      
400 –      os << "  </Snapshot>\n";
401 –      os.flush();
422      } else {
423        int sendBufferLength = buffer.size() + 1;
424        int myturn = 0;
425        for (int i = 1; i < nProc; ++i){
426 <        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
426 >        // wait for the master node to call our number:
427 >        MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
428          if (myturn == worldRank){
429 +          // send the length of our buffer:
430            MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
431 <          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
431 >
432 >          // send our buffer:
433 >          MPI_Send((void *)buffer.c_str(), sendBufferLength,
434 >                   MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
435 >
436          }
437        }
438      }
439 +    
440 +    if (worldRank == masterNode) {      
441 +      os << "    </StuntDoubles>\n";
442 +    }
443  
444 < #endif // is_mpi
444 >    if (doSiteData_) {
445 >      if (worldRank == masterNode) {
446 >        os << "    <SiteData>\n";
447 >      }
448 >      buffer.clear();
449 >      for (mol = info_->beginMolecule(mi); mol != NULL;
450 >           mol = info_->nextMolecule(mi)) {
451 >              
452 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
453 >             sd = mol->nextIntegrableObject(ii)) {      
454 >          
455 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
456 >          // do one for the IO itself
457 >          buffer += prepareSiteLine(sd, ioIndex, 0);
458  
459 +          if (sd->isRigidBody()) {
460 +            
461 +            RigidBody* rb = static_cast<RigidBody*>(sd);
462 +            int siteIndex = 0;
463 +            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
464 +                 atom = rb->nextAtom(ai)) {                                            
465 +              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
466 +              siteIndex++;
467 +            }
468 +          }
469 +        }
470 +      }
471 +
472 +      if (worldRank == masterNode) {    
473 +        os << buffer;
474 +        
475 +        for (int i = 1; i < nProc; ++i) {
476 +          
477 +          // tell processor i to start sending us data:
478 +          MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
479 +          
480 +          // receive the length of the string buffer that was
481 +          // prepared by processor i:        
482 +          int recvLength;
483 +          MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
484 +                   &istatus);
485 +          
486 +          // create a buffer to receive the data
487 +          char* recvBuffer = new char[recvLength];
488 +          if (recvBuffer == NULL) {
489 +          } else {
490 +            // receive the data:
491 +            MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
492 +                     MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
493 +            // send it to the file:
494 +            os << recvBuffer;
495 +            // get rid of the receive buffer:
496 +            delete [] recvBuffer;
497 +          }
498 +        }      
499 +      } else {
500 +        int sendBufferLength = buffer.size() + 1;
501 +        int myturn = 0;
502 +        for (int i = 1; i < nProc; ++i){
503 +          // wait for the master node to call our number:
504 +          MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
505 +          if (myturn == worldRank){
506 +            // send the length of our buffer:
507 +            MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
508 +            // send our buffer:
509 +            MPI_Send((void *)buffer.c_str(), sendBufferLength,
510 +                     MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
511 +          }
512 +        }
513 +      }
514 +      
515 +      if (worldRank == masterNode) {    
516 +        os << "    </SiteData>\n";
517 +      }
518 +    }
519 +    
520 +    if (worldRank == masterNode) {
521 +      os << "  </Snapshot>\n";
522 +      os.flush();
523 +    }
524 +    
525 + #endif // is_mpi
526 +    
527    }
528  
529    std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
# Line 527 | Line 638 | namespace OpenMD {
638    }
639  
640    std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
641 <        
641 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
642  
643      std::string id;
644      std::string type;
# Line 543 | Line 654 | namespace OpenMD {
654      }
655                
656      if (needFlucQ_) {
657 <      type += "cw";
658 <      RealType fqPos = sd->getFlucQPos();
659 <      if (isinf(fqPos) || isnan(fqPos) ) {      
660 <        sprintf( painCave.errMsg,
661 <                 "DumpWriter detected a numerical error writing the"
662 <                 " fluctuating charge for object %s", id.c_str());      
663 <        painCave.isFatal = 1;
664 <        simError();
665 <      }
666 <      sprintf(tempBuffer, " %13e ", fqPos);
667 <      line += tempBuffer;
668 <    
669 <      RealType fqVel = sd->getFlucQVel();
559 <      if (isinf(fqVel) || isnan(fqVel) ) {      
560 <        sprintf( painCave.errMsg,
561 <                 "DumpWriter detected a numerical error writing the"
562 <                 " fluctuating charge velocity for object %s", id.c_str());      
563 <        painCave.isFatal = 1;
564 <        simError();
565 <      }
566 <      sprintf(tempBuffer, " %13e ", fqVel);
567 <      line += tempBuffer;
657 >      if (storageLayout & DataStorage::dslFlucQPosition) {
658 >        type += "c";
659 >        RealType fqPos = sd->getFlucQPos();
660 >        if (isinf(fqPos) || isnan(fqPos) ) {      
661 >          sprintf( painCave.errMsg,
662 >                   "DumpWriter detected a numerical error writing the"
663 >                   " fluctuating charge for object %s", id.c_str());      
664 >          painCave.isFatal = 1;
665 >          simError();
666 >        }
667 >        sprintf(tempBuffer, " %13e ", fqPos);
668 >        line += tempBuffer;
669 >      }
670  
671 <      if (needForceVector_) {
672 <        type += "g";
673 <        RealType fqFrc = sd->getFlucQFrc();        
674 <        if (isinf(fqFrc) || isnan(fqFrc) ) {      
671 >      if (storageLayout & DataStorage::dslFlucQVelocity) {
672 >        type += "w";    
673 >        RealType fqVel = sd->getFlucQVel();
674 >        if (isinf(fqVel) || isnan(fqVel) ) {      
675            sprintf( painCave.errMsg,
676                     "DumpWriter detected a numerical error writing the"
677 <                   " fluctuating charge force for object %s", id.c_str());      
677 >                   " fluctuating charge velocity for object %s", id.c_str());      
678            painCave.isFatal = 1;
679            simError();
680          }
681 <        sprintf(tempBuffer, " %13e ", fqFrc);        
681 >        sprintf(tempBuffer, " %13e ", fqVel);
682          line += tempBuffer;
683        }
582 –    }
684  
685 +      if (needForceVector_) {
686 +        if (storageLayout & DataStorage::dslFlucQForce) {          
687 +          type += "g";
688 +          RealType fqFrc = sd->getFlucQFrc();        
689 +          if (isinf(fqFrc) || isnan(fqFrc) ) {      
690 +            sprintf( painCave.errMsg,
691 +                     "DumpWriter detected a numerical error writing the"
692 +                     " fluctuating charge force for object %s", id.c_str());      
693 +            painCave.isFatal = 1;
694 +            simError();
695 +          }
696 +          sprintf(tempBuffer, " %13e ", fqFrc);        
697 +          line += tempBuffer;
698 +        }
699 +      }
700 +    }
701 +    
702      if (needElectricField_) {
703 <      type += "e";
704 <      Vector3d eField= sd->getElectricField();
705 <      if (isinf(eField[0]) || isnan(eField[0]) ||
706 <          isinf(eField[1]) || isnan(eField[1]) ||
707 <          isinf(eField[2]) || isnan(eField[2]) ) {      
708 <        sprintf( painCave.errMsg,
709 <                 "DumpWriter detected a numerical error writing the electric"
710 <                 " field for object %s", id.c_str());      
711 <        painCave.isFatal = 1;
712 <        simError();
703 >      if (storageLayout & DataStorage::dslElectricField) {
704 >        type += "e";
705 >        Vector3d eField= sd->getElectricField();
706 >        if (isinf(eField[0]) || isnan(eField[0]) ||
707 >            isinf(eField[1]) || isnan(eField[1]) ||
708 >            isinf(eField[2]) || isnan(eField[2]) ) {      
709 >          sprintf( painCave.errMsg,
710 >                   "DumpWriter detected a numerical error writing the electric"
711 >                   " field for object %s", id.c_str());      
712 >          painCave.isFatal = 1;
713 >          simError();
714 >        }
715 >        sprintf(tempBuffer, " %13e %13e %13e",
716 >                eField[0], eField[1], eField[2]);
717 >        line += tempBuffer;
718        }
596 –      sprintf(tempBuffer, " %13e %13e %13e",
597 –              eField[0], eField[1], eField[2]);
598 –      line += tempBuffer;
719      }
720  
721 <
721 >    if (needSitePotential_) {
722 >      if (storageLayout & DataStorage::dslSitePotential) {          
723 >        type += "s";
724 >        RealType sPot = sd->getSitePotential();        
725 >        if (isinf(sPot) || isnan(sPot) ) {      
726 >          sprintf( painCave.errMsg,
727 >                   "DumpWriter detected a numerical error writing the"
728 >                   " site potential for object %s", id.c_str());      
729 >          painCave.isFatal = 1;
730 >          simError();
731 >        }
732 >        sprintf(tempBuffer, " %13e ", sPot);        
733 >        line += tempBuffer;
734 >      }
735 >    }    
736 >    
737      if (needParticlePot_) {
738 <      type += "u";
739 <      RealType particlePot = sd->getParticlePot();
740 <      if (isinf(particlePot) || isnan(particlePot)) {      
741 <        sprintf( painCave.errMsg,
742 <                 "DumpWriter detected a numerical error writing the particle "
743 <                 " potential for object %s", id.c_str());      
744 <        painCave.isFatal = 1;
745 <        simError();
738 >      if (storageLayout & DataStorage::dslParticlePot) {
739 >        type += "u";
740 >        RealType particlePot = sd->getParticlePot();
741 >        if (isinf(particlePot) || isnan(particlePot)) {      
742 >          sprintf( painCave.errMsg,
743 >                   "DumpWriter detected a numerical error writing the particle "
744 >                   " potential for object %s", id.c_str());      
745 >          painCave.isFatal = 1;
746 >          simError();
747 >        }
748 >        sprintf(tempBuffer, " %13e", particlePot);
749 >        line += tempBuffer;
750        }
612 –      sprintf(tempBuffer, " %13e", particlePot);
613 –      line += tempBuffer;
751      }
752 <    
616 <
752 >  
753      sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
754      return std::string(tempBuffer);
755    }
# Line 623 | Line 759 | namespace OpenMD {
759    }
760  
761    void DumpWriter::writeEor() {
762 <    std::ostream* eorStream;
763 <    
762 >
763 >    std::ostream* eorStream = NULL;
764 >
765   #ifdef IS_MPI
766      if (worldRank == 0) {
767   #endif // is_mpi
768 <
768 >      
769        eorStream = createOStream(eorFilename_);
770  
771   #ifdef IS_MPI
772      }
773 < #endif // is_mpi    
774 <
773 > #endif
774 >    
775      writeFrame(*eorStream);
776 <
776 >      
777   #ifdef IS_MPI
778      if (worldRank == 0) {
779 < #endif // is_mpi
779 > #endif
780 >      
781        writeClosing(*eorStream);
782        delete eorStream;
783 +      
784   #ifdef IS_MPI
785      }
786   #endif // is_mpi  
# Line 655 | Line 794 | namespace OpenMD {
794   #ifdef IS_MPI
795      if (worldRank == 0) {
796   #endif // is_mpi
658 –
797        buffers.push_back(dumpFile_->rdbuf());
660 –
798        eorStream = createOStream(eorFilename_);
662 –
799        buffers.push_back(eorStream->rdbuf());
664 –        
800   #ifdef IS_MPI
801      }
802   #endif // is_mpi    
803  
804      TeeBuf tbuf(buffers.begin(), buffers.end());
805      std::ostream os(&tbuf);
671 –
806      writeFrame(os);
807  
808   #ifdef IS_MPI
# Line 678 | Line 812 | namespace OpenMD {
812        delete eorStream;
813   #ifdef IS_MPI
814      }
815 < #endif // is_mpi  
682 <    
815 > #endif // is_mpi      
816    }
817  
818    std::ostream* DumpWriter::createOStream(const std::string& filename) {

Diff Legend

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