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

Diff Legend

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