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 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).          
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 +        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
326 <      
327 <        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        }
338 –      os << "    </StuntDoubles>\n";
339 –      
340 –      os << "  </Snapshot>\n";
341 –      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 <  }
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 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
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) {
530          
531 <    int index = integrableObject->getGlobalIntegrableObjectIndex();
531 >    int index = sd->getGlobalIntegrableObjectIndex();
532      std::string type("pv");
533      std::string line;
534      char tempBuffer[4096];
535  
536      Vector3d pos;
537      Vector3d vel;
538 <    pos = integrableObject->getPos();
538 >    pos = sd->getPos();
539  
540      if (isinf(pos[0]) || isnan(pos[0]) ||
541          isinf(pos[1]) || isnan(pos[1]) ||
# Line 376 | Line 547 | namespace OpenMD {
547        simError();
548      }
549  
550 <    vel = integrableObject->getVel();          
550 >    vel = sd->getVel();        
551  
552      if (isinf(vel[0]) || isnan(vel[0]) ||
553          isinf(vel[1]) || isnan(vel[1]) ||
# Line 393 | Line 564 | namespace OpenMD {
564              vel[0], vel[1], vel[2]);                    
565      line += tempBuffer;
566  
567 <    if (integrableObject->isDirectional()) {
567 >    if (sd->isDirectional()) {
568        type += "qj";
569        Quat4d q;
570        Vector3d ji;
571 <      q = integrableObject->getQ();
571 >      q = sd->getQ();
572  
573        if (isinf(q[0]) || isnan(q[0]) ||
574            isinf(q[1]) || isnan(q[1]) ||
# Line 410 | Line 581 | namespace OpenMD {
581          simError();
582        }
583  
584 <      ji = integrableObject->getJ();
584 >      ji = sd->getJ();
585  
586        if (isinf(ji[0]) || isnan(ji[0]) ||
587            isinf(ji[1]) || isnan(ji[1]) ||
# Line 430 | Line 601 | namespace OpenMD {
601  
602      if (needForceVector_) {
603        type += "f";
604 <      Vector3d frc;
434 <
435 <      frc = integrableObject->getFrc();
436 <
604 >      Vector3d frc = sd->getFrc();
605        if (isinf(frc[0]) || isnan(frc[0]) ||
606            isinf(frc[1]) || isnan(frc[1]) ||
607            isinf(frc[2]) || isnan(frc[2]) ) {      
# Line 447 | Line 615 | namespace OpenMD {
615                frc[0], frc[1], frc[2]);
616        line += tempBuffer;
617        
618 <      if (integrableObject->isDirectional()) {
618 >      if (sd->isDirectional()) {
619          type += "t";
620 <        Vector3d trq;
453 <        
454 <        trq = integrableObject->getTrq();
455 <        
620 >        Vector3d trq = sd->getTrq();        
621          if (isinf(trq[0]) || isnan(trq[0]) ||
622              isinf(trq[1]) || isnan(trq[1]) ||
623              isinf(trq[2]) || isnan(trq[2]) ) {      
# Line 461 | Line 626 | namespace OpenMD {
626                     " for object %d", index);      
627            painCave.isFatal = 1;
628            simError();
629 <        }
465 <        
629 >        }        
630          sprintf(tempBuffer, " %13e %13e %13e",
631                  trq[0], trq[1], trq[2]);
632          line += tempBuffer;
633        }      
634      }
471 –    if (needParticlePot_) {
472 –      type += "u";
473 –      RealType particlePot;
635  
636 <      particlePot = integrableObject->getParticlePot();
636 >    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
637 >    return std::string(tempBuffer);
638 >  }
639  
640 <      if (isinf(particlePot) || isnan(particlePot)) {      
641 <        sprintf( painCave.errMsg,
642 <                 "DumpWriter detected a numerical error writing the particle "
643 <                 " potential for object %d", index);      
644 <        painCave.isFatal = 1;
645 <        simError();
646 <      }
647 <      sprintf(tempBuffer, " %13e", particlePot);
648 <      line += tempBuffer;
640 >  std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
641 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
642 >
643 >    std::string id;
644 >    std::string type;
645 >    std::string line;
646 >    char tempBuffer[4096];
647 >
648 >    if (sd->isRigidBody()) {
649 >      sprintf(tempBuffer, "%10d           ", ioIndex);
650 >      id = std::string(tempBuffer);
651 >    } else {
652 >      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
653 >      id = std::string(tempBuffer);
654      }
655 +              
656 +    if (needFlucQ_) {
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 (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 velocity for object %s", id.c_str());      
678 +          painCave.isFatal = 1;
679 +          simError();
680 +        }
681 +        sprintf(tempBuffer, " %13e ", fqVel);
682 +        line += tempBuffer;
683 +      }
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 <    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
702 >    if (needElectricField_) {
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 >      }
719 >    }
720 >
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 >      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 >      }
751 >    }
752 >  
753 >    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
754      return std::string(tempBuffer);
755    }
756  
# Line 494 | 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 526 | Line 794 | namespace OpenMD {
794   #ifdef IS_MPI
795      if (worldRank == 0) {
796   #endif // is_mpi
529 –
797        buffers.push_back(dumpFile_->rdbuf());
531 –
798        eorStream = createOStream(eorFilename_);
533 –
799        buffers.push_back(eorStream->rdbuf());
535 –        
800   #ifdef IS_MPI
801      }
802   #endif // is_mpi    
803  
804      TeeBuf tbuf(buffers.begin(), buffers.end());
805      std::ostream os(&tbuf);
542 –
806      writeFrame(os);
807  
808   #ifdef IS_MPI
# Line 549 | Line 812 | namespace OpenMD {
812        delete eorStream;
813   #ifdef IS_MPI
814      }
815 < #endif // is_mpi  
553 <    
815 > #endif // is_mpi      
816    }
817  
818    std::ostream* DumpWriter::createOStream(const std::string& filename) {
819  
820      std::ostream* newOStream;
821 < #ifdef HAVE_LIBZ
821 > #ifdef HAVE_ZLIB
822      if (needCompression_) {
823        newOStream = new ogzstream(filename.c_str());
824      } else {
# Line 566 | Line 828 | namespace OpenMD {
828      newOStream = new std::ofstream(filename.c_str());
829   #endif
830      //write out MetaData first
831 <    (*newOStream) << "<OpenMD version=1>" << std::endl;
831 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
832      (*newOStream) << "  <MetaData>" << std::endl;
833      (*newOStream) << info_->getRawMetaData();
834      (*newOStream) << "  </MetaData>" << std::endl;

Diff Legend

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