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 1746 by gezelter, Wed Jun 6 02:18:54 2012 UTC vs.
Revision 1875 by gezelter, Fri May 17 14:41:42 2013 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).          
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   */
# 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 150 | Line 168 | namespace OpenMD {
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_) {
# 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 277 | Line 302 | namespace OpenMD {
302    void DumpWriter::writeFrame(std::ostream& os) {
303  
304   #ifdef IS_MPI
305 <    MPI_Status istatus;
305 >    MPI::Status istatus;
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  
314   #ifndef IS_MPI
315      os << "  <Snapshot>\n";
316  
317      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
292 –
293 –    os << "    <StuntDoubles>\n";
294 –    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
318  
319 +    os << "    <StuntDoubles>\n";
320 +    for (mol = info_->beginMolecule(mi); mol != NULL;
321 +         mol = info_->nextMolecule(mi)) {
322        
323 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
324 <           integrableObject = mol->nextIntegrableObject(ii)) {  
325 <          os << prepareDumpLine(integrableObject);
323 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
324 >           sd = mol->nextIntegrableObject(ii)) {        
325 >          os << prepareDumpLine(sd);
326            
327        }
328      }    
329      os << "    </StuntDoubles>\n";
330 <    
330 >
331 >    if (doSiteData_) {
332 >      os << "    <SiteData>\n";
333 >      for (mol = info_->beginMolecule(mi); mol != NULL;
334 >           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* 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();
360   #else
309 –    //every node prepares the dump lines for integrable objects belong to itself
310 –    std::string buffer;
311 –    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
361  
362 +    const int masterNode = 0;
363 +    int worldRank = MPI::COMM_WORLD.Get_rank();
364 +    int nProc = MPI::COMM_WORLD.Get_size();
365  
366 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
367 <           integrableObject = mol->nextIntegrableObject(ii)) {  
368 <          buffer += prepareDumpLine(integrableObject);
366 >    if (worldRank == masterNode) {      
367 >      os << "  <Snapshot>\n";  
368 >      writeFrameProperties(os,
369 >                           info_->getSnapshotManager()->getCurrentSnapshot());
370 >      os << "    <StuntDoubles>\n";
371 >    }
372 >
373 >    //every node prepares the dump lines for integrable objects belong to itself
374 >    std::string buffer;
375 >    for (mol = info_->beginMolecule(mi); mol != NULL;
376 >         mol = info_->nextMolecule(mi)) {
377 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;
378 >           sd = mol->nextIntegrableObject(ii)) {        
379 >        buffer += prepareDumpLine(sd);
380        }
381      }
382      
320 –    const int masterNode = 0;
321 –    int nProc;
322 –    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
383      if (worldRank == masterNode) {      
324 –      os << "  <Snapshot>\n";  
325 –      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
326 –      os << "    <StuntDoubles>\n";
327 –        
384        os << buffer;
385 <
385 >      
386        for (int i = 1; i < nProc; ++i) {
387 +        // tell processor i to start sending us data:
388 +        MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
389  
390          // receive the length of the string buffer that was
391 <        // prepared by processor i
334 <
335 <        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
391 >        // prepared by processor i:        
392          int recvLength;
393 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
393 >        MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
394 >                             istatus);
395 >
396 >        // create a buffer to receive the data
397          char* recvBuffer = new char[recvLength];
398          if (recvBuffer == NULL) {
399          } else {
400 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
400 >          // receive the data:
401 >          MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
402 >                               MPI::ANY_TAG, istatus);
403 >          // send it to the file:
404            os << recvBuffer;
405 +          // get rid of the receive buffer:
406            delete [] recvBuffer;
407          }
408        }
346 –      os << "    </StuntDoubles>\n";
347 –      
348 –      os << "  </Snapshot>\n";
349 –      os.flush();
409      } else {
410        int sendBufferLength = buffer.size() + 1;
411        int myturn = 0;
412        for (int i = 1; i < nProc; ++i){
413 <        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
413 >        // wait for the master node to call our number:
414 >        MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
415          if (myturn == worldRank){
416 <          MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
417 <          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
416 >          // send the length of our buffer:
417 >          MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
418 >
419 >          // send our buffer:
420 >          MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
421 >                               MPI::CHAR, masterNode, 0);
422          }
423        }
424      }
425 +    
426 +    if (worldRank == masterNode) {      
427 +      os << "    </StuntDoubles>\n";
428 +    }
429  
430 < #endif // is_mpi
430 >    if (doSiteData_) {
431 >      if (worldRank == masterNode) {
432 >        os << "    <SiteData>\n";
433 >      }
434 >      buffer.clear();
435 >      for (mol = info_->beginMolecule(mi); mol != NULL;
436 >           mol = info_->nextMolecule(mi)) {
437 >              
438 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
439 >             sd = mol->nextIntegrableObject(ii)) {      
440 >          
441 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
442 >          // do one for the IO itself
443 >          buffer += prepareSiteLine(sd, ioIndex, 0);
444  
445 +          if (sd->isRigidBody()) {
446 +            
447 +            RigidBody* rb = static_cast<RigidBody*>(sd);
448 +            int siteIndex = 0;
449 +            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
450 +                 atom = rb->nextAtom(ai)) {                                            
451 +              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
452 +              siteIndex++;
453 +            }
454 +          }
455 +        }
456 +      }
457 +
458 +      if (worldRank == masterNode) {    
459 +        os << buffer;
460 +        
461 +        for (int i = 1; i < nProc; ++i) {
462 +          
463 +          // tell processor i to start sending us data:
464 +          MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
465 +          
466 +          // receive the length of the string buffer that was
467 +          // prepared by processor i:        
468 +          int recvLength;
469 +          MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
470 +                               istatus);
471 +          
472 +          // create a buffer to receive the data
473 +          char* recvBuffer = new char[recvLength];
474 +          if (recvBuffer == NULL) {
475 +          } else {
476 +            // receive the data:
477 +            MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
478 +                                 MPI::ANY_TAG, istatus);
479 +            // send it to the file:
480 +            os << recvBuffer;
481 +            // get rid of the receive buffer:
482 +            delete [] recvBuffer;
483 +          }
484 +        }      
485 +      } else {
486 +        int sendBufferLength = buffer.size() + 1;
487 +        int myturn = 0;
488 +        for (int i = 1; i < nProc; ++i){
489 +          // wait for the master node to call our number:
490 +          MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
491 +          if (myturn == worldRank){
492 +            // send the length of our buffer:
493 +            MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
494 +            // send our buffer:
495 +            MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
496 +                                 MPI::CHAR, masterNode, 0);
497 +          }
498 +        }
499 +      }
500 +      
501 +      if (worldRank == masterNode) {    
502 +        os << "    </SiteData>\n";
503 +      }
504 +    }
505 +    
506 +    if (worldRank == masterNode) {
507 +      os << "  </Snapshot>\n";
508 +      os.flush();
509 +    }
510 +    
511 + #endif // is_mpi
512 +    
513    }
514  
515 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
515 >  std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
516          
517 <    int index = integrableObject->getGlobalIntegrableObjectIndex();
517 >    int index = sd->getGlobalIntegrableObjectIndex();
518      std::string type("pv");
519      std::string line;
520      char tempBuffer[4096];
521  
522      Vector3d pos;
523      Vector3d vel;
524 <    pos = integrableObject->getPos();
524 >    pos = sd->getPos();
525  
526      if (isinf(pos[0]) || isnan(pos[0]) ||
527          isinf(pos[1]) || isnan(pos[1]) ||
# Line 384 | Line 533 | namespace OpenMD {
533        simError();
534      }
535  
536 <    vel = integrableObject->getVel();          
536 >    vel = sd->getVel();        
537  
538      if (isinf(vel[0]) || isnan(vel[0]) ||
539          isinf(vel[1]) || isnan(vel[1]) ||
# Line 401 | Line 550 | namespace OpenMD {
550              vel[0], vel[1], vel[2]);                    
551      line += tempBuffer;
552  
553 <    if (integrableObject->isDirectional()) {
553 >    if (sd->isDirectional()) {
554        type += "qj";
555        Quat4d q;
556        Vector3d ji;
557 <      q = integrableObject->getQ();
557 >      q = sd->getQ();
558  
559        if (isinf(q[0]) || isnan(q[0]) ||
560            isinf(q[1]) || isnan(q[1]) ||
# Line 418 | Line 567 | namespace OpenMD {
567          simError();
568        }
569  
570 <      ji = integrableObject->getJ();
570 >      ji = sd->getJ();
571  
572        if (isinf(ji[0]) || isnan(ji[0]) ||
573            isinf(ji[1]) || isnan(ji[1]) ||
# Line 438 | Line 587 | namespace OpenMD {
587  
588      if (needForceVector_) {
589        type += "f";
590 <      Vector3d frc = integrableObject->getFrc();
590 >      Vector3d frc = sd->getFrc();
591        if (isinf(frc[0]) || isnan(frc[0]) ||
592            isinf(frc[1]) || isnan(frc[1]) ||
593            isinf(frc[2]) || isnan(frc[2]) ) {      
# Line 452 | Line 601 | namespace OpenMD {
601                frc[0], frc[1], frc[2]);
602        line += tempBuffer;
603        
604 <      if (integrableObject->isDirectional()) {
604 >      if (sd->isDirectional()) {
605          type += "t";
606 <        Vector3d trq = integrableObject->getTrq();        
606 >        Vector3d trq = sd->getTrq();        
607          if (isinf(trq[0]) || isnan(trq[0]) ||
608              isinf(trq[1]) || isnan(trq[1]) ||
609              isinf(trq[2]) || isnan(trq[2]) ) {      
# Line 470 | Line 619 | namespace OpenMD {
619        }      
620      }
621  
622 <    if (needParticlePot_) {
623 <      type += "u";
624 <      RealType particlePot = integrableObject->getParticlePot();
625 <      if (isinf(particlePot) || isnan(particlePot)) {      
626 <        sprintf( painCave.errMsg,
627 <                 "DumpWriter detected a numerical error writing the particle "
628 <                 " potential for object %d", index);      
629 <        painCave.isFatal = 1;
630 <        simError();
631 <      }
632 <      sprintf(tempBuffer, " %13e", particlePot);
633 <      line += tempBuffer;
622 >    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
623 >    return std::string(tempBuffer);
624 >  }
625 >
626 >  std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
627 >        
628 >
629 >    std::string id;
630 >    std::string type;
631 >    std::string line;
632 >    char tempBuffer[4096];
633 >
634 >    if (sd->isRigidBody()) {
635 >      sprintf(tempBuffer, "%10d           ", ioIndex);
636 >      id = std::string(tempBuffer);
637 >    } else {
638 >      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
639 >      id = std::string(tempBuffer);
640      }
641 <    
641 >              
642      if (needFlucQ_) {
643        type += "cw";
644 <      RealType fqPos = integrableObject->getFlucQPos();
644 >      RealType fqPos = sd->getFlucQPos();
645        if (isinf(fqPos) || isnan(fqPos) ) {      
646          sprintf( painCave.errMsg,
647                   "DumpWriter detected a numerical error writing the"
648 <                 " fluctuating charge for object %d", index);      
648 >                 " fluctuating charge for object %s", id.c_str());      
649          painCave.isFatal = 1;
650          simError();
651        }
652        sprintf(tempBuffer, " %13e ", fqPos);
653        line += tempBuffer;
654      
655 <      RealType fqVel = integrableObject->getFlucQVel();
655 >      RealType fqVel = sd->getFlucQVel();
656        if (isinf(fqVel) || isnan(fqVel) ) {      
657          sprintf( painCave.errMsg,
658                   "DumpWriter detected a numerical error writing the"
659 <                 " fluctuating charge velocity for object %d", index);      
659 >                 " fluctuating charge velocity for object %s", id.c_str());      
660          painCave.isFatal = 1;
661          simError();
662        }
# Line 510 | Line 665 | namespace OpenMD {
665  
666        if (needForceVector_) {
667          type += "g";
668 <        RealType fqFrc = integrableObject->getFlucQFrc();        
668 >        RealType fqFrc = sd->getFlucQFrc();        
669          if (isinf(fqFrc) || isnan(fqFrc) ) {      
670            sprintf( painCave.errMsg,
671                     "DumpWriter detected a numerical error writing the"
672 <                   " fluctuating charge force for object %d", index);      
672 >                   " fluctuating charge force for object %s", id.c_str());      
673            painCave.isFatal = 1;
674            simError();
675          }
# Line 525 | Line 680 | namespace OpenMD {
680  
681      if (needElectricField_) {
682        type += "e";
683 <      Vector3d eField= integrableObject->getElectricField();
683 >      Vector3d eField= sd->getElectricField();
684        if (isinf(eField[0]) || isnan(eField[0]) ||
685            isinf(eField[1]) || isnan(eField[1]) ||
686            isinf(eField[2]) || isnan(eField[2]) ) {      
687          sprintf( painCave.errMsg,
688                   "DumpWriter detected a numerical error writing the electric"
689 <                 " field for object %d", index);      
689 >                 " field for object %s", id.c_str());      
690          painCave.isFatal = 1;
691          simError();
692        }
# Line 540 | Line 695 | namespace OpenMD {
695        line += tempBuffer;
696      }
697  
698 <    sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
698 >
699 >    if (needParticlePot_) {
700 >      type += "u";
701 >      RealType particlePot = sd->getParticlePot();
702 >      if (isinf(particlePot) || isnan(particlePot)) {      
703 >        sprintf( painCave.errMsg,
704 >                 "DumpWriter detected a numerical error writing the particle "
705 >                 " potential for object %s", id.c_str());      
706 >        painCave.isFatal = 1;
707 >        simError();
708 >      }
709 >      sprintf(tempBuffer, " %13e", particlePot);
710 >      line += tempBuffer;
711 >    }
712 >    
713 >
714 >    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
715      return std::string(tempBuffer);
716    }
717  
# Line 554 | Line 725 | namespace OpenMD {
725   #ifdef IS_MPI
726      if (worldRank == 0) {
727   #endif // is_mpi
728 <
728 >      
729        eorStream = createOStream(eorFilename_);
730 <
730 >      writeFrame(*eorStream);
731 >      
732   #ifdef IS_MPI
733      }
562 – #endif // is_mpi    
563 –
564 –    writeFrame(*eorStream);
565 –
566 – #ifdef IS_MPI
734      if (worldRank == 0) {
735   #endif // is_mpi
736 +
737        writeClosing(*eorStream);
738        delete eorStream;
739 +
740   #ifdef IS_MPI
741      }
742   #endif // is_mpi  
# Line 611 | Line 780 | namespace OpenMD {
780    std::ostream* DumpWriter::createOStream(const std::string& filename) {
781  
782      std::ostream* newOStream;
783 < #ifdef HAVE_LIBZ
783 > #ifdef HAVE_ZLIB
784      if (needCompression_) {
785        newOStream = new ogzstream(filename.c_str());
786      } else {

Diff Legend

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