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 1025 by gezelter, Wed Aug 30 20:33:44 2006 UTC vs.
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 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, 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 < namespace oopse {
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 oopse {
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 oopse {
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 oopse {
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 224 | Line 304 | namespace oopse {
304    void DumpWriter::writeFrame(std::ostream& os) {
305  
306   #ifdef IS_MPI
307 <    MPI_Status istatus;
307 >    MPI::Status istatus;
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 oopse {
319      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
320  
321      os << "    <StuntDoubles>\n";
322 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
323 <
324 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
325 <           integrableObject = mol->nextIntegrableObject(ii)) {  
326 <        os << prepareDumpLine(integrableObject);
327 <
322 >    for (mol = info_->beginMolecule(mi); mol != NULL;
323 >         mol = info_->nextMolecule(mi)) {
324 >      
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  
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
363 +
364 +    const int masterNode = 0;
365 +    int worldRank = MPI::COMM_WORLD.Get_rank();
366 +    int nProc = MPI::COMM_WORLD.Get_size();
367 +
368 +    if (worldRank == masterNode) {      
369 +      os << "  <Snapshot>\n";  
370 +      writeFrameProperties(os,
371 +                           info_->getSnapshotManager()->getCurrentSnapshot());
372 +      os << "    <StuntDoubles>\n";
373 +    }
374 +
375      //every node prepares the dump lines for integrable objects belong to itself
376      std::string buffer;
377 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
378 <      
379 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
380 <           integrableObject = mol->nextIntegrableObject(ii)) {  
381 <        buffer += prepareDumpLine(integrableObject);
377 >    for (mol = info_->beginMolecule(mi); mol != NULL;
378 >         mol = info_->nextMolecule(mi)) {
379 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;
380 >           sd = mol->nextIntegrableObject(ii)) {        
381 >        buffer += prepareDumpLine(sd);
382        }
383      }
384      
265    const int masterNode = 0;
266
385      if (worldRank == masterNode) {      
268      os << "  <Snapshot>\n";  
269      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
270      os << "    <StuntDoubles>\n";
271        
386        os << buffer;
387 <
274 <      int nProc;
275 <      MPI_Comm_size(MPI_COMM_WORLD, &nProc);
387 >      
388        for (int i = 1; i < nProc; ++i) {
389 +        // tell processor i to start sending us data:
390 +        MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
391  
392          // receive the length of the string buffer that was
393 <        // prepared by processor i
280 <
393 >        // prepared by processor i:        
394          int recvLength;
395 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
395 >        MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
396 >                             istatus);
397 >
398 >        // create a buffer to receive the data
399          char* recvBuffer = new char[recvLength];
400          if (recvBuffer == NULL) {
401          } else {
402 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
402 >          // receive the data:
403 >          MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
404 >                               MPI::ANY_TAG, istatus);
405 >          // send it to the file:
406            os << recvBuffer;
407 <          delete recvBuffer;
407 >          // get rid of the receive buffer:
408 >          delete [] recvBuffer;
409          }
410        }
291      os << "    </StuntDoubles>\n";
292      
293      os << "  </Snapshot>\n";
294      os.flush();
411      } else {
412        int sendBufferLength = buffer.size() + 1;
413 <      MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
414 <      MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
413 >      int myturn = 0;
414 >      for (int i = 1; i < nProc; ++i){
415 >        // wait for the master node to call our number:
416 >        MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
417 >        if (myturn == worldRank){
418 >          // send the length of our buffer:
419 >          MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
420 >
421 >          // send our buffer:
422 >          MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
423 >                               MPI::CHAR, masterNode, 0);
424 >
425 >        }
426 >      }
427      }
428 +    
429 +    if (worldRank == masterNode) {      
430 +      os << "    </StuntDoubles>\n";
431 +    }
432  
433 < #endif // is_mpi
433 >    if (doSiteData_) {
434 >      if (worldRank == masterNode) {
435 >        os << "    <SiteData>\n";
436 >      }
437 >      buffer.clear();
438 >      for (mol = info_->beginMolecule(mi); mol != NULL;
439 >           mol = info_->nextMolecule(mi)) {
440 >              
441 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
442 >             sd = mol->nextIntegrableObject(ii)) {      
443 >          
444 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
445 >          // do one for the IO itself
446 >          buffer += prepareSiteLine(sd, ioIndex, 0);
447  
448 <  }
448 >          if (sd->isRigidBody()) {
449 >            
450 >            RigidBody* rb = static_cast<RigidBody*>(sd);
451 >            int siteIndex = 0;
452 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
453 >                 atom = rb->nextAtom(ai)) {                                            
454 >              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
455 >              siteIndex++;
456 >            }
457 >          }
458 >        }
459 >      }
460  
461 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
461 >      if (worldRank == masterNode) {    
462 >        os << buffer;
463 >        
464 >        for (int i = 1; i < nProc; ++i) {
465 >          
466 >          // tell processor i to start sending us data:
467 >          MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
468 >          
469 >          // receive the length of the string buffer that was
470 >          // prepared by processor i:        
471 >          int recvLength;
472 >          MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
473 >                               istatus);
474 >          
475 >          // create a buffer to receive the data
476 >          char* recvBuffer = new char[recvLength];
477 >          if (recvBuffer == NULL) {
478 >          } else {
479 >            // receive the data:
480 >            MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
481 >                                 MPI::ANY_TAG, istatus);
482 >            // send it to the file:
483 >            os << recvBuffer;
484 >            // get rid of the receive buffer:
485 >            delete [] recvBuffer;
486 >          }
487 >        }      
488 >      } else {
489 >        int sendBufferLength = buffer.size() + 1;
490 >        int myturn = 0;
491 >        for (int i = 1; i < nProc; ++i){
492 >          // wait for the master node to call our number:
493 >          MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
494 >          if (myturn == worldRank){
495 >            // send the length of our buffer:
496 >            MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
497 >            // send our buffer:
498 >            MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
499 >                                 MPI::CHAR, masterNode, 0);
500 >          }
501 >        }
502 >      }
503 >      
504 >      if (worldRank == masterNode) {    
505 >        os << "    </SiteData>\n";
506 >      }
507 >    }
508 >    
509 >    if (worldRank == masterNode) {
510 >      os << "  </Snapshot>\n";
511 >      os.flush();
512 >    }
513 >    
514 > #endif // is_mpi
515 >    
516 >  }
517 >
518 >  std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
519          
520 <    int index = integrableObject->getGlobalIntegrableObjectIndex();
520 >    int index = sd->getGlobalIntegrableObjectIndex();
521      std::string type("pv");
522      std::string line;
523      char tempBuffer[4096];
524  
525      Vector3d pos;
526      Vector3d vel;
527 <    pos = integrableObject->getPos();
528 <    vel = integrableObject->getVel();          
527 >    pos = sd->getPos();
528 >
529 >    if (isinf(pos[0]) || isnan(pos[0]) ||
530 >        isinf(pos[1]) || isnan(pos[1]) ||
531 >        isinf(pos[2]) || isnan(pos[2]) ) {      
532 >      sprintf( painCave.errMsg,
533 >               "DumpWriter detected a numerical error writing the position"
534 >               " for object %d", index);      
535 >      painCave.isFatal = 1;
536 >      simError();
537 >    }
538 >
539 >    vel = sd->getVel();        
540 >
541 >    if (isinf(vel[0]) || isnan(vel[0]) ||
542 >        isinf(vel[1]) || isnan(vel[1]) ||
543 >        isinf(vel[2]) || isnan(vel[2]) ) {      
544 >      sprintf( painCave.errMsg,
545 >               "DumpWriter detected a numerical error writing the velocity"
546 >               " for object %d", index);      
547 >      painCave.isFatal = 1;
548 >      simError();
549 >    }
550 >
551      sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
552              pos[0], pos[1], pos[2],
553              vel[0], vel[1], vel[2]);                    
554      line += tempBuffer;
555  
556 <    if (integrableObject->isDirectional()) {
556 >    if (sd->isDirectional()) {
557        type += "qj";
558        Quat4d q;
559        Vector3d ji;
560 <      q = integrableObject->getQ();
561 <      ji = integrableObject->getJ();
560 >      q = sd->getQ();
561 >
562 >      if (isinf(q[0]) || isnan(q[0]) ||
563 >          isinf(q[1]) || isnan(q[1]) ||
564 >          isinf(q[2]) || isnan(q[2]) ||
565 >          isinf(q[3]) || isnan(q[3]) ) {      
566 >        sprintf( painCave.errMsg,
567 >                 "DumpWriter detected a numerical error writing the quaternion"
568 >                 " for object %d", index);      
569 >        painCave.isFatal = 1;
570 >        simError();
571 >      }
572 >
573 >      ji = sd->getJ();
574 >
575 >      if (isinf(ji[0]) || isnan(ji[0]) ||
576 >          isinf(ji[1]) || isnan(ji[1]) ||
577 >          isinf(ji[2]) || isnan(ji[2]) ) {      
578 >        sprintf( painCave.errMsg,
579 >                 "DumpWriter detected a numerical error writing the angular"
580 >                 " momentum for object %d", index);      
581 >        painCave.isFatal = 1;
582 >        simError();
583 >      }
584 >
585        sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
586                q[0], q[1], q[2], q[3],
587                ji[0], ji[1], ji[2]);
# Line 331 | Line 589 | namespace oopse {
589      }
590  
591      if (needForceVector_) {
592 <      type += "ft";
593 <      Vector3d frc;
594 <      Vector3d trq;
595 <      frc = integrableObject->getFrc();
596 <      trq = integrableObject->getTrq();
597 <              
598 <      sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e",
599 <              frc[0], frc[1], frc[2],
600 <              trq[0], trq[1], trq[2]);
592 >      type += "f";
593 >      Vector3d frc = sd->getFrc();
594 >      if (isinf(frc[0]) || isnan(frc[0]) ||
595 >          isinf(frc[1]) || isnan(frc[1]) ||
596 >          isinf(frc[2]) || isnan(frc[2]) ) {      
597 >        sprintf( painCave.errMsg,
598 >                 "DumpWriter detected a numerical error writing the force"
599 >                 " for object %d", index);      
600 >        painCave.isFatal = 1;
601 >        simError();
602 >      }
603 >      sprintf(tempBuffer, " %13e %13e %13e",
604 >              frc[0], frc[1], frc[2]);
605        line += tempBuffer;
606 +      
607 +      if (sd->isDirectional()) {
608 +        type += "t";
609 +        Vector3d trq = sd->getTrq();        
610 +        if (isinf(trq[0]) || isnan(trq[0]) ||
611 +            isinf(trq[1]) || isnan(trq[1]) ||
612 +            isinf(trq[2]) || isnan(trq[2]) ) {      
613 +          sprintf( painCave.errMsg,
614 +                   "DumpWriter detected a numerical error writing the torque"
615 +                   " for object %d", index);      
616 +          painCave.isFatal = 1;
617 +          simError();
618 +        }        
619 +        sprintf(tempBuffer, " %13e %13e %13e",
620 +                trq[0], trq[1], trq[2]);
621 +        line += tempBuffer;
622 +      }      
623      }
624 <        
624 >
625      sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
626      return std::string(tempBuffer);
627    }
628  
629 +  std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
630 +    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
631 +
632 +    std::string id;
633 +    std::string type;
634 +    std::string line;
635 +    char tempBuffer[4096];
636 +
637 +    if (sd->isRigidBody()) {
638 +      sprintf(tempBuffer, "%10d           ", ioIndex);
639 +      id = std::string(tempBuffer);
640 +    } else {
641 +      sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
642 +      id = std::string(tempBuffer);
643 +    }
644 +              
645 +    if (needFlucQ_) {
646 +      if (storageLayout & DataStorage::dslFlucQPosition) {
647 +        type += "c";
648 +        RealType fqPos = sd->getFlucQPos();
649 +        if (isinf(fqPos) || isnan(fqPos) ) {      
650 +          sprintf( painCave.errMsg,
651 +                   "DumpWriter detected a numerical error writing the"
652 +                   " fluctuating charge for object %s", id.c_str());      
653 +          painCave.isFatal = 1;
654 +          simError();
655 +        }
656 +        sprintf(tempBuffer, " %13e ", fqPos);
657 +        line += tempBuffer;
658 +      }
659 +
660 +      if (storageLayout & DataStorage::dslFlucQVelocity) {
661 +        type += "w";    
662 +        RealType fqVel = sd->getFlucQVel();
663 +        if (isinf(fqVel) || isnan(fqVel) ) {      
664 +          sprintf( painCave.errMsg,
665 +                   "DumpWriter detected a numerical error writing the"
666 +                   " fluctuating charge velocity for object %s", id.c_str());      
667 +          painCave.isFatal = 1;
668 +          simError();
669 +        }
670 +        sprintf(tempBuffer, " %13e ", fqVel);
671 +        line += tempBuffer;
672 +      }
673 +
674 +      if (needForceVector_) {
675 +        if (storageLayout & DataStorage::dslFlucQForce) {          
676 +          type += "g";
677 +          RealType fqFrc = sd->getFlucQFrc();        
678 +          if (isinf(fqFrc) || isnan(fqFrc) ) {      
679 +            sprintf( painCave.errMsg,
680 +                     "DumpWriter detected a numerical error writing the"
681 +                     " fluctuating charge force for object %s", id.c_str());      
682 +            painCave.isFatal = 1;
683 +            simError();
684 +          }
685 +          sprintf(tempBuffer, " %13e ", fqFrc);        
686 +          line += tempBuffer;
687 +        }
688 +      }
689 +    }
690 +    
691 +    if (needElectricField_) {
692 +      if (storageLayout & DataStorage::dslElectricField) {
693 +        type += "e";
694 +        Vector3d eField= sd->getElectricField();
695 +        if (isinf(eField[0]) || isnan(eField[0]) ||
696 +            isinf(eField[1]) || isnan(eField[1]) ||
697 +            isinf(eField[2]) || isnan(eField[2]) ) {      
698 +          sprintf( painCave.errMsg,
699 +                   "DumpWriter detected a numerical error writing the electric"
700 +                   " field for object %s", id.c_str());      
701 +          painCave.isFatal = 1;
702 +          simError();
703 +        }
704 +        sprintf(tempBuffer, " %13e %13e %13e",
705 +                eField[0], eField[1], eField[2]);
706 +        line += tempBuffer;
707 +      }
708 +    }
709 +
710 +
711 +    if (needParticlePot_) {
712 +      if (storageLayout & DataStorage::dslParticlePot) {
713 +        type += "u";
714 +        RealType particlePot = sd->getParticlePot();
715 +        if (isinf(particlePot) || isnan(particlePot)) {      
716 +          sprintf( painCave.errMsg,
717 +                   "DumpWriter detected a numerical error writing the particle "
718 +                   " potential for object %s", id.c_str());      
719 +          painCave.isFatal = 1;
720 +          simError();
721 +        }
722 +        sprintf(tempBuffer, " %13e", particlePot);
723 +        line += tempBuffer;
724 +      }
725 +    }
726 +  
727 +    sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
728 +    return std::string(tempBuffer);
729 +  }
730 +
731    void DumpWriter::writeDump() {
732      writeFrame(*dumpFile_);
733    }
# Line 357 | Line 738 | namespace oopse {
738   #ifdef IS_MPI
739      if (worldRank == 0) {
740   #endif // is_mpi
741 <
741 >      
742        eorStream = createOStream(eorFilename_);
743  
744   #ifdef IS_MPI
745      }
746 < #endif // is_mpi    
747 <
746 > #endif
747 >    
748      writeFrame(*eorStream);
749 <
749 >      
750   #ifdef IS_MPI
751      if (worldRank == 0) {
752 < #endif // is_mpi
752 > #endif
753 >      
754        writeClosing(*eorStream);
755        delete eorStream;
756 +      
757   #ifdef IS_MPI
758      }
759   #endif // is_mpi  
# Line 414 | Line 797 | namespace oopse {
797    std::ostream* DumpWriter::createOStream(const std::string& filename) {
798  
799      std::ostream* newOStream;
800 < #ifdef HAVE_LIBZ
800 > #ifdef HAVE_ZLIB
801      if (needCompression_) {
802        newOStream = new ogzstream(filename.c_str());
803      } else {
# Line 424 | Line 807 | namespace oopse {
807      newOStream = new std::ofstream(filename.c_str());
808   #endif
809      //write out MetaData first
810 <    (*newOStream) << "<OOPSE version=4>" << std::endl;
810 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
811      (*newOStream) << "  <MetaData>" << std::endl;
812      (*newOStream) << info_->getRawMetaData();
813      (*newOStream) << "  </MetaData>" << std::endl;
# Line 433 | Line 816 | namespace oopse {
816  
817    void DumpWriter::writeClosing(std::ostream& os) {
818  
819 <    os << "</OOPSE>\n";
819 >    os << "</OpenMD>\n";
820      os.flush();
821    }
822  
823 < }//end namespace oopse
823 > }//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.
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines