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 1938 by gezelter, Thu Oct 31 15:32:17 2013 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 52 | Line 58
58   #ifdef _MSC_VER
59   #define isnan(x) _isnan((x))
60   #define isinf(x) (!_finite(x) && !_isnan(x))
55 – #endif
56 –
57 – #ifdef IS_MPI
58 – #include <mpi.h>
61   #endif
62  
63   using namespace std;
# Line 302 | Line 304 | namespace OpenMD {
304    void DumpWriter::writeFrame(std::ostream& os) {
305  
306   #ifdef IS_MPI
307 <    MPI_Status istatus;
307 >    MPI::Status istatus;
308   #endif
309  
310      Molecule* mol;
# Line 310 | Line 312 | namespace OpenMD {
312      SimInfo::MoleculeIterator mi;
313      Molecule::IntegrableObjectIterator ii;
314      RigidBody::AtomIterator ai;
313 –    Atom* atom;
315  
316   #ifndef IS_MPI
317      os << "  <Snapshot>\n";
# Line 318 | Line 319 | namespace OpenMD {
319      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
320  
321      os << "    <StuntDoubles>\n";
322 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
323 <
322 >    for (mol = info_->beginMolecule(mi); mol != NULL;
323 >         mol = info_->nextMolecule(mi)) {
324        
325        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
326             sd = mol->nextIntegrableObject(ii)) {        
# Line 331 | Line 332 | namespace OpenMD {
332  
333      if (doSiteData_) {
334        os << "    <SiteData>\n";
335 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
335 >      for (mol = info_->beginMolecule(mi); mol != NULL;
336 >           mol = info_->nextMolecule(mi)) {
337                
338          for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
339             sd = mol->nextIntegrableObject(ii)) {        
# Line 344 | Line 346 | namespace OpenMD {
346              
347              RigidBody* rb = static_cast<RigidBody*>(sd);
348              int siteIndex = 0;
349 <            for (atom = rb->beginAtom(ai); atom != NULL;  
349 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
350                   atom = rb->nextAtom(ai)) {                                            
351                os << prepareSiteLine(atom, ioIndex, siteIndex);
352                siteIndex++;
# Line 358 | Line 360 | namespace OpenMD {
360  
361      os.flush();
362   #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)) {
363  
364 +    const int masterNode = 0;
365 +    int worldRank = MPI::COMM_WORLD.Get_rank();
366 +    int nProc = MPI::COMM_WORLD.Get_size();
367  
368 +    if (worldRank == masterNode) {      
369 +      os << "  <Snapshot>\n";  
370 +      writeFrameProperties(os,
371 +                           info_->getSnapshotManager()->getCurrentSnapshot());
372 +      os << "    <StuntDoubles>\n";
373 +    }
374 +
375 +    //every node prepares the dump lines for integrable objects belong to itself
376 +    std::string buffer;
377 +    for (mol = info_->beginMolecule(mi); mol != NULL;
378 +         mol = info_->nextMolecule(mi)) {
379        for (sd = mol->beginIntegrableObject(ii); sd != NULL;
380             sd = mol->nextIntegrableObject(ii)) {        
381 <          buffer += prepareDumpLine(sd);
381 >        buffer += prepareDumpLine(sd);
382        }
383      }
384      
372 –    const int masterNode = 0;
373 –    int nProc;
374 –    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
385      if (worldRank == masterNode) {      
376 –      os << "  <Snapshot>\n";  
377 –      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
378 –      os << "    <StuntDoubles>\n";
379 –        
386        os << buffer;
387 <
387 >      
388        for (int i = 1; i < nProc; ++i) {
389 +        // tell processor i to start sending us data:
390 +        MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
391  
392          // receive the length of the string buffer that was
393 <        // prepared by processor i
386 <
387 <        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
393 >        // prepared by processor i:        
394          int recvLength;
395 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
395 >        MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
396 >                             istatus);
397 >
398 >        // create a buffer to receive the data
399          char* recvBuffer = new char[recvLength];
400          if (recvBuffer == NULL) {
401          } else {
402 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
402 >          // receive the data:
403 >          MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
404 >                               MPI::ANY_TAG, istatus);
405 >          // send it to the file:
406            os << recvBuffer;
407 +          // get rid of the receive buffer:
408            delete [] recvBuffer;
409          }
410        }
398 –      os << "    </StuntDoubles>\n";
399 –      
400 –      os << "  </Snapshot>\n";
401 –      os.flush();
411      } else {
412        int sendBufferLength = buffer.size() + 1;
413        int myturn = 0;
414        for (int i = 1; i < nProc; ++i){
415 <        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
415 >        // wait for the master node to call our number:
416 >        MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
417          if (myturn == worldRank){
418 <          MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
419 <          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
418 >          // send the length of our buffer:
419 >          MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
420 >
421 >          // send our buffer:
422 >          MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
423 >                               MPI::CHAR, masterNode, 0);
424 >
425          }
426        }
427      }
428 +    
429 +    if (worldRank == masterNode) {      
430 +      os << "    </StuntDoubles>\n";
431 +    }
432  
433 < #endif // is_mpi
433 >    if (doSiteData_) {
434 >      if (worldRank == masterNode) {
435 >        os << "    <SiteData>\n";
436 >      }
437 >      buffer.clear();
438 >      for (mol = info_->beginMolecule(mi); mol != NULL;
439 >           mol = info_->nextMolecule(mi)) {
440 >              
441 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
442 >             sd = mol->nextIntegrableObject(ii)) {      
443 >          
444 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
445 >          // do one for the IO itself
446 >          buffer += prepareSiteLine(sd, ioIndex, 0);
447  
448 +          if (sd->isRigidBody()) {
449 +            
450 +            RigidBody* rb = static_cast<RigidBody*>(sd);
451 +            int siteIndex = 0;
452 +            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
453 +                 atom = rb->nextAtom(ai)) {                                            
454 +              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
455 +              siteIndex++;
456 +            }
457 +          }
458 +        }
459 +      }
460 +
461 +      if (worldRank == masterNode) {    
462 +        os << buffer;
463 +        
464 +        for (int i = 1; i < nProc; ++i) {
465 +          
466 +          // tell processor i to start sending us data:
467 +          MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
468 +          
469 +          // receive the length of the string buffer that was
470 +          // prepared by processor i:        
471 +          int recvLength;
472 +          MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
473 +                               istatus);
474 +          
475 +          // create a buffer to receive the data
476 +          char* recvBuffer = new char[recvLength];
477 +          if (recvBuffer == NULL) {
478 +          } else {
479 +            // receive the data:
480 +            MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
481 +                                 MPI::ANY_TAG, istatus);
482 +            // send it to the file:
483 +            os << recvBuffer;
484 +            // get rid of the receive buffer:
485 +            delete [] recvBuffer;
486 +          }
487 +        }      
488 +      } else {
489 +        int sendBufferLength = buffer.size() + 1;
490 +        int myturn = 0;
491 +        for (int i = 1; i < nProc; ++i){
492 +          // wait for the master node to call our number:
493 +          MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
494 +          if (myturn == worldRank){
495 +            // send the length of our buffer:
496 +            MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
497 +            // send our buffer:
498 +            MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
499 +                                 MPI::CHAR, masterNode, 0);
500 +          }
501 +        }
502 +      }
503 +      
504 +      if (worldRank == masterNode) {    
505 +        os << "    </SiteData>\n";
506 +      }
507 +    }
508 +    
509 +    if (worldRank == masterNode) {
510 +      os << "  </Snapshot>\n";
511 +      os.flush();
512 +    }
513 +    
514 + #endif // is_mpi
515 +    
516    }
517  
518    std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
# Line 527 | Line 627 | namespace OpenMD {
627    }
628  
629    std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
630 <        
630 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
631  
632      std::string id;
633      std::string type;
# Line 543 | Line 643 | namespace OpenMD {
643      }
644                
645      if (needFlucQ_) {
646 <      type += "cw";
647 <      RealType fqPos = sd->getFlucQPos();
648 <      if (isinf(fqPos) || isnan(fqPos) ) {      
649 <        sprintf( painCave.errMsg,
650 <                 "DumpWriter detected a numerical error writing the"
651 <                 " fluctuating charge for object %s", id.c_str());      
652 <        painCave.isFatal = 1;
653 <        simError();
654 <      }
655 <      sprintf(tempBuffer, " %13e ", fqPos);
656 <      line += tempBuffer;
657 <    
658 <      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;
646 >      if (storageLayout & DataStorage::dslFlucQPosition) {
647 >        type += "c";
648 >        RealType fqPos = sd->getFlucQPos();
649 >        if (isinf(fqPos) || isnan(fqPos) ) {      
650 >          sprintf( painCave.errMsg,
651 >                   "DumpWriter detected a numerical error writing the"
652 >                   " fluctuating charge for object %s", id.c_str());      
653 >          painCave.isFatal = 1;
654 >          simError();
655 >        }
656 >        sprintf(tempBuffer, " %13e ", fqPos);
657 >        line += tempBuffer;
658 >      }
659  
660 <      if (needForceVector_) {
661 <        type += "g";
662 <        RealType fqFrc = sd->getFlucQFrc();        
663 <        if (isinf(fqFrc) || isnan(fqFrc) ) {      
660 >      if (storageLayout & DataStorage::dslFlucQVelocity) {
661 >        type += "w";    
662 >        RealType fqVel = sd->getFlucQVel();
663 >        if (isinf(fqVel) || isnan(fqVel) ) {      
664            sprintf( painCave.errMsg,
665                     "DumpWriter detected a numerical error writing the"
666 <                   " fluctuating charge force for object %s", id.c_str());      
666 >                   " fluctuating charge velocity for object %s", id.c_str());      
667            painCave.isFatal = 1;
668            simError();
669          }
670 <        sprintf(tempBuffer, " %13e ", fqFrc);        
670 >        sprintf(tempBuffer, " %13e ", fqVel);
671          line += tempBuffer;
672        }
582 –    }
673  
674 +      if (needForceVector_) {
675 +        if (storageLayout & DataStorage::dslFlucQForce) {          
676 +          type += "g";
677 +          RealType fqFrc = sd->getFlucQFrc();        
678 +          if (isinf(fqFrc) || isnan(fqFrc) ) {      
679 +            sprintf( painCave.errMsg,
680 +                     "DumpWriter detected a numerical error writing the"
681 +                     " fluctuating charge force for object %s", id.c_str());      
682 +            painCave.isFatal = 1;
683 +            simError();
684 +          }
685 +          sprintf(tempBuffer, " %13e ", fqFrc);        
686 +          line += tempBuffer;
687 +        }
688 +      }
689 +    }
690 +    
691      if (needElectricField_) {
692 <      type += "e";
693 <      Vector3d eField= sd->getElectricField();
694 <      if (isinf(eField[0]) || isnan(eField[0]) ||
695 <          isinf(eField[1]) || isnan(eField[1]) ||
696 <          isinf(eField[2]) || isnan(eField[2]) ) {      
697 <        sprintf( painCave.errMsg,
698 <                 "DumpWriter detected a numerical error writing the electric"
699 <                 " field for object %s", id.c_str());      
700 <        painCave.isFatal = 1;
701 <        simError();
692 >      if (storageLayout & DataStorage::dslElectricField) {
693 >        type += "e";
694 >        Vector3d eField= sd->getElectricField();
695 >        if (isinf(eField[0]) || isnan(eField[0]) ||
696 >            isinf(eField[1]) || isnan(eField[1]) ||
697 >            isinf(eField[2]) || isnan(eField[2]) ) {      
698 >          sprintf( painCave.errMsg,
699 >                   "DumpWriter detected a numerical error writing the electric"
700 >                   " field for object %s", id.c_str());      
701 >          painCave.isFatal = 1;
702 >          simError();
703 >        }
704 >        sprintf(tempBuffer, " %13e %13e %13e",
705 >                eField[0], eField[1], eField[2]);
706 >        line += tempBuffer;
707        }
596 –      sprintf(tempBuffer, " %13e %13e %13e",
597 –              eField[0], eField[1], eField[2]);
598 –      line += tempBuffer;
708      }
709  
710  
711      if (needParticlePot_) {
712 <      type += "u";
713 <      RealType particlePot = sd->getParticlePot();
714 <      if (isinf(particlePot) || isnan(particlePot)) {      
715 <        sprintf( painCave.errMsg,
716 <                 "DumpWriter detected a numerical error writing the particle "
717 <                 " potential for object %s", id.c_str());      
718 <        painCave.isFatal = 1;
719 <        simError();
712 >      if (storageLayout & DataStorage::dslParticlePot) {
713 >        type += "u";
714 >        RealType particlePot = sd->getParticlePot();
715 >        if (isinf(particlePot) || isnan(particlePot)) {      
716 >          sprintf( painCave.errMsg,
717 >                   "DumpWriter detected a numerical error writing the particle "
718 >                   " potential for object %s", id.c_str());      
719 >          painCave.isFatal = 1;
720 >          simError();
721 >        }
722 >        sprintf(tempBuffer, " %13e", particlePot);
723 >        line += tempBuffer;
724        }
612 –      sprintf(tempBuffer, " %13e", particlePot);
613 –      line += tempBuffer;
725      }
726 <    
616 <
726 >  
727      sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
728      return std::string(tempBuffer);
729    }
# Line 628 | Line 738 | namespace OpenMD {
738   #ifdef IS_MPI
739      if (worldRank == 0) {
740   #endif // is_mpi
741 <
741 >      
742        eorStream = createOStream(eorFilename_);
743  
744   #ifdef IS_MPI
745      }
746 < #endif // is_mpi    
747 <
746 > #endif
747 >    
748      writeFrame(*eorStream);
749 <
749 >      
750   #ifdef IS_MPI
751      if (worldRank == 0) {
752 < #endif // is_mpi
752 > #endif
753 >      
754        writeClosing(*eorStream);
755        delete eorStream;
756 +      
757   #ifdef IS_MPI
758      }
759   #endif // is_mpi  

Diff Legend

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