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 1610 by gezelter, Fri Aug 12 14:37:25 2011 UTC vs.
Revision 2020 by gezelter, Mon Sep 22 19:18:35 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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
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"
51   #include "utils/simError.h"
52   #include "io/basic_teebuf.hpp"
53 + #ifdef HAVE_ZLIB
54   #include "io/gzstream.hpp"
55 + #endif
56   #include "io/Globals.hpp"
57  
58 + #ifdef _MSC_VER
59 + #define isnan(x) _isnan((x))
60 + #define isinf(x) (!_finite(x) && !_isnan(x))
61 + #endif
62  
50 – #ifdef IS_MPI
51 – #include <mpi.h>
52 – #endif //is_mpi
53 –
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();
72 <    needForceVector_ = simParams->getOutputForceVector();
73 <    needParticlePot_ = simParams->getOutputParticlePotential();
71 >    needCompression_   = simParams->getCompressDumpFile();
72 >    needForceVector_   = simParams->getOutputForceVector();
73 >    needParticlePot_   = simParams->getOutputParticlePotential();
74 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
75 >    needElectricField_ = simParams->getOutputElectricField();
76 >    needSitePotential_ = simParams->getOutputSitePotential();
77 >
78 >    if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
79 >        needSitePotential_) {
80 >      doSiteData_ = true;
81 >    } else {
82 >      doSiteData_ = false;
83 >    }
84 >
85      createDumpFile_ = true;
86   #ifdef HAVE_LIBZ
87      if (needCompression_) {
# Line 98 | Line 119 | namespace OpenMD {
119      Globals* simParams = info->getSimParams();
120      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
121  
122 <    needCompression_ = simParams->getCompressDumpFile();
123 <    needForceVector_ = simParams->getOutputForceVector();
124 <    needParticlePot_ = simParams->getOutputParticlePotential();
122 >    needCompression_   = simParams->getCompressDumpFile();
123 >    needForceVector_   = simParams->getOutputForceVector();
124 >    needParticlePot_   = simParams->getOutputParticlePotential();
125 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
126 >    needElectricField_ = simParams->getOutputElectricField();
127 >    needSitePotential_ = simParams->getOutputSitePotential();
128 >
129 >    if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
130 >        needSitePotential_) {
131 >      doSiteData_ = true;
132 >    } else {
133 >      doSiteData_ = false;
134 >    }
135 >
136      createDumpFile_ = true;
137   #ifdef HAVE_LIBZ
138      if (needCompression_) {
# Line 138 | Line 170 | namespace OpenMD {
170      Globals* simParams = info->getSimParams();
171      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
172      
173 <    needCompression_ = simParams->getCompressDumpFile();
174 <    needForceVector_ = simParams->getOutputForceVector();
175 <    needParticlePot_ = simParams->getOutputParticlePotential();
176 <    
173 >    needCompression_   = simParams->getCompressDumpFile();
174 >    needForceVector_   = simParams->getOutputForceVector();
175 >    needParticlePot_   = simParams->getOutputParticlePotential();
176 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
177 >    needElectricField_ = simParams->getOutputElectricField();
178 >    needSitePotential_ = simParams->getOutputSitePotential();
179 >
180 >    if (needParticlePot_ || needFlucQ_ || needElectricField_ ||
181 >        needSitePotential_) {
182 >      doSiteData_ = true;
183 >    } else {
184 >      doSiteData_ = false;
185 >    }
186 >
187   #ifdef HAVE_LIBZ
188      if (needCompression_) {
189        filename_ += ".gz";
# Line 230 | Line 272 | namespace OpenMD {
272              hmat(0, 2), hmat(1, 2), hmat(2, 2));
273      os << buffer;
274  
275 <    RealType chi = s->getChi();
276 <    RealType integralOfChiDt = s->getIntegralOfChiDt();
277 <    if (isinf(chi) || isnan(chi) ||
278 <        isinf(integralOfChiDt) || isnan(integralOfChiDt)) {      
275 >    pair<RealType, RealType> thermostat = s->getThermostat();
276 >
277 >    if (isinf(thermostat.first)  || isnan(thermostat.first) ||
278 >        isinf(thermostat.second) || isnan(thermostat.second)) {      
279        sprintf( painCave.errMsg,
280                 "DumpWriter detected a numerical error writing the thermostat");
281        painCave.isFatal = 1;
282        simError();
283      }
284 <    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
284 >    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", thermostat.first,
285 >            thermostat.second);
286      os << buffer;
287  
288      Mat3x3d eta;
289 <    eta = s->getEta();
289 >    eta = s->getBarostat();
290  
291      for (unsigned int i = 0; i < 3; i++) {
292        for (unsigned int j = 0; j < 3; j++) {
# Line 272 | Line 315 | namespace OpenMD {
315   #endif
316  
317      Molecule* mol;
318 <    StuntDouble* integrableObject;
318 >    StuntDouble* sd;
319      SimInfo::MoleculeIterator mi;
320      Molecule::IntegrableObjectIterator ii;
321 +    RigidBody::AtomIterator ai;
322  
323   #ifndef IS_MPI
324      os << "  <Snapshot>\n";
# Line 282 | 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 (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
333 <           integrableObject = mol->nextIntegrableObject(ii)) {  
334 <          os << prepareDumpLine(integrableObject);
332 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
333 >           sd = mol->nextIntegrableObject(ii)) {        
334 >          os << prepareDumpLine(sd);
335            
336        }
337      }    
338      os << "    </StuntDoubles>\n";
339 <    
339 >
340 >    if (doSiteData_) {
341 >      os << "    <SiteData>\n";
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)) {        
347 >
348 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
349 >          // do one for the IO itself
350 >          os << prepareSiteLine(sd, ioIndex, 0);
351 >
352 >          if (sd->isRigidBody()) {
353 >            
354 >            RigidBody* rb = static_cast<RigidBody*>(sd);
355 >            int siteIndex = 0;
356 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
357 >                 atom = rb->nextAtom(ai)) {                                            
358 >              os << prepareSiteLine(atom, ioIndex, siteIndex);
359 >              siteIndex++;
360 >            }
361 >          }
362 >        }
363 >      }    
364 >      os << "    </SiteData>\n";
365 >    }
366      os << "  </Snapshot>\n";
367  
368      os.flush();
369   #else
300 –    //every node prepares the dump lines for integrable objects belong to itself
301 –    std::string buffer;
302 –    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
370  
304 –
305 –      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
306 –           integrableObject = mol->nextIntegrableObject(ii)) {  
307 –          buffer += prepareDumpLine(integrableObject);
308 –      }
309 –    }
310 –    
371      const int masterNode = 0;
372 +    int worldRank;
373      int nProc;
374 <    MPI_Comm_size(MPI_COMM_WORLD, &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, info_->getSnapshotManager()->getCurrentSnapshot());
381 >      writeFrameProperties(os,
382 >                           info_->getSnapshotManager()->getCurrentSnapshot());
383        os << "    <StuntDoubles>\n";
384 <        
319 <      os << buffer;
384 >    }
385  
386 <    
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);
393 >      }
394 >    }
395 >    
396 >    if (worldRank == masterNode) {      
397 >      os << buffer;
398 >      
399        for (int i = 1; i < nProc; ++i) {
400 +        // tell processor i to start sending us data:
401  
402 +        MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
403 +
404          // receive the length of the string buffer that was
405 <        // prepared by processor i
326 <      
327 <        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
405 >        // prepared by processor i:        
406          int recvLength;
407 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
407 >        MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
408 >                 &istatus);
409 >
410 >        // create a buffer to receive the data
411          char* recvBuffer = new char[recvLength];
412          if (recvBuffer == NULL) {
413          } else {
414 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
414 >          // receive the data:
415 >          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
416 >                               MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
417 >          // send it to the file:
418            os << recvBuffer;
419 +          // get rid of the receive buffer:
420            delete [] recvBuffer;
421          }
422        }
338 –      os << "    </StuntDoubles>\n";
339 –      
340 –      os << "  </Snapshot>\n";
341 –      os.flush();
423      } else {
424        int sendBufferLength = buffer.size() + 1;
425        int myturn = 0;
426        for (int i = 1; i < nProc; ++i){
427 <        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
427 >        // wait for the master node to call our number:
428 >        MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
429          if (myturn == worldRank){
430 +          // send the length of our buffer:
431 +
432            MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
433 <          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
433 >
434 >          // send our buffer:
435 >          MPI_Send((void *)buffer.c_str(), sendBufferLength,
436 >                   MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
437 >
438          }
439        }
440      }
441 +    
442 +    if (worldRank == masterNode) {      
443 +      os << "    </StuntDoubles>\n";
444 +    }
445  
446 < #endif // is_mpi
446 >    if (doSiteData_) {
447 >      if (worldRank == masterNode) {
448 >        os << "    <SiteData>\n";
449 >      }
450 >      buffer.clear();
451 >      for (mol = info_->beginMolecule(mi); mol != NULL;
452 >           mol = info_->nextMolecule(mi)) {
453 >              
454 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
455 >             sd = mol->nextIntegrableObject(ii)) {      
456 >          
457 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
458 >          // do one for the IO itself
459 >          buffer += prepareSiteLine(sd, ioIndex, 0);
460  
461 <  }
461 >          if (sd->isRigidBody()) {
462 >            
463 >            RigidBody* rb = static_cast<RigidBody*>(sd);
464 >            int siteIndex = 0;
465 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
466 >                 atom = rb->nextAtom(ai)) {                                            
467 >              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
468 >              siteIndex++;
469 >            }
470 >          }
471 >        }
472 >      }
473  
474 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
474 >      if (worldRank == masterNode) {    
475 >        os << buffer;
476 >        
477 >        for (int i = 1; i < nProc; ++i) {
478 >          
479 >          // tell processor i to start sending us data:
480 >          MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
481 >          
482 >          // receive the length of the string buffer that was
483 >          // prepared by processor i:        
484 >          int recvLength;
485 >          MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
486 >                   &istatus);
487 >          
488 >          // create a buffer to receive the data
489 >          char* recvBuffer = new char[recvLength];
490 >          if (recvBuffer == NULL) {
491 >          } else {
492 >            // receive the data:
493 >            MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
494 >                     MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
495 >            // send it to the file:
496 >            os << recvBuffer;
497 >            // get rid of the receive buffer:
498 >            delete [] recvBuffer;
499 >          }
500 >        }      
501 >      } else {
502 >        int sendBufferLength = buffer.size() + 1;
503 >        int myturn = 0;
504 >        for (int i = 1; i < nProc; ++i){
505 >          // wait for the master node to call our number:
506 >          MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
507 >          if (myturn == worldRank){
508 >            // send the length of our buffer:
509 >            MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
510 >            // send our buffer:
511 >            MPI_Send((void *)buffer.c_str(), sendBufferLength,
512 >                     MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
513 >          }
514 >        }
515 >      }
516 >      
517 >      if (worldRank == masterNode) {    
518 >        os << "    </SiteData>\n";
519 >      }
520 >    }
521 >    
522 >    if (worldRank == masterNode) {
523 >      os << "  </Snapshot>\n";
524 >      os.flush();
525 >    }
526 >    
527 > #endif // is_mpi
528 >    
529 >  }
530 >
531 >  std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
532          
533 <    int index = integrableObject->getGlobalIntegrableObjectIndex();
533 >    int index = sd->getGlobalIntegrableObjectIndex();
534      std::string type("pv");
535      std::string line;
536      char tempBuffer[4096];
537  
538      Vector3d pos;
539      Vector3d vel;
540 <    pos = integrableObject->getPos();
540 >    pos = sd->getPos();
541  
542      if (isinf(pos[0]) || isnan(pos[0]) ||
543          isinf(pos[1]) || isnan(pos[1]) ||
# Line 376 | Line 549 | namespace OpenMD {
549        simError();
550      }
551  
552 <    vel = integrableObject->getVel();          
552 >    vel = sd->getVel();        
553  
554      if (isinf(vel[0]) || isnan(vel[0]) ||
555          isinf(vel[1]) || isnan(vel[1]) ||
# Line 393 | Line 566 | namespace OpenMD {
566              vel[0], vel[1], vel[2]);                    
567      line += tempBuffer;
568  
569 <    if (integrableObject->isDirectional()) {
569 >    if (sd->isDirectional()) {
570        type += "qj";
571        Quat4d q;
572        Vector3d ji;
573 <      q = integrableObject->getQ();
573 >      q = sd->getQ();
574  
575        if (isinf(q[0]) || isnan(q[0]) ||
576            isinf(q[1]) || isnan(q[1]) ||
# Line 410 | Line 583 | namespace OpenMD {
583          simError();
584        }
585  
586 <      ji = integrableObject->getJ();
586 >      ji = sd->getJ();
587  
588        if (isinf(ji[0]) || isnan(ji[0]) ||
589            isinf(ji[1]) || isnan(ji[1]) ||
# Line 430 | Line 603 | namespace OpenMD {
603  
604      if (needForceVector_) {
605        type += "f";
606 <      Vector3d frc;
434 <
435 <      frc = integrableObject->getFrc();
436 <
606 >      Vector3d frc = sd->getFrc();
607        if (isinf(frc[0]) || isnan(frc[0]) ||
608            isinf(frc[1]) || isnan(frc[1]) ||
609            isinf(frc[2]) || isnan(frc[2]) ) {      
# Line 447 | Line 617 | namespace OpenMD {
617                frc[0], frc[1], frc[2]);
618        line += tempBuffer;
619        
620 <      if (integrableObject->isDirectional()) {
620 >      if (sd->isDirectional()) {
621          type += "t";
622 <        Vector3d trq;
453 <        
454 <        trq = integrableObject->getTrq();
455 <        
622 >        Vector3d trq = sd->getTrq();        
623          if (isinf(trq[0]) || isnan(trq[0]) ||
624              isinf(trq[1]) || isnan(trq[1]) ||
625              isinf(trq[2]) || isnan(trq[2]) ) {      
# Line 461 | Line 628 | namespace OpenMD {
628                     " for object %d", index);      
629            painCave.isFatal = 1;
630            simError();
631 <        }
465 <        
631 >        }        
632          sprintf(tempBuffer, " %13e %13e %13e",
633                  trq[0], trq[1], trq[2]);
634          line += tempBuffer;
635        }      
636      }
471 –    if (needParticlePot_) {
472 –      type += "u";
473 –      RealType particlePot;
637  
638 <      particlePot = integrableObject->getParticlePot();
638 >    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
639 >    return std::string(tempBuffer);
640 >  }
641  
642 <      if (isinf(particlePot) || isnan(particlePot)) {      
643 <        sprintf( painCave.errMsg,
644 <                 "DumpWriter detected a numerical error writing the particle "
645 <                 " potential for object %d", index);      
646 <        painCave.isFatal = 1;
647 <        simError();
648 <      }
649 <      sprintf(tempBuffer, " %13e", particlePot);
650 <      line += tempBuffer;
642 >  std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
643 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
644 >
645 >    std::string id;
646 >    std::string type;
647 >    std::string line;
648 >    char tempBuffer[4096];
649 >
650 >    if (sd->isRigidBody()) {
651 >      sprintf(tempBuffer, "%10d           ", ioIndex);
652 >      id = std::string(tempBuffer);
653 >    } else {
654 >      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
655 >      id = std::string(tempBuffer);
656 >    }
657 >              
658 >    if (needFlucQ_) {
659 >      if (storageLayout & DataStorage::dslFlucQPosition) {
660 >        type += "c";
661 >        RealType fqPos = sd->getFlucQPos();
662 >        if (isinf(fqPos) || isnan(fqPos) ) {      
663 >          sprintf( painCave.errMsg,
664 >                   "DumpWriter detected a numerical error writing the"
665 >                   " fluctuating charge for object %s", id.c_str());      
666 >          painCave.isFatal = 1;
667 >          simError();
668 >        }
669 >        sprintf(tempBuffer, " %13e ", fqPos);
670 >        line += tempBuffer;
671 >      }
672 >
673 >      if (storageLayout & DataStorage::dslFlucQVelocity) {
674 >        type += "w";    
675 >        RealType fqVel = sd->getFlucQVel();
676 >        if (isinf(fqVel) || isnan(fqVel) ) {      
677 >          sprintf( painCave.errMsg,
678 >                   "DumpWriter detected a numerical error writing the"
679 >                   " fluctuating charge velocity for object %s", id.c_str());      
680 >          painCave.isFatal = 1;
681 >          simError();
682 >        }
683 >        sprintf(tempBuffer, " %13e ", fqVel);
684 >        line += tempBuffer;
685 >      }
686 >
687 >      if (needForceVector_) {
688 >        if (storageLayout & DataStorage::dslFlucQForce) {          
689 >          type += "g";
690 >          RealType fqFrc = sd->getFlucQFrc();        
691 >          if (isinf(fqFrc) || isnan(fqFrc) ) {      
692 >            sprintf( painCave.errMsg,
693 >                     "DumpWriter detected a numerical error writing the"
694 >                     " fluctuating charge force for object %s", id.c_str());      
695 >            painCave.isFatal = 1;
696 >            simError();
697 >          }
698 >          sprintf(tempBuffer, " %13e ", fqFrc);        
699 >          line += tempBuffer;
700 >        }
701 >      }
702      }
703      
704 <    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
704 >    if (needElectricField_) {
705 >      if (storageLayout & DataStorage::dslElectricField) {
706 >        type += "e";
707 >        Vector3d eField= sd->getElectricField();
708 >        if (isinf(eField[0]) || isnan(eField[0]) ||
709 >            isinf(eField[1]) || isnan(eField[1]) ||
710 >            isinf(eField[2]) || isnan(eField[2]) ) {      
711 >          sprintf( painCave.errMsg,
712 >                   "DumpWriter detected a numerical error writing the electric"
713 >                   " field for object %s", id.c_str());      
714 >          painCave.isFatal = 1;
715 >          simError();
716 >        }
717 >        sprintf(tempBuffer, " %13e %13e %13e",
718 >                eField[0], eField[1], eField[2]);
719 >        line += tempBuffer;
720 >      }
721 >    }
722 >
723 >    if (needSitePotential_) {
724 >      if (storageLayout & DataStorage::dslSitePotential) {          
725 >        type += "s";
726 >        RealType sPot = sd->getSitePotential();        
727 >        if (isinf(sPot) || isnan(sPot) ) {      
728 >          sprintf( painCave.errMsg,
729 >                   "DumpWriter detected a numerical error writing the"
730 >                   " site potential for object %s", id.c_str());      
731 >          painCave.isFatal = 1;
732 >          simError();
733 >        }
734 >        sprintf(tempBuffer, " %13e ", sPot);        
735 >        line += tempBuffer;
736 >      }
737 >    }    
738 >    
739 >    if (needParticlePot_) {
740 >      if (storageLayout & DataStorage::dslParticlePot) {
741 >        type += "u";
742 >        RealType particlePot = sd->getParticlePot();
743 >        if (isinf(particlePot) || isnan(particlePot)) {      
744 >          sprintf( painCave.errMsg,
745 >                   "DumpWriter detected a numerical error writing the particle "
746 >                   " potential for object %s", id.c_str());      
747 >          painCave.isFatal = 1;
748 >          simError();
749 >        }
750 >        sprintf(tempBuffer, " %13e", particlePot);
751 >        line += tempBuffer;
752 >      }
753 >    }
754 >  
755 >    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
756      return std::string(tempBuffer);
757    }
758  
# Line 494 | Line 761 | namespace OpenMD {
761    }
762  
763    void DumpWriter::writeEor() {
764 <    std::ostream* eorStream;
765 <    
764 >
765 >    std::ostream* eorStream = NULL;
766 >
767   #ifdef IS_MPI
768      if (worldRank == 0) {
769   #endif // is_mpi
770 <
770 >      
771        eorStream = createOStream(eorFilename_);
772  
773   #ifdef IS_MPI
774      }
775 < #endif // is_mpi    
776 <
775 > #endif
776 >    
777      writeFrame(*eorStream);
778 <
778 >      
779   #ifdef IS_MPI
780      if (worldRank == 0) {
781 < #endif // is_mpi
781 > #endif
782 >      
783        writeClosing(*eorStream);
784        delete eorStream;
785 +      
786   #ifdef IS_MPI
787      }
788   #endif // is_mpi  
# Line 526 | Line 796 | namespace OpenMD {
796   #ifdef IS_MPI
797      if (worldRank == 0) {
798   #endif // is_mpi
529 –
799        buffers.push_back(dumpFile_->rdbuf());
531 –
800        eorStream = createOStream(eorFilename_);
533 –
801        buffers.push_back(eorStream->rdbuf());
535 –        
802   #ifdef IS_MPI
803      }
804   #endif // is_mpi    
805  
806      TeeBuf tbuf(buffers.begin(), buffers.end());
807      std::ostream os(&tbuf);
542 –
808      writeFrame(os);
809  
810   #ifdef IS_MPI
# Line 549 | Line 814 | namespace OpenMD {
814        delete eorStream;
815   #ifdef IS_MPI
816      }
817 < #endif // is_mpi  
553 <    
817 > #endif // is_mpi      
818    }
819  
820    std::ostream* DumpWriter::createOStream(const std::string& filename) {
821  
822      std::ostream* newOStream;
823 < #ifdef HAVE_LIBZ
823 > #ifdef HAVE_ZLIB
824      if (needCompression_) {
825        newOStream = new ogzstream(filename.c_str());
826      } else {
# Line 566 | Line 830 | namespace OpenMD {
830      newOStream = new std::ofstream(filename.c_str());
831   #endif
832      //write out MetaData first
833 <    (*newOStream) << "<OpenMD version=1>" << std::endl;
833 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
834      (*newOStream) << "  <MetaData>" << std::endl;
835      (*newOStream) << info_->getRawMetaData();
836      (*newOStream) << "  </MetaData>" << std::endl;

Diff Legend

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