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

Comparing branches/development/src/io/DumpWriter.cpp (file contents):
Revision 1714 by gezelter, Sat May 19 18:12:46 2012 UTC vs.
Revision 1769 by gezelter, Mon Jul 9 14:15:52 2012 UTC

# Line 44 | Line 44
44   #include "primitives/Molecule.hpp"
45   #include "utils/simError.h"
46   #include "io/basic_teebuf.hpp"
47 + #ifdef HAVE_ZLIB
48   #include "io/gzstream.hpp"
49 + #endif
50   #include "io/Globals.hpp"
51  
52 + #ifdef _MSC_VER
53 + #define isnan(x) _isnan((x))
54 + #define isinf(x) (!_finite(x) && !_isnan(x))
55 + #endif
56  
57   #ifdef IS_MPI
58   #include <mpi.h>
59 < #endif //is_mpi
59 > #endif
60  
61   using namespace std;
62   namespace OpenMD {
# Line 65 | Line 71 | namespace OpenMD {
71      needFlucQ_         = simParams->getOutputFluctuatingCharges();
72      needElectricField_ = simParams->getOutputElectricField();
73  
74 +    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
75 +      doSiteData_ = true;
76 +    } else {
77 +      doSiteData_ = false;
78 +    }
79 +
80      createDumpFile_ = true;
81   #ifdef HAVE_LIBZ
82      if (needCompression_) {
# Line 108 | Line 120 | namespace OpenMD {
120      needFlucQ_         = simParams->getOutputFluctuatingCharges();
121      needElectricField_ = simParams->getOutputElectricField();
122  
123 +    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
124 +      doSiteData_ = true;
125 +    } else {
126 +      doSiteData_ = false;
127 +    }
128 +
129      createDumpFile_ = true;
130   #ifdef HAVE_LIBZ
131      if (needCompression_) {
# Line 151 | Line 169 | namespace OpenMD {
169      needFlucQ_         = simParams->getOutputFluctuatingCharges();
170      needElectricField_ = simParams->getOutputElectricField();
171  
172 +    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
173 +      doSiteData_ = true;
174 +    } else {
175 +      doSiteData_ = false;
176 +    }
177 +
178   #ifdef HAVE_LIBZ
179      if (needCompression_) {
180        filename_ += ".gz";
# Line 239 | Line 263 | namespace OpenMD {
263              hmat(0, 2), hmat(1, 2), hmat(2, 2));
264      os << buffer;
265  
266 <    RealType chi = s->getChi();
267 <    RealType integralOfChiDt = s->getIntegralOfChiDt();
268 <    if (isinf(chi) || isnan(chi) ||
269 <        isinf(integralOfChiDt) || isnan(integralOfChiDt)) {      
266 >    pair<RealType, RealType> thermostat = s->getThermostat();
267 >
268 >    if (isinf(thermostat.first)  || isnan(thermostat.first) ||
269 >        isinf(thermostat.second) || isnan(thermostat.second)) {      
270        sprintf( painCave.errMsg,
271                 "DumpWriter detected a numerical error writing the thermostat");
272        painCave.isFatal = 1;
273        simError();
274      }
275 <    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
275 >    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", thermostat.first,
276 >            thermostat.second);
277      os << buffer;
278  
279      Mat3x3d eta;
280 <    eta = s->getEta();
280 >    eta = s->getBarostat();
281  
282      for (unsigned int i = 0; i < 3; i++) {
283        for (unsigned int j = 0; j < 3; j++) {
# Line 281 | Line 306 | namespace OpenMD {
306   #endif
307  
308      Molecule* mol;
309 <    StuntDouble* integrableObject;
309 >    StuntDouble* sd;
310      SimInfo::MoleculeIterator mi;
311      Molecule::IntegrableObjectIterator ii;
312 +    RigidBody::AtomIterator ai;
313 +    Atom* atom;
314  
315   #ifndef IS_MPI
316      os << "  <Snapshot>\n";
# Line 294 | Line 321 | namespace OpenMD {
321      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
322  
323        
324 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
325 <           integrableObject = mol->nextIntegrableObject(ii)) {  
326 <          os << prepareDumpLine(integrableObject);
324 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
325 >           sd = mol->nextIntegrableObject(ii)) {        
326 >          os << prepareDumpLine(sd);
327            
328        }
329      }    
330      os << "    </StuntDoubles>\n";
331 <    
331 >
332 >    if (doSiteData_) {
333 >      os << "    <SiteData>\n";
334 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
335 >              
336 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
337 >           sd = mol->nextIntegrableObject(ii)) {        
338 >
339 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
340 >          // do one for the IO itself
341 >          os << prepareSiteLine(sd, ioIndex, 0);
342 >
343 >          if (sd->isRigidBody()) {
344 >            
345 >            RigidBody* rb = static_cast<RigidBody*>(sd);
346 >            int siteIndex = 0;
347 >            for (atom = rb->beginAtom(ai); atom != NULL;  
348 >                 atom = rb->nextAtom(ai)) {                                            
349 >              os << prepareSiteLine(atom, ioIndex, siteIndex);
350 >              siteIndex++;
351 >            }
352 >          }
353 >        }
354 >      }    
355 >      os << "    </SiteData>\n";
356 >    }
357      os << "  </Snapshot>\n";
358  
359      os.flush();
# Line 311 | Line 363 | namespace OpenMD {
363      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
364  
365  
366 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
367 <           integrableObject = mol->nextIntegrableObject(ii)) {  
368 <          buffer += prepareDumpLine(integrableObject);
366 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;
367 >           sd = mol->nextIntegrableObject(ii)) {        
368 >          buffer += prepareDumpLine(sd);
369        }
370      }
371      
# Line 363 | Line 415 | namespace OpenMD {
415  
416    }
417  
418 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
418 >  std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
419          
420 <    int index = integrableObject->getGlobalIntegrableObjectIndex();
420 >    int index = sd->getGlobalIntegrableObjectIndex();
421      std::string type("pv");
422      std::string line;
423      char tempBuffer[4096];
424  
425      Vector3d pos;
426      Vector3d vel;
427 <    pos = integrableObject->getPos();
427 >    pos = sd->getPos();
428  
429      if (isinf(pos[0]) || isnan(pos[0]) ||
430          isinf(pos[1]) || isnan(pos[1]) ||
# Line 384 | Line 436 | namespace OpenMD {
436        simError();
437      }
438  
439 <    vel = integrableObject->getVel();          
439 >    vel = sd->getVel();        
440  
441      if (isinf(vel[0]) || isnan(vel[0]) ||
442          isinf(vel[1]) || isnan(vel[1]) ||
# Line 401 | Line 453 | namespace OpenMD {
453              vel[0], vel[1], vel[2]);                    
454      line += tempBuffer;
455  
456 <    if (integrableObject->isDirectional()) {
456 >    if (sd->isDirectional()) {
457        type += "qj";
458        Quat4d q;
459        Vector3d ji;
460 <      q = integrableObject->getQ();
460 >      q = sd->getQ();
461  
462        if (isinf(q[0]) || isnan(q[0]) ||
463            isinf(q[1]) || isnan(q[1]) ||
# Line 418 | Line 470 | namespace OpenMD {
470          simError();
471        }
472  
473 <      ji = integrableObject->getJ();
473 >      ji = sd->getJ();
474  
475        if (isinf(ji[0]) || isnan(ji[0]) ||
476            isinf(ji[1]) || isnan(ji[1]) ||
# Line 438 | Line 490 | namespace OpenMD {
490  
491      if (needForceVector_) {
492        type += "f";
493 <      Vector3d frc = integrableObject->getFrc();
493 >      Vector3d frc = sd->getFrc();
494        if (isinf(frc[0]) || isnan(frc[0]) ||
495            isinf(frc[1]) || isnan(frc[1]) ||
496            isinf(frc[2]) || isnan(frc[2]) ) {      
# Line 452 | Line 504 | namespace OpenMD {
504                frc[0], frc[1], frc[2]);
505        line += tempBuffer;
506        
507 <      if (integrableObject->isDirectional()) {
507 >      if (sd->isDirectional()) {
508          type += "t";
509 <        Vector3d trq = integrableObject->getTrq();        
509 >        Vector3d trq = sd->getTrq();        
510          if (isinf(trq[0]) || isnan(trq[0]) ||
511              isinf(trq[1]) || isnan(trq[1]) ||
512              isinf(trq[2]) || isnan(trq[2]) ) {      
# Line 470 | Line 522 | namespace OpenMD {
522        }      
523      }
524  
525 <    if (needParticlePot_) {
526 <      type += "u";
527 <      RealType particlePot = integrableObject->getParticlePot();
528 <      if (isinf(particlePot) || isnan(particlePot)) {      
529 <        sprintf( painCave.errMsg,
530 <                 "DumpWriter detected a numerical error writing the particle "
531 <                 " potential for object %d", index);      
532 <        painCave.isFatal = 1;
533 <        simError();
534 <      }
535 <      sprintf(tempBuffer, " %13e", particlePot);
536 <      line += tempBuffer;
525 >    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
526 >    return std::string(tempBuffer);
527 >  }
528 >
529 >  std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
530 >        
531 >
532 >    std::string id;
533 >    std::string type;
534 >    std::string line;
535 >    char tempBuffer[4096];
536 >
537 >    if (sd->isRigidBody()) {
538 >      sprintf(tempBuffer, "%10d           ", ioIndex);
539 >      id = std::string(tempBuffer);
540 >    } else {
541 >      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
542 >      id = std::string(tempBuffer);
543      }
544 <    
544 >              
545      if (needFlucQ_) {
546        type += "cw";
547 <      RealType fqPos = integrableObject->getFlucQPos();
547 >      RealType fqPos = sd->getFlucQPos();
548        if (isinf(fqPos) || isnan(fqPos) ) {      
549          sprintf( painCave.errMsg,
550                   "DumpWriter detected a numerical error writing the"
551 <                 " fluctuating charge for object %d", index);      
551 >                 " fluctuating charge for object %s", id.c_str());      
552          painCave.isFatal = 1;
553          simError();
554        }
555        sprintf(tempBuffer, " %13e ", fqPos);
556        line += tempBuffer;
557      
558 <      RealType fqVel = integrableObject->getFlucQVel();
558 >      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 %d", index);      
562 >                 " fluctuating charge velocity for object %s", id.c_str());      
563          painCave.isFatal = 1;
564          simError();
565        }
# Line 510 | Line 568 | namespace OpenMD {
568  
569        if (needForceVector_) {
570          type += "g";
571 <        RealType fqFrc = integrableObject->getFlucQFrc();        
571 >        RealType fqFrc = sd->getFlucQFrc();        
572          if (isinf(fqFrc) || isnan(fqFrc) ) {      
573            sprintf( painCave.errMsg,
574                     "DumpWriter detected a numerical error writing the"
575 <                   " fluctuating charge force for object %d", index);      
575 >                   " fluctuating charge force for object %s", id.c_str());      
576            painCave.isFatal = 1;
577            simError();
578          }
# Line 525 | Line 583 | namespace OpenMD {
583  
584      if (needElectricField_) {
585        type += "e";
586 <      Vector3d eField= integrableObject->getElectricField();
586 >      Vector3d eField= sd->getElectricField();
587        if (isinf(eField[0]) || isnan(eField[0]) ||
588            isinf(eField[1]) || isnan(eField[1]) ||
589            isinf(eField[2]) || isnan(eField[2]) ) {      
590          sprintf( painCave.errMsg,
591                   "DumpWriter detected a numerical error writing the electric"
592 <                 " field for object %d", index);      
592 >                 " field for object %s", id.c_str());      
593          painCave.isFatal = 1;
594          simError();
595        }
# Line 540 | Line 598 | namespace OpenMD {
598        line += tempBuffer;
599      }
600  
601 <    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
601 >
602 >    if (needParticlePot_) {
603 >      type += "u";
604 >      RealType particlePot = sd->getParticlePot();
605 >      if (isinf(particlePot) || isnan(particlePot)) {      
606 >        sprintf( painCave.errMsg,
607 >                 "DumpWriter detected a numerical error writing the particle "
608 >                 " potential for object %s", id.c_str());      
609 >        painCave.isFatal = 1;
610 >        simError();
611 >      }
612 >      sprintf(tempBuffer, " %13e", particlePot);
613 >      line += tempBuffer;
614 >    }
615 >    
616 >
617 >    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
618      return std::string(tempBuffer);
619    }
620  
# Line 611 | Line 685 | namespace OpenMD {
685    std::ostream* DumpWriter::createOStream(const std::string& filename) {
686  
687      std::ostream* newOStream;
688 < #ifdef HAVE_LIBZ
688 > #ifdef HAVE_ZLIB
689      if (needCompression_) {
690        newOStream = new ogzstream(filename.c_str());
691      } else {
# Line 621 | Line 695 | namespace OpenMD {
695      newOStream = new std::ofstream(filename.c_str());
696   #endif
697      //write out MetaData first
698 <    (*newOStream) << "<OpenMD version=1>" << std::endl;
698 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
699      (*newOStream) << "  <MetaData>" << std::endl;
700      (*newOStream) << info_->getRawMetaData();
701      (*newOStream) << "  </MetaData>" << std::endl;

Diff Legend

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