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 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1983 by gezelter, Tue Apr 15 20:36:19 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 IS_MPI
59 < #include <mpi.h>
60 < #endif //is_mpi
58 > #ifdef _MSC_VER
59 > #define isnan(x) _isnan((x))
60 > #define isinf(x) (!_finite(x) && !_isnan(x))
61 > #endif
62  
63 + using namespace std;
64   namespace OpenMD {
65  
66    DumpWriter::DumpWriter(SimInfo* info)
67      : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
68  
69      Globals* simParams = info->getSimParams();
70 <    needCompression_ = simParams->getCompressDumpFile();
71 <    needForceVector_ = simParams->getOutputForceVector();
70 >    needCompression_   = simParams->getCompressDumpFile();
71 >    needForceVector_   = simParams->getOutputForceVector();
72 >    needParticlePot_   = simParams->getOutputParticlePotential();
73 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
74 >    needElectricField_ = simParams->getOutputElectricField();
75 >
76 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
77 >      doSiteData_ = true;
78 >    } else {
79 >      doSiteData_ = false;
80 >    }
81 >
82      createDumpFile_ = true;
83   #ifdef HAVE_LIBZ
84      if (needCompression_) {
# Line 95 | Line 116 | namespace OpenMD {
116      Globals* simParams = info->getSimParams();
117      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
118  
119 <    needCompression_ = simParams->getCompressDumpFile();
120 <    needForceVector_ = simParams->getOutputForceVector();
119 >    needCompression_   = simParams->getCompressDumpFile();
120 >    needForceVector_   = simParams->getOutputForceVector();
121 >    needParticlePot_   = simParams->getOutputParticlePotential();
122 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
123 >    needElectricField_ = simParams->getOutputElectricField();
124 >
125 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
126 >      doSiteData_ = true;
127 >    } else {
128 >      doSiteData_ = false;
129 >    }
130 >
131      createDumpFile_ = true;
132   #ifdef HAVE_LIBZ
133      if (needCompression_) {
# Line 134 | Line 165 | namespace OpenMD {
165      Globals* simParams = info->getSimParams();
166      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
167      
168 <    needCompression_ = simParams->getCompressDumpFile();
169 <    needForceVector_ = simParams->getOutputForceVector();
170 <    
168 >    needCompression_   = simParams->getCompressDumpFile();
169 >    needForceVector_   = simParams->getOutputForceVector();
170 >    needParticlePot_   = simParams->getOutputParticlePotential();
171 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
172 >    needElectricField_ = simParams->getOutputElectricField();
173 >
174 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
175 >      doSiteData_ = true;
176 >    } else {
177 >      doSiteData_ = false;
178 >    }
179 >
180   #ifdef HAVE_LIBZ
181      if (needCompression_) {
182        filename_ += ".gz";
# Line 194 | Line 234 | namespace OpenMD {
234      os << "    <FrameData>\n";
235  
236      RealType currentTime = s->getTime();
237 +
238 +    if (isinf(currentTime) || isnan(currentTime)) {      
239 +      sprintf( painCave.errMsg,
240 +               "DumpWriter detected a numerical error writing the time");      
241 +      painCave.isFatal = 1;
242 +      simError();
243 +    }
244 +    
245      sprintf(buffer, "        Time: %.10g\n", currentTime);
246      os << buffer;
247  
248      Mat3x3d hmat;
249      hmat = s->getHmat();
250 +
251 +    for (unsigned int i = 0; i < 3; i++) {
252 +      for (unsigned int j = 0; j < 3; j++) {
253 +        if (isinf(hmat(i,j)) || isnan(hmat(i,j))) {      
254 +          sprintf( painCave.errMsg,
255 +                   "DumpWriter detected a numerical error writing the box");
256 +          painCave.isFatal = 1;
257 +          simError();
258 +        }        
259 +      }
260 +    }
261 +    
262      sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
263              hmat(0, 0), hmat(1, 0), hmat(2, 0),
264              hmat(0, 1), hmat(1, 1), hmat(2, 1),
265              hmat(0, 2), hmat(1, 2), hmat(2, 2));
266      os << buffer;
267  
268 <    RealType chi = s->getChi();
269 <    RealType integralOfChiDt = s->getIntegralOfChiDt();
270 <    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
268 >    pair<RealType, RealType> thermostat = s->getThermostat();
269 >
270 >    if (isinf(thermostat.first)  || isnan(thermostat.first) ||
271 >        isinf(thermostat.second) || isnan(thermostat.second)) {      
272 >      sprintf( painCave.errMsg,
273 >               "DumpWriter detected a numerical error writing the thermostat");
274 >      painCave.isFatal = 1;
275 >      simError();
276 >    }
277 >    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", thermostat.first,
278 >            thermostat.second);
279      os << buffer;
280  
281      Mat3x3d eta;
282 <    eta = s->getEta();
282 >    eta = s->getBarostat();
283 >
284 >    for (unsigned int i = 0; i < 3; i++) {
285 >      for (unsigned int j = 0; j < 3; j++) {
286 >        if (isinf(eta(i,j)) || isnan(eta(i,j))) {      
287 >          sprintf( painCave.errMsg,
288 >                   "DumpWriter detected a numerical error writing the barostat");
289 >          painCave.isFatal = 1;
290 >          simError();
291 >        }        
292 >      }
293 >    }
294 >
295      sprintf(buffer, "    Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
296              eta(0, 0), eta(1, 0), eta(2, 0),
297              eta(0, 1), eta(1, 1), eta(2, 1),
# Line 228 | Line 308 | namespace OpenMD {
308   #endif
309  
310      Molecule* mol;
311 <    StuntDouble* integrableObject;
311 >    StuntDouble* sd;
312      SimInfo::MoleculeIterator mi;
313      Molecule::IntegrableObjectIterator ii;
314 +    RigidBody::AtomIterator ai;
315  
316   #ifndef IS_MPI
317      os << "  <Snapshot>\n";
# Line 238 | 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 (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
326 <           integrableObject = mol->nextIntegrableObject(ii)) {  
327 <          os << prepareDumpLine(integrableObject);
325 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
326 >           sd = mol->nextIntegrableObject(ii)) {        
327 >          os << prepareDumpLine(sd);
328            
329        }
330      }    
331      os << "    </StuntDoubles>\n";
332 <    
332 >
333 >    if (doSiteData_) {
334 >      os << "    <SiteData>\n";
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)) {        
340 >
341 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
342 >          // do one for the IO itself
343 >          os << prepareSiteLine(sd, ioIndex, 0);
344 >
345 >          if (sd->isRigidBody()) {
346 >            
347 >            RigidBody* rb = static_cast<RigidBody*>(sd);
348 >            int siteIndex = 0;
349 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
350 >                 atom = rb->nextAtom(ai)) {                                            
351 >              os << prepareSiteLine(atom, ioIndex, siteIndex);
352 >              siteIndex++;
353 >            }
354 >          }
355 >        }
356 >      }    
357 >      os << "    </SiteData>\n";
358 >    }
359      os << "  </Snapshot>\n";
360  
361      os.flush();
362   #else
256    //every node prepares the dump lines for integrable objects belong to itself
257    std::string buffer;
258    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
363  
260
261      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
262           integrableObject = mol->nextIntegrableObject(ii)) {  
263          buffer += prepareDumpLine(integrableObject);
264      }
265    }
266    
364      const int masterNode = 0;
365 +    int worldRank;
366 +    int nProc;
367  
368 +    MPI_Comm_size( MPI_COMM_WORLD, &nProc);
369 +    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
370 +
371 +
372      if (worldRank == masterNode) {      
373        os << "  <Snapshot>\n";  
374 <      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
374 >      writeFrameProperties(os,
375 >                           info_->getSnapshotManager()->getCurrentSnapshot());
376        os << "    <StuntDoubles>\n";
377 <        
274 <      os << buffer;
377 >    }
378  
379 <      int nProc;
380 <      MPI_Comm_size(MPI_COMM_WORLD, &nProc);
379 >    //every node prepares the dump lines for integrable objects belong to itself
380 >    std::string buffer;
381 >    for (mol = info_->beginMolecule(mi); mol != NULL;
382 >         mol = info_->nextMolecule(mi)) {
383 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;
384 >           sd = mol->nextIntegrableObject(ii)) {        
385 >        buffer += prepareDumpLine(sd);
386 >      }
387 >    }
388 >    
389 >    if (worldRank == masterNode) {      
390 >      os << buffer;
391 >      
392        for (int i = 1; i < nProc; ++i) {
393 +        // tell processor i to start sending us data:
394 +        MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
395  
396          // receive the length of the string buffer that was
397 <        // prepared by processor i
282 <
397 >        // prepared by processor i:        
398          int recvLength;
399 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
399 >        MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
400 >                 &istatus);
401 >
402 >        // create a buffer to receive the data
403          char* recvBuffer = new char[recvLength];
404          if (recvBuffer == NULL) {
405          } else {
406 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
406 >          // receive the data:
407 >          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
408 >                               MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
409 >          // send it to the file:
410            os << recvBuffer;
411 +          // get rid of the receive buffer:
412            delete [] recvBuffer;
413          }
414        }
293      os << "    </StuntDoubles>\n";
294      
295      os << "  </Snapshot>\n";
296      os.flush();
415      } else {
416        int sendBufferLength = buffer.size() + 1;
417 <      MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
418 <      MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
417 >      int myturn = 0;
418 >      for (int i = 1; i < nProc; ++i){
419 >        // wait for the master node to call our number:
420 >        MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
421 >        if (myturn == worldRank){
422 >          // send the length of our buffer:
423 >          MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
424 >
425 >          // send our buffer:
426 >          MPI_Send((void *)buffer.c_str(), sendBufferLength,
427 >                   MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
428 >
429 >        }
430 >      }
431      }
432 +    
433 +    if (worldRank == masterNode) {      
434 +      os << "    </StuntDoubles>\n";
435 +    }
436  
437 < #endif // is_mpi
437 >    if (doSiteData_) {
438 >      if (worldRank == masterNode) {
439 >        os << "    <SiteData>\n";
440 >      }
441 >      buffer.clear();
442 >      for (mol = info_->beginMolecule(mi); mol != NULL;
443 >           mol = info_->nextMolecule(mi)) {
444 >              
445 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
446 >             sd = mol->nextIntegrableObject(ii)) {      
447 >          
448 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
449 >          // do one for the IO itself
450 >          buffer += prepareSiteLine(sd, ioIndex, 0);
451  
452 <  }
452 >          if (sd->isRigidBody()) {
453 >            
454 >            RigidBody* rb = static_cast<RigidBody*>(sd);
455 >            int siteIndex = 0;
456 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
457 >                 atom = rb->nextAtom(ai)) {                                            
458 >              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
459 >              siteIndex++;
460 >            }
461 >          }
462 >        }
463 >      }
464  
465 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
465 >      if (worldRank == masterNode) {    
466 >        os << buffer;
467 >        
468 >        for (int i = 1; i < nProc; ++i) {
469 >          
470 >          // tell processor i to start sending us data:
471 >          MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
472 >          
473 >          // receive the length of the string buffer that was
474 >          // prepared by processor i:        
475 >          int recvLength;
476 >          MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
477 >                   &istatus);
478 >          
479 >          // create a buffer to receive the data
480 >          char* recvBuffer = new char[recvLength];
481 >          if (recvBuffer == NULL) {
482 >          } else {
483 >            // receive the data:
484 >            MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
485 >                     MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
486 >            // send it to the file:
487 >            os << recvBuffer;
488 >            // get rid of the receive buffer:
489 >            delete [] recvBuffer;
490 >          }
491 >        }      
492 >      } else {
493 >        int sendBufferLength = buffer.size() + 1;
494 >        int myturn = 0;
495 >        for (int i = 1; i < nProc; ++i){
496 >          // wait for the master node to call our number:
497 >          MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
498 >          if (myturn == worldRank){
499 >            // send the length of our buffer:
500 >            MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
501 >            // send our buffer:
502 >            MPI_Send((void *)buffer.c_str(), sendBufferLength,
503 >                     MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
504 >          }
505 >        }
506 >      }
507 >      
508 >      if (worldRank == masterNode) {    
509 >        os << "    </SiteData>\n";
510 >      }
511 >    }
512 >    
513 >    if (worldRank == masterNode) {
514 >      os << "  </Snapshot>\n";
515 >      os.flush();
516 >    }
517 >    
518 > #endif // is_mpi
519 >    
520 >  }
521 >
522 >  std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
523          
524 <    int index = integrableObject->getGlobalIntegrableObjectIndex();
524 >    int index = sd->getGlobalIntegrableObjectIndex();
525      std::string type("pv");
526      std::string line;
527      char tempBuffer[4096];
528  
529      Vector3d pos;
530      Vector3d vel;
531 <    pos = integrableObject->getPos();
532 <    vel = integrableObject->getVel();          
531 >    pos = sd->getPos();
532 >
533 >    if (isinf(pos[0]) || isnan(pos[0]) ||
534 >        isinf(pos[1]) || isnan(pos[1]) ||
535 >        isinf(pos[2]) || isnan(pos[2]) ) {      
536 >      sprintf( painCave.errMsg,
537 >               "DumpWriter detected a numerical error writing the position"
538 >               " for object %d", index);      
539 >      painCave.isFatal = 1;
540 >      simError();
541 >    }
542 >
543 >    vel = sd->getVel();        
544 >
545 >    if (isinf(vel[0]) || isnan(vel[0]) ||
546 >        isinf(vel[1]) || isnan(vel[1]) ||
547 >        isinf(vel[2]) || isnan(vel[2]) ) {      
548 >      sprintf( painCave.errMsg,
549 >               "DumpWriter detected a numerical error writing the velocity"
550 >               " for object %d", index);      
551 >      painCave.isFatal = 1;
552 >      simError();
553 >    }
554 >
555      sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
556              pos[0], pos[1], pos[2],
557              vel[0], vel[1], vel[2]);                    
558      line += tempBuffer;
559  
560 <    if (integrableObject->isDirectional()) {
560 >    if (sd->isDirectional()) {
561        type += "qj";
562        Quat4d q;
563        Vector3d ji;
564 <      q = integrableObject->getQ();
565 <      ji = integrableObject->getJ();
564 >      q = sd->getQ();
565 >
566 >      if (isinf(q[0]) || isnan(q[0]) ||
567 >          isinf(q[1]) || isnan(q[1]) ||
568 >          isinf(q[2]) || isnan(q[2]) ||
569 >          isinf(q[3]) || isnan(q[3]) ) {      
570 >        sprintf( painCave.errMsg,
571 >                 "DumpWriter detected a numerical error writing the quaternion"
572 >                 " for object %d", index);      
573 >        painCave.isFatal = 1;
574 >        simError();
575 >      }
576 >
577 >      ji = sd->getJ();
578 >
579 >      if (isinf(ji[0]) || isnan(ji[0]) ||
580 >          isinf(ji[1]) || isnan(ji[1]) ||
581 >          isinf(ji[2]) || isnan(ji[2]) ) {      
582 >        sprintf( painCave.errMsg,
583 >                 "DumpWriter detected a numerical error writing the angular"
584 >                 " momentum for object %d", index);      
585 >        painCave.isFatal = 1;
586 >        simError();
587 >      }
588 >
589        sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
590                q[0], q[1], q[2], q[3],
591                ji[0], ji[1], ji[2]);
# Line 333 | Line 593 | namespace OpenMD {
593      }
594  
595      if (needForceVector_) {
596 <      type += "ft";
597 <      Vector3d frc;
598 <      Vector3d trq;
599 <      frc = integrableObject->getFrc();
600 <      trq = integrableObject->getTrq();
601 <              
602 <      sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e",
603 <              frc[0], frc[1], frc[2],
604 <              trq[0], trq[1], trq[2]);
596 >      type += "f";
597 >      Vector3d frc = sd->getFrc();
598 >      if (isinf(frc[0]) || isnan(frc[0]) ||
599 >          isinf(frc[1]) || isnan(frc[1]) ||
600 >          isinf(frc[2]) || isnan(frc[2]) ) {      
601 >        sprintf( painCave.errMsg,
602 >                 "DumpWriter detected a numerical error writing the force"
603 >                 " for object %d", index);      
604 >        painCave.isFatal = 1;
605 >        simError();
606 >      }
607 >      sprintf(tempBuffer, " %13e %13e %13e",
608 >              frc[0], frc[1], frc[2]);
609        line += tempBuffer;
610 +      
611 +      if (sd->isDirectional()) {
612 +        type += "t";
613 +        Vector3d trq = sd->getTrq();        
614 +        if (isinf(trq[0]) || isnan(trq[0]) ||
615 +            isinf(trq[1]) || isnan(trq[1]) ||
616 +            isinf(trq[2]) || isnan(trq[2]) ) {      
617 +          sprintf( painCave.errMsg,
618 +                   "DumpWriter detected a numerical error writing the torque"
619 +                   " for object %d", index);      
620 +          painCave.isFatal = 1;
621 +          simError();
622 +        }        
623 +        sprintf(tempBuffer, " %13e %13e %13e",
624 +                trq[0], trq[1], trq[2]);
625 +        line += tempBuffer;
626 +      }      
627      }
628 <        
628 >
629      sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
630      return std::string(tempBuffer);
631    }
632  
633 +  std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
634 +    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
635 +
636 +    std::string id;
637 +    std::string type;
638 +    std::string line;
639 +    char tempBuffer[4096];
640 +
641 +    if (sd->isRigidBody()) {
642 +      sprintf(tempBuffer, "%10d           ", ioIndex);
643 +      id = std::string(tempBuffer);
644 +    } else {
645 +      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
646 +      id = std::string(tempBuffer);
647 +    }
648 +              
649 +    if (needFlucQ_) {
650 +      if (storageLayout & DataStorage::dslFlucQPosition) {
651 +        type += "c";
652 +        RealType fqPos = sd->getFlucQPos();
653 +        if (isinf(fqPos) || isnan(fqPos) ) {      
654 +          sprintf( painCave.errMsg,
655 +                   "DumpWriter detected a numerical error writing the"
656 +                   " fluctuating charge for object %s", id.c_str());      
657 +          painCave.isFatal = 1;
658 +          simError();
659 +        }
660 +        sprintf(tempBuffer, " %13e ", fqPos);
661 +        line += tempBuffer;
662 +      }
663 +
664 +      if (storageLayout & DataStorage::dslFlucQVelocity) {
665 +        type += "w";    
666 +        RealType fqVel = sd->getFlucQVel();
667 +        if (isinf(fqVel) || isnan(fqVel) ) {      
668 +          sprintf( painCave.errMsg,
669 +                   "DumpWriter detected a numerical error writing the"
670 +                   " fluctuating charge velocity for object %s", id.c_str());      
671 +          painCave.isFatal = 1;
672 +          simError();
673 +        }
674 +        sprintf(tempBuffer, " %13e ", fqVel);
675 +        line += tempBuffer;
676 +      }
677 +
678 +      if (needForceVector_) {
679 +        if (storageLayout & DataStorage::dslFlucQForce) {          
680 +          type += "g";
681 +          RealType fqFrc = sd->getFlucQFrc();        
682 +          if (isinf(fqFrc) || isnan(fqFrc) ) {      
683 +            sprintf( painCave.errMsg,
684 +                     "DumpWriter detected a numerical error writing the"
685 +                     " fluctuating charge force for object %s", id.c_str());      
686 +            painCave.isFatal = 1;
687 +            simError();
688 +          }
689 +          sprintf(tempBuffer, " %13e ", fqFrc);        
690 +          line += tempBuffer;
691 +        }
692 +      }
693 +    }
694 +    
695 +    if (needElectricField_) {
696 +      if (storageLayout & DataStorage::dslElectricField) {
697 +        type += "e";
698 +        Vector3d eField= sd->getElectricField();
699 +        if (isinf(eField[0]) || isnan(eField[0]) ||
700 +            isinf(eField[1]) || isnan(eField[1]) ||
701 +            isinf(eField[2]) || isnan(eField[2]) ) {      
702 +          sprintf( painCave.errMsg,
703 +                   "DumpWriter detected a numerical error writing the electric"
704 +                   " field for object %s", id.c_str());      
705 +          painCave.isFatal = 1;
706 +          simError();
707 +        }
708 +        sprintf(tempBuffer, " %13e %13e %13e",
709 +                eField[0], eField[1], eField[2]);
710 +        line += tempBuffer;
711 +      }
712 +    }
713 +
714 +
715 +    if (needParticlePot_) {
716 +      if (storageLayout & DataStorage::dslParticlePot) {
717 +        type += "u";
718 +        RealType particlePot = sd->getParticlePot();
719 +        if (isinf(particlePot) || isnan(particlePot)) {      
720 +          sprintf( painCave.errMsg,
721 +                   "DumpWriter detected a numerical error writing the particle "
722 +                   " potential for object %s", id.c_str());      
723 +          painCave.isFatal = 1;
724 +          simError();
725 +        }
726 +        sprintf(tempBuffer, " %13e", particlePot);
727 +        line += tempBuffer;
728 +      }
729 +    }
730 +  
731 +    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
732 +    return std::string(tempBuffer);
733 +  }
734 +
735    void DumpWriter::writeDump() {
736      writeFrame(*dumpFile_);
737    }
738  
739    void DumpWriter::writeEor() {
740 <    std::ostream* eorStream;
741 <    
740 >
741 >    std::ostream* eorStream = NULL;
742 >
743   #ifdef IS_MPI
744      if (worldRank == 0) {
745   #endif // is_mpi
746 <
746 >      
747        eorStream = createOStream(eorFilename_);
748  
749   #ifdef IS_MPI
750      }
751 < #endif // is_mpi    
752 <
751 > #endif
752 >    
753      writeFrame(*eorStream);
754 <
754 >      
755   #ifdef IS_MPI
756      if (worldRank == 0) {
757 < #endif // is_mpi
757 > #endif
758 >      
759        writeClosing(*eorStream);
760        delete eorStream;
761 +      
762   #ifdef IS_MPI
763      }
764   #endif // is_mpi  
# Line 386 | Line 772 | namespace OpenMD {
772   #ifdef IS_MPI
773      if (worldRank == 0) {
774   #endif // is_mpi
389
775        buffers.push_back(dumpFile_->rdbuf());
391
776        eorStream = createOStream(eorFilename_);
393
777        buffers.push_back(eorStream->rdbuf());
395        
778   #ifdef IS_MPI
779      }
780   #endif // is_mpi    
781  
782      TeeBuf tbuf(buffers.begin(), buffers.end());
783      std::ostream os(&tbuf);
402
784      writeFrame(os);
785  
786   #ifdef IS_MPI
# Line 409 | Line 790 | namespace OpenMD {
790        delete eorStream;
791   #ifdef IS_MPI
792      }
793 < #endif // is_mpi  
413 <    
793 > #endif // is_mpi      
794    }
795  
796    std::ostream* DumpWriter::createOStream(const std::string& filename) {
797  
798      std::ostream* newOStream;
799 < #ifdef HAVE_LIBZ
799 > #ifdef HAVE_ZLIB
800      if (needCompression_) {
801        newOStream = new ogzstream(filename.c_str());
802      } else {
# Line 426 | Line 806 | namespace OpenMD {
806      newOStream = new std::ofstream(filename.c_str());
807   #endif
808      //write out MetaData first
809 <    (*newOStream) << "<OpenMD version=1>" << std::endl;
809 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
810      (*newOStream) << "  <MetaData>" << std::endl;
811      (*newOStream) << info_->getRawMetaData();
812      (*newOStream) << "  </MetaData>" << std::endl;

Comparing trunk/src/io/DumpWriter.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1983 by gezelter, Tue Apr 15 20:36:19 2014 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines