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

Comparing branches/development/src/io/DumpReader.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1794 by gezelter, Thu Sep 6 19:44:06 2012 UTC

# Line 36 | Line 36
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).                        
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
# Line 57 | Line 58
58   #include "utils/simError.h"
59   #include "utils/MemoryUtils.hpp"
60   #include "utils/StringTokenizer.hpp"
61 + #include "brains/Thermo.hpp"
62  
63   #ifdef IS_MPI
62 –
64   #include <mpi.h>
65 < #define TAKE_THIS_TAG_CHAR 0
65 < #define TAKE_THIS_TAG_INT 1
65 > #endif
66  
67 – #endif // is_mpi
67  
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 142 | Line 141 | namespace OpenMD {
141        prevPos = currPos;
142        bool foundOpenSnapshotTag = false;
143        bool foundClosedSnapshotTag = false;
144 +      bool foundOpenSiteDataTag = false;
145        while(inFile_->getline(buffer, bufferSize)) {
146          ++lineNo;
147          
# Line 227 | Line 227 | namespace OpenMD {
227        needVel_ = false;
228      }
229      
230 <    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
230 >    if (storageLayout & DataStorage::dslAmat ||
231 >        storageLayout & DataStorage::dslDipole ||
232 >        storageLayout & DataStorage::dslQuadrupole) {
233        needQuaternion_ = true;
234      } else {
235        needQuaternion_ = false;
# Line 243 | Line 245 | namespace OpenMD {
245  
246      if (needCOMprops_) {
247        Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
248 +      Thermo thermo(info_);
249        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 –    }
250  
251 +      if (needPos_ && needVel_) {
252 +        Vector3d comvel;
253 +        Vector3d comw;
254 +        thermo.getComAll(com, comvel);
255 +        comw = thermo.getAngularMomentum();
256 +      } else {
257 +        com = thermo.getCom();
258 +      }                    
259 +    }
260    }
261    
262    void DumpReader::readSet(int whichFrame) {    
# Line 324 | Line 324 | namespace OpenMD {
324  
325      inputStream.getline(buffer, bufferSize);
326      line = buffer;
327 <    if (line.find("</Snapshot>") == std::string::npos) {
328 <      sprintf(painCave.errMsg,
329 <              "DumpReader Error: can not find </Snapshot>\n");
330 <      painCave.isFatal = 1;
331 <      simError();
332 <    }        
333 <  
327 >
328 >    if (line.find("<SiteData>") != std::string::npos) {
329 >      //read SiteData
330 >      readSiteData(inputStream);        
331 >    } else {
332 >      if (line.find("</Snapshot>") == std::string::npos) {
333 >        sprintf(painCave.errMsg,
334 >                "DumpReader Error: can not find </Snapshot>\n");
335 >        painCave.isFatal = 1;
336 >        simError();
337 >      }        
338 >    }
339    }
340    
341    void DumpReader::parseDumpLine(const std::string& line) {
# Line 350 | Line 355 | namespace OpenMD {
355  
356      int index = tokenizer.nextTokenAsInt();
357  
358 <    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
358 >    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
359  
360 <    if (integrableObject == NULL) {
360 >    if (sd == NULL) {
361        return;
362      }
363      std::string type = tokenizer.nextToken();
# Line 372 | Line 377 | namespace OpenMD {
377        }
378      }
379      
380 <    if (integrableObject->isDirectional()) {
380 >    if (sd->isDirectional()) {
381        if (needQuaternion_) {
382          found = type.find("q");      
383          if (found == std::string::npos) {
# Line 395 | Line 400 | namespace OpenMD {
400              pos[1] = tokenizer.nextTokenAsDouble();
401              pos[2] = tokenizer.nextTokenAsDouble();
402              if (needPos_) {
403 <              integrableObject->setPos(pos);
403 >              sd->setPos(pos);
404              }            
405              break;
406          }
# Line 405 | Line 410 | namespace OpenMD {
410              vel[1] = tokenizer.nextTokenAsDouble();
411              vel[2] = tokenizer.nextTokenAsDouble();
412              if (needVel_) {
413 <              integrableObject->setVel(vel);
413 >              sd->setVel(vel);
414              }
415              break;
416          }
417  
418          case 'q' : {
419             Quat4d q;
420 <           if (integrableObject->isDirectional()) {
420 >           if (sd->isDirectional()) {
421                
422               q[0] = tokenizer.nextTokenAsDouble();
423               q[1] = tokenizer.nextTokenAsDouble();
# Line 431 | Line 436 | namespace OpenMD {
436                
437               q.normalize();
438               if (needQuaternion_) {            
439 <               integrableObject->setQ(q);
439 >               sd->setQ(q);
440               }              
441             }            
442             break;
443          }  
444          case 'j' : {
445            Vector3d ji;
446 <          if (integrableObject->isDirectional()) {
446 >          if (sd->isDirectional()) {
447               ji[0] = tokenizer.nextTokenAsDouble();
448               ji[1] = tokenizer.nextTokenAsDouble();
449               ji[2] = tokenizer.nextTokenAsDouble();
450               if (needAngMom_) {
451 <               integrableObject->setJ(ji);
451 >               sd->setJ(ji);
452               }
453            }
454            break;
# Line 454 | Line 459 | namespace OpenMD {
459            force[0] = tokenizer.nextTokenAsDouble();
460            force[1] = tokenizer.nextTokenAsDouble();
461            force[2] = tokenizer.nextTokenAsDouble();          
462 <          integrableObject->setFrc(force);
462 >          sd->setFrc(force);
463            break;
464          }
465          case 't' : {
# Line 463 | Line 468 | namespace OpenMD {
468             torque[0] = tokenizer.nextTokenAsDouble();
469             torque[1] = tokenizer.nextTokenAsDouble();
470             torque[2] = tokenizer.nextTokenAsDouble();          
471 <           integrableObject->setTrq(torque);          
471 >           sd->setTrq(torque);          
472             break;
473          }
474 +        case 'u' : {
475 +
476 +           RealType particlePot;
477 +           particlePot = tokenizer.nextTokenAsDouble();
478 +           sd->setParticlePot(particlePot);          
479 +           break;
480 +        }
481 +        case 'c' : {
482 +
483 +           RealType flucQPos;
484 +           flucQPos = tokenizer.nextTokenAsDouble();
485 +           sd->setFlucQPos(flucQPos);          
486 +           break;
487 +        }
488 +        case 'w' : {
489 +
490 +           RealType flucQVel;
491 +           flucQVel = tokenizer.nextTokenAsDouble();
492 +           sd->setFlucQVel(flucQVel);          
493 +           break;
494 +        }
495 +        case 'g' : {
496 +
497 +           RealType flucQFrc;
498 +           flucQFrc = tokenizer.nextTokenAsDouble();
499 +           sd->setFlucQFrc(flucQFrc);          
500 +           break;
501 +        }
502 +        case 'e' : {
503 +
504 +           Vector3d eField;
505 +           eField[0] = tokenizer.nextTokenAsDouble();
506 +           eField[1] = tokenizer.nextTokenAsDouble();
507 +           eField[2] = tokenizer.nextTokenAsDouble();          
508 +           sd->setElectricField(eField);          
509 +           break;
510 +        }
511          default: {
512                 sprintf(painCave.errMsg,
513                         "DumpReader Error: %s is an unrecognized type\n", type.c_str());
# Line 480 | Line 522 | namespace OpenMD {
522    }
523    
524  
525 <  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
525 >  void DumpReader::parseSiteLine(const std::string& line) {
526 >
527 >    StringTokenizer tokenizer(line);
528 >    int nTokens;
529 >    
530 >    nTokens = tokenizer.countTokens();
531 >    
532 >    if (nTokens < 2) {  
533 >      sprintf(painCave.errMsg,
534 >              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
535 >      painCave.isFatal = 1;
536 >      simError();
537 >    }
538  
539 +    /**
540 +     * The first token is the global integrable object index.
541 +     */
542 +
543 +    int index = tokenizer.nextTokenAsInt();
544 +    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
545 +    if (sd == NULL) {
546 +      return;
547 +    }
548 +
549 +    /**
550 +     * Test to see if the next token is an integer or not.  If not,
551 +     * we've got data on the integrable object itself.  If there is an
552 +     * integer, we're parsing data for a site on a rigid body.
553 +     */
554 +
555 +    std::string indexTest = tokenizer.peekNextToken();
556 +    std::istringstream i(indexTest);
557 +    int siteIndex;
558 +    if (i >> siteIndex) {
559 +      // chew up this token and parse as an int:
560 +      siteIndex = tokenizer.nextTokenAsInt();
561 +      RigidBody* rb = static_cast<RigidBody*>(sd);
562 +      sd = rb->getAtoms()[siteIndex];
563 +    }
564 +
565 +    /**
566 +     * The next token contains information on what follows.
567 +     */
568 +    std::string type = tokenizer.nextToken();
569 +    int size = type.size();
570 +    
571 +    for(int i = 0; i < size; ++i) {
572 +      switch(type[i]) {
573 +        
574 +      case 'u' : {
575 +        
576 +        RealType particlePot;
577 +        particlePot = tokenizer.nextTokenAsDouble();
578 +        sd->setParticlePot(particlePot);
579 +        break;
580 +      }
581 +      case 'c' : {
582 +        
583 +        RealType flucQPos;
584 +        flucQPos = tokenizer.nextTokenAsDouble();
585 +        sd->setFlucQPos(flucQPos);
586 +        break;
587 +      }
588 +      case 'w' : {
589 +        
590 +        RealType flucQVel;
591 +        flucQVel = tokenizer.nextTokenAsDouble();
592 +        sd->setFlucQVel(flucQVel);
593 +        break;
594 +      }
595 +      case 'g' : {
596 +        
597 +        RealType flucQFrc;
598 +        flucQFrc = tokenizer.nextTokenAsDouble();
599 +        sd->setFlucQFrc(flucQFrc);
600 +        break;
601 +      }
602 +      case 'e' : {
603 +        
604 +        Vector3d eField;
605 +        eField[0] = tokenizer.nextTokenAsDouble();
606 +        eField[1] = tokenizer.nextTokenAsDouble();
607 +        eField[2] = tokenizer.nextTokenAsDouble();  
608 +        sd->setElectricField(eField);          
609 +        break;
610 +      }
611 +      default: {
612 +        sprintf(painCave.errMsg,
613 +                "DumpReader Error: %s is an unrecognized type\n", type.c_str());
614 +        painCave.isFatal = 1;
615 +        simError();
616 +        break;  
617 +      }
618 +      }
619 +    }    
620 +  }
621 +  
622 +  
623 +  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
624 +    
625      inputStream.getline(buffer, bufferSize);
626      std::string line(buffer);
627      
# Line 504 | Line 644 | namespace OpenMD {
644    
645    }
646  
647 +  void  DumpReader::readSiteData(std::istream& inputStream) {
648 +
649 +    inputStream.getline(buffer, bufferSize);
650 +    std::string line(buffer);
651 +    
652 +    if (line.find("<SiteData>") == std::string::npos) {
653 +      // site data isn't required for a simulation, so skip
654 +      return;
655 +    }
656 +
657 +    while(inputStream.getline(buffer, bufferSize)) {
658 +      line = buffer;
659 +      
660 +      if(line.find("</SiteData>") != std::string::npos) {
661 +        break;
662 +      }
663 +
664 +      parseSiteLine(line);
665 +    }
666 +  
667 +  }
668 +
669    void DumpReader::readFrameProperties(std::istream& inputStream) {
670  
671      Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 549 | Line 711 | namespace OpenMD {
711          hmat(2, 2) = tokenizer.nextTokenAsDouble();
712          s->setHmat(hmat);      
713        } else if (propertyName == "Thermostat") {
714 <        RealType chi = tokenizer.nextTokenAsDouble();
715 <        RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
716 <        s->setChi(chi);
717 <        s->setIntegralOfChiDt(integralOfChiDt);        
714 >        pair<RealType, RealType> thermostat;
715 >        thermostat.first = tokenizer.nextTokenAsDouble();
716 >        thermostat.second = tokenizer.nextTokenAsDouble();
717 >        s->setThermostat(thermostat);
718       } else if (propertyName == "Barostat") {
719          Mat3x3d eta;
720          eta(0, 0) = tokenizer.nextTokenAsDouble();
# Line 564 | Line 726 | namespace OpenMD {
726          eta(2, 0) = tokenizer.nextTokenAsDouble();
727          eta(2, 1) = tokenizer.nextTokenAsDouble();
728          eta(2, 2) = tokenizer.nextTokenAsDouble();
729 <        s->setEta(eta);
729 >        s->setBarostat(eta);
730        } else {
731          sprintf(painCave.errMsg,
732                  "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());

Diff Legend

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