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

Comparing trunk/src/io/DumpReader.cpp (file contents):
Revision 1610 by gezelter, Fri Aug 12 14:37:25 2011 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
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   #define _LARGEFILE_SOURCE64
44   #define _FILE_OFFSET_BITS 64
45  
46 + #ifdef IS_MPI
47 + #include <mpi.h>
48 + #endif
49 +
50   #include <sys/types.h>
51   #include <sys/stat.h>
52  
# Line 57 | Line 62
62   #include "utils/simError.h"
63   #include "utils/MemoryUtils.hpp"
64   #include "utils/StringTokenizer.hpp"
65 + #include "brains/Thermo.hpp"
66  
61 – #ifdef IS_MPI
67  
63 – #include <mpi.h>
64 – #define TAKE_THIS_TAG_CHAR 0
65 – #define TAKE_THIS_TAG_INT 1
66 –
67 – #endif // is_mpi
68 –
69 –
68   namespace OpenMD {
69    
70    DumpReader::DumpReader(SimInfo* info, const std::string& filename)
# Line 77 | Line 75 | namespace OpenMD {
75      if (worldRank == 0) {
76   #endif
77        
78 <      inFile_ = new std::ifstream(filename_.c_str());
78 >      inFile_ = new std::ifstream(filename_.c_str(),  
79 >                                  ifstream::in | ifstream::binary);
80        
81        if (inFile_->fail()) {
82          sprintf(painCave.errMsg,
# Line 129 | Line 128 | namespace OpenMD {
128    }
129    
130    void DumpReader::scanFile(void) {
131 <    int lineNo = 0;
131 >
132      std::streampos prevPos;
133      std::streampos  currPos;
134      
# Line 142 | Line 141 | namespace OpenMD {
141        prevPos = currPos;
142        bool foundOpenSnapshotTag = false;
143        bool foundClosedSnapshotTag = false;
144 +
145 +      int lineNo = 0;
146        while(inFile_->getline(buffer, bufferSize)) {
147          ++lineNo;
148          
# Line 227 | Line 228 | namespace OpenMD {
228        needVel_ = false;
229      }
230      
231 <    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
231 >    if (storageLayout & DataStorage::dslAmat ||
232 >        storageLayout & DataStorage::dslDipole ||
233 >        storageLayout & DataStorage::dslQuadrupole) {
234        needQuaternion_ = true;
235      } else {
236        needQuaternion_ = false;
# Line 243 | Line 246 | namespace OpenMD {
246  
247      if (needCOMprops_) {
248        Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
249 +      Thermo thermo(info_);
250        Vector3d com;
247 –      Vector3d comvel;
248 –      Vector3d comw;
249 –      if (needPos_ && needVel_){
250 –        info_->getComAll(com, comvel);
251 –        comw = info_->getAngularMomentum();
252 –      }else{
253 –        com = info_->getCom();
254 –        comvel = 0.0;
255 –        comw   = 0.0;
256 –      }
257 –      s->setCOMprops(com, comvel, comw);      
258 –    }
251  
252 +      if (needPos_ && needVel_) {
253 +        Vector3d comvel;
254 +        Vector3d comw;
255 +        thermo.getComAll(com, comvel);
256 +        comw = thermo.getAngularMomentum();
257 +      } else {
258 +        com = thermo.getCom();
259 +      }                    
260 +    }
261    }
262    
263    void DumpReader::readSet(int whichFrame) {    
# Line 289 | Line 290 | namespace OpenMD {
290  
291        int sendBufferSize = sendBuffer.size();
292        MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
293 <      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
293 >      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize,
294 >                MPI_CHAR, masterNode, MPI_COMM_WORLD);    
295        
296        sstream.str(sendBuffer);
297      } else {
298        int sendBufferSize;
299 <      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
299 >      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
300        char * recvBuffer = new char[sendBufferSize+1];
301        assert(recvBuffer);
302        recvBuffer[sendBufferSize] = '\0';
303 <      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
303 >      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
304        sstream.str(recvBuffer);
305        delete [] recvBuffer;
306      }      
# Line 324 | Line 326 | namespace OpenMD {
326  
327      inputStream.getline(buffer, bufferSize);
328      line = buffer;
329 <    if (line.find("</Snapshot>") == std::string::npos) {
330 <      sprintf(painCave.errMsg,
331 <              "DumpReader Error: can not find </Snapshot>\n");
332 <      painCave.isFatal = 1;
333 <      simError();
334 <    }        
335 <  
329 >
330 >    if (line.find("<SiteData>") != std::string::npos) {
331 >      //read SiteData
332 >      readSiteData(inputStream);        
333 >    } else {
334 >      if (line.find("</Snapshot>") == std::string::npos) {
335 >        sprintf(painCave.errMsg,
336 >                "DumpReader Error: can not find </Snapshot>\n");
337 >        painCave.isFatal = 1;
338 >        simError();
339 >      }        
340 >    }
341    }
342    
343    void DumpReader::parseDumpLine(const std::string& line) {
# Line 350 | Line 357 | namespace OpenMD {
357  
358      int index = tokenizer.nextTokenAsInt();
359  
360 <    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
360 >    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
361  
362 <    if (integrableObject == NULL) {
362 >    if (sd == NULL) {
363        return;
364      }
365      std::string type = tokenizer.nextToken();
# Line 372 | Line 379 | namespace OpenMD {
379        }
380      }
381      
382 <    if (integrableObject->isDirectional()) {
382 >    if (sd->isDirectional()) {
383        if (needQuaternion_) {
384          found = type.find("q");      
385          if (found == std::string::npos) {
# Line 395 | Line 402 | namespace OpenMD {
402              pos[1] = tokenizer.nextTokenAsDouble();
403              pos[2] = tokenizer.nextTokenAsDouble();
404              if (needPos_) {
405 <              integrableObject->setPos(pos);
405 >              sd->setPos(pos);
406              }            
407              break;
408          }
# Line 405 | Line 412 | namespace OpenMD {
412              vel[1] = tokenizer.nextTokenAsDouble();
413              vel[2] = tokenizer.nextTokenAsDouble();
414              if (needVel_) {
415 <              integrableObject->setVel(vel);
415 >              sd->setVel(vel);
416              }
417              break;
418          }
419  
420          case 'q' : {
421             Quat4d q;
422 <           if (integrableObject->isDirectional()) {
422 >           if (sd->isDirectional()) {
423                
424               q[0] = tokenizer.nextTokenAsDouble();
425               q[1] = tokenizer.nextTokenAsDouble();
# Line 431 | Line 438 | namespace OpenMD {
438                
439               q.normalize();
440               if (needQuaternion_) {            
441 <               integrableObject->setQ(q);
441 >               sd->setQ(q);
442               }              
443             }            
444             break;
445          }  
446          case 'j' : {
447            Vector3d ji;
448 <          if (integrableObject->isDirectional()) {
448 >          if (sd->isDirectional()) {
449               ji[0] = tokenizer.nextTokenAsDouble();
450               ji[1] = tokenizer.nextTokenAsDouble();
451               ji[2] = tokenizer.nextTokenAsDouble();
452               if (needAngMom_) {
453 <               integrableObject->setJ(ji);
453 >               sd->setJ(ji);
454               }
455            }
456            break;
# Line 454 | Line 461 | namespace OpenMD {
461            force[0] = tokenizer.nextTokenAsDouble();
462            force[1] = tokenizer.nextTokenAsDouble();
463            force[2] = tokenizer.nextTokenAsDouble();          
464 <          integrableObject->setFrc(force);
464 >          sd->setFrc(force);
465            break;
466          }
467          case 't' : {
# Line 463 | Line 470 | namespace OpenMD {
470             torque[0] = tokenizer.nextTokenAsDouble();
471             torque[1] = tokenizer.nextTokenAsDouble();
472             torque[2] = tokenizer.nextTokenAsDouble();          
473 <           integrableObject->setTrq(torque);          
473 >           sd->setTrq(torque);          
474             break;
475          }
476          case 'u' : {
477  
478             RealType particlePot;
479             particlePot = tokenizer.nextTokenAsDouble();
480 <           integrableObject->setParticlePot(particlePot);          
480 >           sd->setParticlePot(particlePot);          
481             break;
482          }
483 +        case 'c' : {
484 +
485 +           RealType flucQPos;
486 +           flucQPos = tokenizer.nextTokenAsDouble();
487 +           sd->setFlucQPos(flucQPos);          
488 +           break;
489 +        }
490 +        case 'w' : {
491 +
492 +           RealType flucQVel;
493 +           flucQVel = tokenizer.nextTokenAsDouble();
494 +           sd->setFlucQVel(flucQVel);          
495 +           break;
496 +        }
497 +        case 'g' : {
498 +
499 +           RealType flucQFrc;
500 +           flucQFrc = tokenizer.nextTokenAsDouble();
501 +           sd->setFlucQFrc(flucQFrc);          
502 +           break;
503 +        }
504 +        case 'e' : {
505 +
506 +           Vector3d eField;
507 +           eField[0] = tokenizer.nextTokenAsDouble();
508 +           eField[1] = tokenizer.nextTokenAsDouble();
509 +           eField[2] = tokenizer.nextTokenAsDouble();          
510 +           sd->setElectricField(eField);          
511 +           break;
512 +        }
513          default: {
514                 sprintf(painCave.errMsg,
515                         "DumpReader Error: %s is an unrecognized type\n", type.c_str());
# Line 483 | Line 520 | namespace OpenMD {
520  
521        }
522      }
486 –    
523    }
524    
525  
526 <  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
526 >  void DumpReader::parseSiteLine(const std::string& line) {
527  
528 +    StringTokenizer tokenizer(line);
529 +    int nTokens;
530 +        
531 +    nTokens = tokenizer.countTokens();
532 +    
533 +    if (nTokens < 2) {  
534 +      sprintf(painCave.errMsg,
535 +              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
536 +      painCave.isFatal = 1;
537 +      simError();
538 +    }
539 +
540 +    /**
541 +     * The first token is the global integrable object index.
542 +     */
543 +
544 +    int index = tokenizer.nextTokenAsInt();
545 +    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
546 +    if (sd == NULL) {
547 +      return;
548 +    }
549 +
550 +    /**
551 +     * Test to see if the next token is an integer or not.  If not,
552 +     * we've got data on the integrable object itself.  If there is an
553 +     * integer, we're parsing data for a site on a rigid body.
554 +     */
555 +
556 +    std::string indexTest = tokenizer.peekNextToken();
557 +    std::istringstream i(indexTest);
558 +    int siteIndex;
559 +    if (i >> siteIndex) {
560 +      // chew up this token and parse as an int:
561 +      siteIndex = tokenizer.nextTokenAsInt();
562 +
563 +      if (sd->isRigidBody()) {
564 +        RigidBody* rb = static_cast<RigidBody*>(sd);
565 +        sd = rb->getAtoms()[siteIndex];
566 +      }
567 +    }
568 +
569 +    /**
570 +     * The next token contains information on what follows.
571 +     */
572 +    std::string type = tokenizer.nextToken();
573 +    int size = type.size();
574 +    
575 +    for(int i = 0; i < size; ++i) {
576 +      switch(type[i]) {
577 +        
578 +      case 'u' : {
579 +        
580 +        RealType particlePot;
581 +        particlePot = tokenizer.nextTokenAsDouble();
582 +        sd->setParticlePot(particlePot);
583 +        break;
584 +      }
585 +      case 'c' : {
586 +        
587 +        RealType flucQPos;
588 +        flucQPos = tokenizer.nextTokenAsDouble();
589 +        sd->setFlucQPos(flucQPos);
590 +        break;
591 +      }
592 +      case 'w' : {
593 +        
594 +        RealType flucQVel;
595 +        flucQVel = tokenizer.nextTokenAsDouble();
596 +        sd->setFlucQVel(flucQVel);
597 +        break;
598 +      }
599 +      case 'g' : {
600 +        
601 +        RealType flucQFrc;
602 +        flucQFrc = tokenizer.nextTokenAsDouble();
603 +        sd->setFlucQFrc(flucQFrc);
604 +        break;
605 +      }
606 +      case 'e' : {
607 +        
608 +        Vector3d eField;
609 +        eField[0] = tokenizer.nextTokenAsDouble();
610 +        eField[1] = tokenizer.nextTokenAsDouble();
611 +        eField[2] = tokenizer.nextTokenAsDouble();  
612 +        sd->setElectricField(eField);          
613 +        break;
614 +      }
615 +      default: {
616 +        sprintf(painCave.errMsg,
617 +                "DumpReader Error: %s is an unrecognized type\n", type.c_str());
618 +        painCave.isFatal = 1;
619 +        simError();
620 +        break;  
621 +      }
622 +      }
623 +    }    
624 +  }
625 +  
626 +  
627 +  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
628 +    
629      inputStream.getline(buffer, bufferSize);
630      std::string line(buffer);
631      
# Line 511 | Line 648 | namespace OpenMD {
648    
649    }
650  
651 +  void  DumpReader::readSiteData(std::istream& inputStream) {
652 +
653 +    std::string line(buffer);
654 +
655 +    // We already found the starting <SiteData> tag or we wouldn't be
656 +    // here, so just start parsing until we get to the ending
657 +    // </SiteData> tag:
658 +    
659 +    while(inputStream.getline(buffer, bufferSize)) {
660 +      line = buffer;
661 +      
662 +      if(line.find("</SiteData>") != std::string::npos) {
663 +        break;
664 +      }
665 +
666 +      parseSiteLine(line);
667 +    }
668 +  
669 +  }
670 +
671    void DumpReader::readFrameProperties(std::istream& inputStream) {
672  
673      Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 556 | Line 713 | namespace OpenMD {
713          hmat(2, 2) = tokenizer.nextTokenAsDouble();
714          s->setHmat(hmat);      
715        } else if (propertyName == "Thermostat") {
716 <        RealType chi = tokenizer.nextTokenAsDouble();
717 <        RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
718 <        s->setChi(chi);
719 <        s->setIntegralOfChiDt(integralOfChiDt);        
716 >        pair<RealType, RealType> thermostat;
717 >        thermostat.first = tokenizer.nextTokenAsDouble();
718 >        thermostat.second = tokenizer.nextTokenAsDouble();
719 >        s->setThermostat(thermostat);
720       } else if (propertyName == "Barostat") {
721          Mat3x3d eta;
722          eta(0, 0) = tokenizer.nextTokenAsDouble();
# Line 571 | Line 728 | namespace OpenMD {
728          eta(2, 0) = tokenizer.nextTokenAsDouble();
729          eta(2, 1) = tokenizer.nextTokenAsDouble();
730          eta(2, 2) = tokenizer.nextTokenAsDouble();
731 <        s->setEta(eta);
731 >        s->setBarostat(eta);
732        } else {
733          sprintf(painCave.errMsg,
734                  "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());

Diff Legend

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