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 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

# Line 36 | Line 36
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).                        
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 "io/DumpWriter.hpp"
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 58 | Line 65 | namespace OpenMD {
65      : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
66  
67      Globals* simParams = info->getSimParams();
68 <    needCompression_ = simParams->getCompressDumpFile();
69 <    needForceVector_ = simParams->getOutputForceVector();
70 <    needParticlePot_ = simParams->getOutputParticlePotential();
68 >    needCompression_   = simParams->getCompressDumpFile();
69 >    needForceVector_   = simParams->getOutputForceVector();
70 >    needParticlePot_   = simParams->getOutputParticlePotential();
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 98 | Line 114 | namespace OpenMD {
114      Globals* simParams = info->getSimParams();
115      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
116  
117 <    needCompression_ = simParams->getCompressDumpFile();
118 <    needForceVector_ = simParams->getOutputForceVector();
119 <    needParticlePot_ = simParams->getOutputParticlePotential();
117 >    needCompression_   = simParams->getCompressDumpFile();
118 >    needForceVector_   = simParams->getOutputForceVector();
119 >    needParticlePot_   = simParams->getOutputParticlePotential();
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 138 | Line 163 | namespace OpenMD {
163      Globals* simParams = info->getSimParams();
164      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
165      
166 <    needCompression_ = simParams->getCompressDumpFile();
167 <    needForceVector_ = simParams->getOutputForceVector();
168 <    needParticlePot_ = simParams->getOutputParticlePotential();
169 <    
166 >    needCompression_   = simParams->getCompressDumpFile();
167 >    needForceVector_   = simParams->getOutputForceVector();
168 >    needParticlePot_   = simParams->getOutputParticlePotential();
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 230 | 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 272 | 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 285 | 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 302 | 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 318 | Line 379 | namespace OpenMD {
379          
380        os << buffer;
381  
321 –    
382        for (int i = 1; i < nProc; ++i) {
383  
384          // receive the length of the string buffer that was
385          // prepared by processor i
386 <      
386 >
387          MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
388          int recvLength;
389          MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
# Line 355 | 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 376 | 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 393 | 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 410 | 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 430 | Line 490 | namespace OpenMD {
490  
491      if (needForceVector_) {
492        type += "f";
493 <      Vector3d frc;
434 <
435 <      frc = integrableObject->getFrc();
436 <
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 447 | 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;
453 <        
454 <        trq = integrableObject->getTrq();
455 <        
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 461 | Line 515 | namespace OpenMD {
515                     " for object %d", index);      
516            painCave.isFatal = 1;
517            simError();
518 <        }
465 <        
518 >        }        
519          sprintf(tempBuffer, " %13e %13e %13e",
520                  trq[0], trq[1], trq[2]);
521          line += tempBuffer;
522        }      
523      }
471 –    if (needParticlePot_) {
472 –      type += "u";
473 –      RealType particlePot;
524  
525 <      particlePot = integrableObject->getParticlePot();
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 +              
545 +    if (needFlucQ_) {
546 +      type += "cw";
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 %s", id.c_str());      
552 +        painCave.isFatal = 1;
553 +        simError();
554 +      }
555 +      sprintf(tempBuffer, " %13e ", fqPos);
556 +      line += tempBuffer;
557 +    
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 %s", id.c_str());      
563 +        painCave.isFatal = 1;
564 +        simError();
565 +      }
566 +      sprintf(tempBuffer, " %13e ", fqVel);
567 +      line += tempBuffer;
568 +
569 +      if (needForceVector_) {
570 +        type += "g";
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 %s", id.c_str());      
576 +          painCave.isFatal = 1;
577 +          simError();
578 +        }
579 +        sprintf(tempBuffer, " %13e ", fqFrc);        
580 +        line += tempBuffer;
581 +      }
582 +    }
583 +
584 +    if (needElectricField_) {
585 +      type += "e";
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 %s", id.c_str());      
593 +        painCave.isFatal = 1;
594 +        simError();
595 +      }
596 +      sprintf(tempBuffer, " %13e %13e %13e",
597 +              eField[0], eField[1], eField[2]);
598 +      line += tempBuffer;
599 +    }
600 +
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 %d", index);      
608 >                 " potential for object %s", id.c_str());      
609          painCave.isFatal = 1;
610          simError();
611        }
# Line 485 | Line 613 | namespace OpenMD {
613        line += tempBuffer;
614      }
615      
616 <    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
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 556 | 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 566 | 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