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

Comparing:
trunk/src/io/DumpWriter.cpp (file contents), Revision 1025 by gezelter, Wed Aug 30 20:33:44 2006 UTC vs.
branches/development/src/io/DumpWriter.cpp (file contents), Revision 1752 by gezelter, Sun Jun 10 14:05:02 2012 UTC

# Line 1 | Line 1
1   /*
2 < * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
2 > * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
5   * non-exclusive, royalty free, license to use, modify and
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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]  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"
# Line 46 | Line 47
47   #include "io/gzstream.hpp"
48   #include "io/Globals.hpp"
49  
50 +
51   #ifdef IS_MPI
52   #include <mpi.h>
53   #endif //is_mpi
54  
55 < namespace oopse {
55 > using namespace std;
56 > namespace OpenMD {
57  
58    DumpWriter::DumpWriter(SimInfo* info)
59      : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
60  
61      Globals* simParams = info->getSimParams();
62 <    needCompression_ = simParams->getCompressDumpFile();
63 <    needForceVector_ = simParams->getOutputForceVector();
62 >    needCompression_   = simParams->getCompressDumpFile();
63 >    needForceVector_   = simParams->getOutputForceVector();
64 >    needParticlePot_   = simParams->getOutputParticlePotential();
65 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
66 >    needElectricField_ = simParams->getOutputElectricField();
67 >
68 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
69 >      doSiteData_ = true;
70 >    } else {
71 >      doSiteData_ = false;
72 >    }
73 >
74      createDumpFile_ = true;
75   #ifdef HAVE_LIBZ
76      if (needCompression_) {
# Line 95 | Line 108 | namespace oopse {
108      Globals* simParams = info->getSimParams();
109      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
110  
111 <    needCompression_ = simParams->getCompressDumpFile();
112 <    needForceVector_ = simParams->getOutputForceVector();
111 >    needCompression_   = simParams->getCompressDumpFile();
112 >    needForceVector_   = simParams->getOutputForceVector();
113 >    needParticlePot_   = simParams->getOutputParticlePotential();
114 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
115 >    needElectricField_ = simParams->getOutputElectricField();
116 >
117 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
118 >      doSiteData_ = true;
119 >    } else {
120 >      doSiteData_ = false;
121 >    }
122 >
123      createDumpFile_ = true;
124   #ifdef HAVE_LIBZ
125      if (needCompression_) {
# Line 134 | Line 157 | namespace oopse {
157      Globals* simParams = info->getSimParams();
158      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
159      
160 <    needCompression_ = simParams->getCompressDumpFile();
161 <    needForceVector_ = simParams->getOutputForceVector();
162 <    
160 >    needCompression_   = simParams->getCompressDumpFile();
161 >    needForceVector_   = simParams->getOutputForceVector();
162 >    needParticlePot_   = simParams->getOutputParticlePotential();
163 >    needFlucQ_         = simParams->getOutputFluctuatingCharges();
164 >    needElectricField_ = simParams->getOutputElectricField();
165 >
166 >    if (needParticlePot_ || needFlucQ_ || needElectricField_) {
167 >      doSiteData_ = true;
168 >    } else {
169 >      doSiteData_ = false;
170 >    }
171 >
172   #ifdef HAVE_LIBZ
173      if (needCompression_) {
174        filename_ += ".gz";
# Line 194 | Line 226 | namespace oopse {
226      os << "    <FrameData>\n";
227  
228      RealType currentTime = s->getTime();
229 +
230 +    if (isinf(currentTime) || isnan(currentTime)) {      
231 +      sprintf( painCave.errMsg,
232 +               "DumpWriter detected a numerical error writing the time");      
233 +      painCave.isFatal = 1;
234 +      simError();
235 +    }
236 +    
237      sprintf(buffer, "        Time: %.10g\n", currentTime);
238      os << buffer;
239  
240      Mat3x3d hmat;
241      hmat = s->getHmat();
242 +
243 +    for (unsigned int i = 0; i < 3; i++) {
244 +      for (unsigned int j = 0; j < 3; j++) {
245 +        if (isinf(hmat(i,j)) || isnan(hmat(i,j))) {      
246 +          sprintf( painCave.errMsg,
247 +                   "DumpWriter detected a numerical error writing the box");
248 +          painCave.isFatal = 1;
249 +          simError();
250 +        }        
251 +      }
252 +    }
253 +    
254      sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
255              hmat(0, 0), hmat(1, 0), hmat(2, 0),
256              hmat(0, 1), hmat(1, 1), hmat(2, 1),
# Line 207 | Line 259 | namespace oopse {
259  
260      RealType chi = s->getChi();
261      RealType integralOfChiDt = s->getIntegralOfChiDt();
262 +    if (isinf(chi) || isnan(chi) ||
263 +        isinf(integralOfChiDt) || isnan(integralOfChiDt)) {      
264 +      sprintf( painCave.errMsg,
265 +               "DumpWriter detected a numerical error writing the thermostat");
266 +      painCave.isFatal = 1;
267 +      simError();
268 +    }
269      sprintf(buffer, "  Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
270      os << buffer;
271  
272      Mat3x3d eta;
273      eta = s->getEta();
274 +
275 +    for (unsigned int i = 0; i < 3; i++) {
276 +      for (unsigned int j = 0; j < 3; j++) {
277 +        if (isinf(eta(i,j)) || isnan(eta(i,j))) {      
278 +          sprintf( painCave.errMsg,
279 +                   "DumpWriter detected a numerical error writing the barostat");
280 +          painCave.isFatal = 1;
281 +          simError();
282 +        }        
283 +      }
284 +    }
285 +
286      sprintf(buffer, "    Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
287              eta(0, 0), eta(1, 0), eta(2, 0),
288              eta(0, 1), eta(1, 1), eta(2, 1),
# Line 231 | Line 302 | namespace oopse {
302      StuntDouble* integrableObject;
303      SimInfo::MoleculeIterator mi;
304      Molecule::IntegrableObjectIterator ii;
305 +    RigidBody::AtomIterator ai;
306 +    Atom* atom;
307  
308   #ifndef IS_MPI
309      os << "  <Snapshot>\n";
# Line 240 | Line 313 | namespace oopse {
313      os << "    <StuntDoubles>\n";
314      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
315  
316 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
316 >      
317 >      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
318             integrableObject = mol->nextIntegrableObject(ii)) {  
319 <        os << prepareDumpLine(integrableObject);
320 <
319 >          os << prepareDumpLine(integrableObject);
320 >          
321        }
322      }    
323      os << "    </StuntDoubles>\n";
324  
325 +    if (doSiteData_) {
326 +      os << "    <SiteData>\n";
327 +      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
328 +              
329 +        for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
330 +           integrableObject = mol->nextIntegrableObject(ii)) {  
331 +
332 +          int ioIndex = integrableObject->getGlobalIntegrableObjectIndex();
333 +          // do one for the IO itself
334 +          os << prepareSiteLine(integrableObject, ioIndex, 0);
335 +
336 +          if (integrableObject->isRigidBody()) {
337 +            
338 +            RigidBody* rb = static_cast<RigidBody*>(integrableObject);
339 +            int siteIndex = 0;
340 +            for (atom = rb->beginAtom(ai); atom != NULL;  
341 +                 atom = rb->nextAtom(ai)) {                                            
342 +              os << prepareSiteLine(atom, ioIndex, siteIndex);
343 +              siteIndex++;
344 +            }
345 +          }
346 +        }
347 +      }    
348 +      os << "    </SiteData>\n";
349 +    }
350      os << "  </Snapshot>\n";
351  
352      os.flush();
# Line 255 | Line 354 | namespace oopse {
354      //every node prepares the dump lines for integrable objects belong to itself
355      std::string buffer;
356      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
357 <      
357 >
358 >
359        for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
360             integrableObject = mol->nextIntegrableObject(ii)) {  
361 <        buffer += prepareDumpLine(integrableObject);
361 >          buffer += prepareDumpLine(integrableObject);
362        }
363      }
364      
365      const int masterNode = 0;
366 <
366 >    int nProc;
367 >    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
368      if (worldRank == masterNode) {      
369        os << "  <Snapshot>\n";  
370        writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
# Line 271 | Line 372 | namespace oopse {
372          
373        os << buffer;
374  
274      int nProc;
275      MPI_Comm_size(MPI_COMM_WORLD, &nProc);
375        for (int i = 1; i < nProc; ++i) {
376  
377          // receive the length of the string buffer that was
378          // prepared by processor i
379  
380 +        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
381          int recvLength;
382          MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
383          char* recvBuffer = new char[recvLength];
# Line 285 | Line 385 | namespace oopse {
385          } else {
386            MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
387            os << recvBuffer;
388 <          delete recvBuffer;
388 >          delete [] recvBuffer;
389          }
390        }
391        os << "    </StuntDoubles>\n";
# Line 294 | Line 394 | namespace oopse {
394        os.flush();
395      } else {
396        int sendBufferLength = buffer.size() + 1;
397 <      MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
398 <      MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
397 >      int myturn = 0;
398 >      for (int i = 1; i < nProc; ++i){
399 >        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
400 >        if (myturn == worldRank){
401 >          MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
402 >          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
403 >        }
404 >      }
405      }
406  
407   #endif // is_mpi
# Line 312 | Line 418 | namespace oopse {
418      Vector3d pos;
419      Vector3d vel;
420      pos = integrableObject->getPos();
421 +
422 +    if (isinf(pos[0]) || isnan(pos[0]) ||
423 +        isinf(pos[1]) || isnan(pos[1]) ||
424 +        isinf(pos[2]) || isnan(pos[2]) ) {      
425 +      sprintf( painCave.errMsg,
426 +               "DumpWriter detected a numerical error writing the position"
427 +               " for object %d", index);      
428 +      painCave.isFatal = 1;
429 +      simError();
430 +    }
431 +
432      vel = integrableObject->getVel();          
433 +
434 +    if (isinf(vel[0]) || isnan(vel[0]) ||
435 +        isinf(vel[1]) || isnan(vel[1]) ||
436 +        isinf(vel[2]) || isnan(vel[2]) ) {      
437 +      sprintf( painCave.errMsg,
438 +               "DumpWriter detected a numerical error writing the velocity"
439 +               " for object %d", index);      
440 +      painCave.isFatal = 1;
441 +      simError();
442 +    }
443 +
444      sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
445              pos[0], pos[1], pos[2],
446              vel[0], vel[1], vel[2]);                    
# Line 323 | Line 451 | namespace oopse {
451        Quat4d q;
452        Vector3d ji;
453        q = integrableObject->getQ();
454 +
455 +      if (isinf(q[0]) || isnan(q[0]) ||
456 +          isinf(q[1]) || isnan(q[1]) ||
457 +          isinf(q[2]) || isnan(q[2]) ||
458 +          isinf(q[3]) || isnan(q[3]) ) {      
459 +        sprintf( painCave.errMsg,
460 +                 "DumpWriter detected a numerical error writing the quaternion"
461 +                 " for object %d", index);      
462 +        painCave.isFatal = 1;
463 +        simError();
464 +      }
465 +
466        ji = integrableObject->getJ();
467 +
468 +      if (isinf(ji[0]) || isnan(ji[0]) ||
469 +          isinf(ji[1]) || isnan(ji[1]) ||
470 +          isinf(ji[2]) || isnan(ji[2]) ) {      
471 +        sprintf( painCave.errMsg,
472 +                 "DumpWriter detected a numerical error writing the angular"
473 +                 " momentum for object %d", index);      
474 +        painCave.isFatal = 1;
475 +        simError();
476 +      }
477 +
478        sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
479                q[0], q[1], q[2], q[3],
480                ji[0], ji[1], ji[2]);
# Line 331 | Line 482 | namespace oopse {
482      }
483  
484      if (needForceVector_) {
485 <      type += "ft";
486 <      Vector3d frc;
487 <      Vector3d trq;
488 <      frc = integrableObject->getFrc();
489 <      trq = integrableObject->getTrq();
490 <              
491 <      sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e",
492 <              frc[0], frc[1], frc[2],
493 <              trq[0], trq[1], trq[2]);
485 >      type += "f";
486 >      Vector3d frc = integrableObject->getFrc();
487 >      if (isinf(frc[0]) || isnan(frc[0]) ||
488 >          isinf(frc[1]) || isnan(frc[1]) ||
489 >          isinf(frc[2]) || isnan(frc[2]) ) {      
490 >        sprintf( painCave.errMsg,
491 >                 "DumpWriter detected a numerical error writing the force"
492 >                 " for object %d", index);      
493 >        painCave.isFatal = 1;
494 >        simError();
495 >      }
496 >      sprintf(tempBuffer, " %13e %13e %13e",
497 >              frc[0], frc[1], frc[2]);
498        line += tempBuffer;
499 +      
500 +      if (integrableObject->isDirectional()) {
501 +        type += "t";
502 +        Vector3d trq = integrableObject->getTrq();        
503 +        if (isinf(trq[0]) || isnan(trq[0]) ||
504 +            isinf(trq[1]) || isnan(trq[1]) ||
505 +            isinf(trq[2]) || isnan(trq[2]) ) {      
506 +          sprintf( painCave.errMsg,
507 +                   "DumpWriter detected a numerical error writing the torque"
508 +                   " for object %d", index);      
509 +          painCave.isFatal = 1;
510 +          simError();
511 +        }        
512 +        sprintf(tempBuffer, " %13e %13e %13e",
513 +                trq[0], trq[1], trq[2]);
514 +        line += tempBuffer;
515 +      }      
516      }
517 <        
517 >
518      sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
519      return std::string(tempBuffer);
520    }
521  
522 +  std::string DumpWriter::prepareSiteLine(StuntDouble* integrableObject, int ioIndex, int siteIndex) {
523 +        
524 +
525 +    std::string id;
526 +    std::string type;
527 +    std::string line;
528 +    char tempBuffer[4096];
529 +
530 +    if (integrableObject->isRigidBody()) {
531 +      sprintf(tempBuffer, "%10d           ", ioIndex);
532 +      id = std::string(tempBuffer);
533 +    } else {
534 +      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
535 +      id = std::string(tempBuffer);
536 +    }
537 +              
538 +    if (needFlucQ_) {
539 +      type += "cw";
540 +      RealType fqPos = integrableObject->getFlucQPos();
541 +      if (isinf(fqPos) || isnan(fqPos) ) {      
542 +        sprintf( painCave.errMsg,
543 +                 "DumpWriter detected a numerical error writing the"
544 +                 " fluctuating charge for object %s", id.c_str());      
545 +        painCave.isFatal = 1;
546 +        simError();
547 +      }
548 +      sprintf(tempBuffer, " %13e ", fqPos);
549 +      line += tempBuffer;
550 +    
551 +      RealType fqVel = integrableObject->getFlucQVel();
552 +      if (isinf(fqVel) || isnan(fqVel) ) {      
553 +        sprintf( painCave.errMsg,
554 +                 "DumpWriter detected a numerical error writing the"
555 +                 " fluctuating charge velocity for object %s", id.c_str());      
556 +        painCave.isFatal = 1;
557 +        simError();
558 +      }
559 +      sprintf(tempBuffer, " %13e ", fqVel);
560 +      line += tempBuffer;
561 +
562 +      if (needForceVector_) {
563 +        type += "g";
564 +        RealType fqFrc = integrableObject->getFlucQFrc();        
565 +        if (isinf(fqFrc) || isnan(fqFrc) ) {      
566 +          sprintf( painCave.errMsg,
567 +                   "DumpWriter detected a numerical error writing the"
568 +                   " fluctuating charge force for object %s", id.c_str());      
569 +          painCave.isFatal = 1;
570 +          simError();
571 +        }
572 +        sprintf(tempBuffer, " %13e ", fqFrc);        
573 +        line += tempBuffer;
574 +      }
575 +    }
576 +
577 +    if (needElectricField_) {
578 +      type += "e";
579 +      Vector3d eField= integrableObject->getElectricField();
580 +      if (isinf(eField[0]) || isnan(eField[0]) ||
581 +          isinf(eField[1]) || isnan(eField[1]) ||
582 +          isinf(eField[2]) || isnan(eField[2]) ) {      
583 +        sprintf( painCave.errMsg,
584 +                 "DumpWriter detected a numerical error writing the electric"
585 +                 " field for object %s", id.c_str());      
586 +        painCave.isFatal = 1;
587 +        simError();
588 +      }
589 +      sprintf(tempBuffer, " %13e %13e %13e",
590 +              eField[0], eField[1], eField[2]);
591 +      line += tempBuffer;
592 +    }
593 +
594 +
595 +    if (needParticlePot_) {
596 +      type += "u";
597 +      RealType particlePot = integrableObject->getParticlePot();
598 +      if (isinf(particlePot) || isnan(particlePot)) {      
599 +        sprintf( painCave.errMsg,
600 +                 "DumpWriter detected a numerical error writing the particle "
601 +                 " potential for object %s", id.c_str());      
602 +        painCave.isFatal = 1;
603 +        simError();
604 +      }
605 +      sprintf(tempBuffer, " %13e", particlePot);
606 +      line += tempBuffer;
607 +    }
608 +    
609 +
610 +    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
611 +    return std::string(tempBuffer);
612 +  }
613 +
614    void DumpWriter::writeDump() {
615      writeFrame(*dumpFile_);
616    }
# Line 424 | Line 688 | namespace oopse {
688      newOStream = new std::ofstream(filename.c_str());
689   #endif
690      //write out MetaData first
691 <    (*newOStream) << "<OOPSE version=4>" << std::endl;
691 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
692      (*newOStream) << "  <MetaData>" << std::endl;
693      (*newOStream) << info_->getRawMetaData();
694      (*newOStream) << "  </MetaData>" << std::endl;
# Line 433 | Line 697 | namespace oopse {
697  
698    void DumpWriter::writeClosing(std::ostream& os) {
699  
700 <    os << "</OOPSE>\n";
700 >    os << "</OpenMD>\n";
701      os.flush();
702    }
703  
704 < }//end namespace oopse
704 > }//end namespace OpenMD

Comparing:
trunk/src/io/DumpWriter.cpp (property svn:keywords), Revision 1025 by gezelter, Wed Aug 30 20:33:44 2006 UTC vs.
branches/development/src/io/DumpWriter.cpp (property svn:keywords), Revision 1752 by gezelter, Sun Jun 10 14:05:02 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines