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 1217 by xsun, Wed Jan 23 21:23:02 2008 UTC vs.
branches/development/src/io/DumpWriter.cpp (file contents), Revision 1875 by gezelter, Fri May 17 14:41:42 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 "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 < namespace oopse {
61 > using namespace std;
62 > namespace OpenMD {
63  
64    DumpWriter::DumpWriter(SimInfo* info)
65      : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
66  
67      Globals* simParams = info->getSimParams();
68 <    needCompression_ = simParams->getCompressDumpFile();
69 <    needForceVector_ = simParams->getOutputForceVector();
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 95 | Line 114 | namespace oopse {
114      Globals* simParams = info->getSimParams();
115      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
116  
117 <    needCompression_ = simParams->getCompressDumpFile();
118 <    needForceVector_ = simParams->getOutputForceVector();
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 134 | Line 163 | namespace oopse {
163      Globals* simParams = info->getSimParams();
164      eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";    
165      
166 <    needCompression_ = simParams->getCompressDumpFile();
167 <    needForceVector_ = simParams->getOutputForceVector();
168 <    
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 194 | Line 232 | namespace oopse {
232      os << "    <FrameData>\n";
233  
234      RealType currentTime = s->getTime();
235 +
236 +    if (isinf(currentTime) || isnan(currentTime)) {      
237 +      sprintf( painCave.errMsg,
238 +               "DumpWriter detected a numerical error writing the time");      
239 +      painCave.isFatal = 1;
240 +      simError();
241 +    }
242 +    
243      sprintf(buffer, "        Time: %.10g\n", currentTime);
244      os << buffer;
245  
246      Mat3x3d hmat;
247      hmat = s->getHmat();
248 +
249 +    for (unsigned int i = 0; i < 3; i++) {
250 +      for (unsigned int j = 0; j < 3; j++) {
251 +        if (isinf(hmat(i,j)) || isnan(hmat(i,j))) {      
252 +          sprintf( painCave.errMsg,
253 +                   "DumpWriter detected a numerical error writing the box");
254 +          painCave.isFatal = 1;
255 +          simError();
256 +        }        
257 +      }
258 +    }
259 +    
260      sprintf(buffer, "        Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
261              hmat(0, 0), hmat(1, 0), hmat(2, 0),
262              hmat(0, 1), hmat(1, 1), hmat(2, 1),
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 <    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", chi, 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", 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++) {
284 >        if (isinf(eta(i,j)) || isnan(eta(i,j))) {      
285 >          sprintf( painCave.errMsg,
286 >                   "DumpWriter detected a numerical error writing the barostat");
287 >          painCave.isFatal = 1;
288 >          simError();
289 >        }        
290 >      }
291 >    }
292 >
293      sprintf(buffer, "    Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
294              eta(0, 0), eta(1, 0), eta(2, 0),
295              eta(0, 1), eta(1, 1), eta(2, 1),
# Line 224 | Line 302 | namespace oopse {
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";
# Line 238 | Line 317 | namespace oopse {
317      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
318  
319      os << "    <StuntDoubles>\n";
320 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
321 <
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
256    //every node prepares the dump lines for integrable objects belong to itself
257    std::string buffer;
258    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
361  
260
261      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
262           integrableObject = mol->nextIntegrableObject(ii)) {  
263          buffer += prepareDumpLine(integrableObject);
264      }
265    }
266    
362      const int masterNode = 0;
363 +    int worldRank = MPI::COMM_WORLD.Get_rank();
364 +    int nProc = MPI::COMM_WORLD.Get_size();
365  
366      if (worldRank == masterNode) {      
367        os << "  <Snapshot>\n";  
368 <      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
368 >      writeFrameProperties(os,
369 >                           info_->getSnapshotManager()->getCurrentSnapshot());
370        os << "    <StuntDoubles>\n";
371 <        
274 <      os << buffer;
371 >    }
372  
373 <      int nProc;
374 <      MPI_Comm_size(MPI_COMM_WORLD, &nProc);
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 >    
383 >    if (worldRank == masterNode) {      
384 >      os << buffer;
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
282 <
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 <          delete recvBuffer;
405 >          // get rid of the receive buffer:
406 >          delete [] recvBuffer;
407          }
408        }
293      os << "    </StuntDoubles>\n";
294      
295      os << "  </Snapshot>\n";
296      os.flush();
409      } else {
410        int sendBufferLength = buffer.size() + 1;
411 <      MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
412 <      MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
411 >      int myturn = 0;
412 >      for (int i = 1; i < nProc; ++i){
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 >          // 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 <  }
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 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
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* 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();
525 <    vel = integrableObject->getVel();          
524 >    pos = sd->getPos();
525 >
526 >    if (isinf(pos[0]) || isnan(pos[0]) ||
527 >        isinf(pos[1]) || isnan(pos[1]) ||
528 >        isinf(pos[2]) || isnan(pos[2]) ) {      
529 >      sprintf( painCave.errMsg,
530 >               "DumpWriter detected a numerical error writing the position"
531 >               " for object %d", index);      
532 >      painCave.isFatal = 1;
533 >      simError();
534 >    }
535 >
536 >    vel = sd->getVel();        
537 >
538 >    if (isinf(vel[0]) || isnan(vel[0]) ||
539 >        isinf(vel[1]) || isnan(vel[1]) ||
540 >        isinf(vel[2]) || isnan(vel[2]) ) {      
541 >      sprintf( painCave.errMsg,
542 >               "DumpWriter detected a numerical error writing the velocity"
543 >               " for object %d", index);      
544 >      painCave.isFatal = 1;
545 >      simError();
546 >    }
547 >
548      sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
549              pos[0], pos[1], pos[2],
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();
558 <      ji = integrableObject->getJ();
557 >      q = sd->getQ();
558 >
559 >      if (isinf(q[0]) || isnan(q[0]) ||
560 >          isinf(q[1]) || isnan(q[1]) ||
561 >          isinf(q[2]) || isnan(q[2]) ||
562 >          isinf(q[3]) || isnan(q[3]) ) {      
563 >        sprintf( painCave.errMsg,
564 >                 "DumpWriter detected a numerical error writing the quaternion"
565 >                 " for object %d", index);      
566 >        painCave.isFatal = 1;
567 >        simError();
568 >      }
569 >
570 >      ji = sd->getJ();
571 >
572 >      if (isinf(ji[0]) || isnan(ji[0]) ||
573 >          isinf(ji[1]) || isnan(ji[1]) ||
574 >          isinf(ji[2]) || isnan(ji[2]) ) {      
575 >        sprintf( painCave.errMsg,
576 >                 "DumpWriter detected a numerical error writing the angular"
577 >                 " momentum for object %d", index);      
578 >        painCave.isFatal = 1;
579 >        simError();
580 >      }
581 >
582        sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
583                q[0], q[1], q[2], q[3],
584                ji[0], ji[1], ji[2]);
# Line 333 | Line 586 | namespace oopse {
586      }
587  
588      if (needForceVector_) {
589 <      type += "ft";
590 <      Vector3d frc;
591 <      Vector3d trq;
592 <      frc = integrableObject->getFrc();
593 <      trq = integrableObject->getTrq();
594 <              
595 <      sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e",
596 <              frc[0], frc[1], frc[2],
597 <              trq[0], trq[1], trq[2]);
589 >      type += "f";
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]) ) {      
594 >        sprintf( painCave.errMsg,
595 >                 "DumpWriter detected a numerical error writing the force"
596 >                 " for object %d", index);      
597 >        painCave.isFatal = 1;
598 >        simError();
599 >      }
600 >      sprintf(tempBuffer, " %13e %13e %13e",
601 >              frc[0], frc[1], frc[2]);
602        line += tempBuffer;
603 +      
604 +      if (sd->isDirectional()) {
605 +        type += "t";
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]) ) {      
610 +          sprintf( painCave.errMsg,
611 +                   "DumpWriter detected a numerical error writing the torque"
612 +                   " for object %d", index);      
613 +          painCave.isFatal = 1;
614 +          simError();
615 +        }        
616 +        sprintf(tempBuffer, " %13e %13e %13e",
617 +                trq[0], trq[1], trq[2]);
618 +        line += tempBuffer;
619 +      }      
620      }
621 <        
621 >
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 +              
642 +    if (needFlucQ_) {
643 +      type += "cw";
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 %s", id.c_str());      
649 +        painCave.isFatal = 1;
650 +        simError();
651 +      }
652 +      sprintf(tempBuffer, " %13e ", fqPos);
653 +      line += tempBuffer;
654 +    
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 %s", id.c_str());      
660 +        painCave.isFatal = 1;
661 +        simError();
662 +      }
663 +      sprintf(tempBuffer, " %13e ", fqVel);
664 +      line += tempBuffer;
665 +
666 +      if (needForceVector_) {
667 +        type += "g";
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 %s", id.c_str());      
673 +          painCave.isFatal = 1;
674 +          simError();
675 +        }
676 +        sprintf(tempBuffer, " %13e ", fqFrc);        
677 +        line += tempBuffer;
678 +      }
679 +    }
680 +
681 +    if (needElectricField_) {
682 +      type += "e";
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 %s", id.c_str());      
690 +        painCave.isFatal = 1;
691 +        simError();
692 +      }
693 +      sprintf(tempBuffer, " %13e %13e %13e",
694 +              eField[0], eField[1], eField[2]);
695 +      line += tempBuffer;
696 +    }
697 +
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 +
718    void DumpWriter::writeDump() {
719      writeFrame(*dumpFile_);
720    }
# Line 359 | Line 725 | namespace oopse {
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      }
367 #endif // is_mpi    
368
369    writeFrame(*eorStream);
370
371 #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 416 | Line 780 | namespace oopse {
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 {
# Line 426 | Line 790 | namespace oopse {
790      newOStream = new std::ofstream(filename.c_str());
791   #endif
792      //write out MetaData first
793 <    (*newOStream) << "<OOPSE version=4>" << std::endl;
793 >    (*newOStream) << "<OpenMD version=2>" << std::endl;
794      (*newOStream) << "  <MetaData>" << std::endl;
795      (*newOStream) << info_->getRawMetaData();
796      (*newOStream) << "  </MetaData>" << std::endl;
# Line 435 | Line 799 | namespace oopse {
799  
800    void DumpWriter::writeClosing(std::ostream& os) {
801  
802 <    os << "</OOPSE>\n";
802 >    os << "</OpenMD>\n";
803      os.flush();
804    }
805  
806 < }//end namespace oopse
806 > }//end namespace OpenMD

Comparing:
trunk/src/io/DumpWriter.cpp (property svn:keywords), Revision 1217 by xsun, Wed Jan 23 21:23:02 2008 UTC vs.
branches/development/src/io/DumpWriter.cpp (property svn:keywords), Revision 1875 by gezelter, Fri May 17 14:41:42 2013 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines