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 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 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 243 | Line 242 | namespace OpenMD {
242  
243      if (needCOMprops_) {
244        Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
245 +      Thermo thermo(info_);
246        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);      
247 >
248 >      if (needPos_ && needVel_) {
249 >        Vector3d comvel;
250 >        Vector3d comw;
251 >        thermo.getComAll(com, comvel);
252 >        comw = thermo.getAngularMomentum();
253 >      } else {
254 >        com = thermo.getCom();
255 >      }                    
256      }
259 –
257    }
258    
259    void DumpReader::readSet(int whichFrame) {    
# Line 324 | Line 321 | namespace OpenMD {
321  
322      inputStream.getline(buffer, bufferSize);
323      line = buffer;
324 <    if (line.find("</Snapshot>") == std::string::npos) {
325 <      sprintf(painCave.errMsg,
326 <              "DumpReader Error: can not find </Snapshot>\n");
327 <      painCave.isFatal = 1;
328 <      simError();
329 <    }        
330 <  
324 >
325 >    if (line.find("<SiteData>") != std::string::npos) {
326 >      //read SiteData
327 >      readSiteData(inputStream);        
328 >    } else {
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 >    }
336    }
337    
338    void DumpReader::parseDumpLine(const std::string& line) {
# Line 341 | Line 343 | namespace OpenMD {
343      
344      nTokens = tokenizer.countTokens();
345      
346 <    if (nTokens < 2) {
346 >    if (nTokens < 2) {  
347        sprintf(painCave.errMsg,
348                "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
349        painCave.isFatal = 1;
# Line 350 | Line 352 | namespace OpenMD {
352  
353      int index = tokenizer.nextTokenAsInt();
354  
355 <    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
355 >    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
356  
357 <    if (integrableObject == NULL) {
357 >    if (sd == NULL) {
358        return;
359      }
360      std::string type = tokenizer.nextToken();
361      int size = type.size();
362  
363 +    size_t found;
364 +    
365 +    if (needPos_) {
366 +      found = type.find("p");      
367 +      if (found == std::string::npos) {
368 +        sprintf(painCave.errMsg,
369 +                "DumpReader Error: StuntDouble %d has no Position\n"
370 +                "\tField (\"p\") specified.\n%s\n", index,
371 +                line.c_str());  
372 +        painCave.isFatal = 1;
373 +        simError();
374 +      }
375 +    }
376 +    
377 +    if (sd->isDirectional()) {
378 +      if (needQuaternion_) {
379 +        found = type.find("q");      
380 +        if (found == std::string::npos) {
381 +          sprintf(painCave.errMsg,
382 +                  "DumpReader Error: Directional StuntDouble %d has no\n"
383 +                  "\tQuaternion Field (\"q\") specified.\n%s\n", index,
384 +                  line.c_str());  
385 +          painCave.isFatal = 1;
386 +          simError();
387 +        }
388 +      }      
389 +    }
390 +
391      for(int i = 0; i < size; ++i) {
392        switch(type[i]) {
393          
# Line 367 | Line 397 | namespace OpenMD {
397              pos[1] = tokenizer.nextTokenAsDouble();
398              pos[2] = tokenizer.nextTokenAsDouble();
399              if (needPos_) {
400 <              integrableObject->setPos(pos);
400 >              sd->setPos(pos);
401              }            
402              break;
403          }
# Line 377 | Line 407 | namespace OpenMD {
407              vel[1] = tokenizer.nextTokenAsDouble();
408              vel[2] = tokenizer.nextTokenAsDouble();
409              if (needVel_) {
410 <              integrableObject->setVel(vel);
410 >              sd->setVel(vel);
411              }
412              break;
413          }
414  
415          case 'q' : {
416             Quat4d q;
417 <           if (integrableObject->isDirectional()) {
417 >           if (sd->isDirectional()) {
418                
419               q[0] = tokenizer.nextTokenAsDouble();
420               q[1] = tokenizer.nextTokenAsDouble();
# Line 403 | Line 433 | namespace OpenMD {
433                
434               q.normalize();
435               if (needQuaternion_) {            
436 <               integrableObject->setQ(q);
436 >               sd->setQ(q);
437               }              
438             }            
439             break;
440          }  
441          case 'j' : {
442            Vector3d ji;
443 <          if (integrableObject->isDirectional()) {
443 >          if (sd->isDirectional()) {
444               ji[0] = tokenizer.nextTokenAsDouble();
445               ji[1] = tokenizer.nextTokenAsDouble();
446               ji[2] = tokenizer.nextTokenAsDouble();
447               if (needAngMom_) {
448 <               integrableObject->setJ(ji);
448 >               sd->setJ(ji);
449               }
450            }
451            break;
# Line 426 | Line 456 | namespace OpenMD {
456            force[0] = tokenizer.nextTokenAsDouble();
457            force[1] = tokenizer.nextTokenAsDouble();
458            force[2] = tokenizer.nextTokenAsDouble();          
459 <          integrableObject->setFrc(force);
459 >          sd->setFrc(force);
460            break;
461          }
462          case 't' : {
# Line 435 | Line 465 | namespace OpenMD {
465             torque[0] = tokenizer.nextTokenAsDouble();
466             torque[1] = tokenizer.nextTokenAsDouble();
467             torque[2] = tokenizer.nextTokenAsDouble();          
468 <           integrableObject->setTrq(torque);          
468 >           sd->setTrq(torque);          
469             break;
470          }
471 +        case 'u' : {
472 +
473 +           RealType particlePot;
474 +           particlePot = tokenizer.nextTokenAsDouble();
475 +           sd->setParticlePot(particlePot);          
476 +           break;
477 +        }
478 +        case 'c' : {
479 +
480 +           RealType flucQPos;
481 +           flucQPos = tokenizer.nextTokenAsDouble();
482 +           sd->setFlucQPos(flucQPos);          
483 +           break;
484 +        }
485 +        case 'w' : {
486 +
487 +           RealType flucQVel;
488 +           flucQVel = tokenizer.nextTokenAsDouble();
489 +           sd->setFlucQVel(flucQVel);          
490 +           break;
491 +        }
492 +        case 'g' : {
493 +
494 +           RealType flucQFrc;
495 +           flucQFrc = tokenizer.nextTokenAsDouble();
496 +           sd->setFlucQFrc(flucQFrc);          
497 +           break;
498 +        }
499 +        case 'e' : {
500 +
501 +           Vector3d eField;
502 +           eField[0] = tokenizer.nextTokenAsDouble();
503 +           eField[1] = tokenizer.nextTokenAsDouble();
504 +           eField[2] = tokenizer.nextTokenAsDouble();          
505 +           sd->setElectricField(eField);          
506 +           break;
507 +        }
508          default: {
509                 sprintf(painCave.errMsg,
510                         "DumpReader Error: %s is an unrecognized type\n", type.c_str());
# Line 452 | Line 519 | namespace OpenMD {
519    }
520    
521  
522 <  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
522 >  void DumpReader::parseSiteLine(const std::string& line) {
523  
524 +    StringTokenizer tokenizer(line);
525 +    int nTokens;
526 +    
527 +    nTokens = tokenizer.countTokens();
528 +    
529 +    if (nTokens < 2) {  
530 +      sprintf(painCave.errMsg,
531 +              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
532 +      painCave.isFatal = 1;
533 +      simError();
534 +    }
535 +
536 +    /**
537 +     * The first token is the global integrable object index.
538 +     */
539 +
540 +    int index = tokenizer.nextTokenAsInt();
541 +    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
542 +    if (sd == NULL) {
543 +      return;
544 +    }
545 +
546 +    /**
547 +     * Test to see if the next token is an integer or not.  If not,
548 +     * we've got data on the integrable object itself.  If there is an
549 +     * integer, we're parsing data for a site on a rigid body.
550 +     */
551 +
552 +    std::string indexTest = tokenizer.peekNextToken();
553 +    std::istringstream i(indexTest);
554 +    int siteIndex;
555 +    if (i >> siteIndex) {
556 +      // chew up this token and parse as an int:
557 +      siteIndex = tokenizer.nextTokenAsInt();
558 +      RigidBody* rb = static_cast<RigidBody*>(sd);
559 +      sd = rb->getAtoms()[siteIndex];
560 +    }
561 +
562 +    /**
563 +     * The next token contains information on what follows.
564 +     */
565 +    std::string type = tokenizer.nextToken();
566 +    int size = type.size();
567 +    
568 +    for(int i = 0; i < size; ++i) {
569 +      switch(type[i]) {
570 +        
571 +      case 'u' : {
572 +        
573 +        RealType particlePot;
574 +        particlePot = tokenizer.nextTokenAsDouble();
575 +        sd->setParticlePot(particlePot);
576 +        break;
577 +      }
578 +      case 'c' : {
579 +        
580 +        RealType flucQPos;
581 +        flucQPos = tokenizer.nextTokenAsDouble();
582 +        sd->setFlucQPos(flucQPos);
583 +        break;
584 +      }
585 +      case 'w' : {
586 +        
587 +        RealType flucQVel;
588 +        flucQVel = tokenizer.nextTokenAsDouble();
589 +        sd->setFlucQVel(flucQVel);
590 +        break;
591 +      }
592 +      case 'g' : {
593 +        
594 +        RealType flucQFrc;
595 +        flucQFrc = tokenizer.nextTokenAsDouble();
596 +        sd->setFlucQFrc(flucQFrc);
597 +        break;
598 +      }
599 +      case 'e' : {
600 +        
601 +        Vector3d eField;
602 +        eField[0] = tokenizer.nextTokenAsDouble();
603 +        eField[1] = tokenizer.nextTokenAsDouble();
604 +        eField[2] = tokenizer.nextTokenAsDouble();  
605 +        sd->setElectricField(eField);          
606 +        break;
607 +      }
608 +      default: {
609 +        sprintf(painCave.errMsg,
610 +                "DumpReader Error: %s is an unrecognized type\n", type.c_str());
611 +        painCave.isFatal = 1;
612 +        simError();
613 +        break;  
614 +      }
615 +      }
616 +    }    
617 +  }
618 +  
619 +  
620 +  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
621 +    
622      inputStream.getline(buffer, bufferSize);
623      std::string line(buffer);
624      
# Line 476 | Line 641 | namespace OpenMD {
641    
642    }
643  
644 +  void  DumpReader::readSiteData(std::istream& inputStream) {
645 +
646 +    inputStream.getline(buffer, bufferSize);
647 +    std::string line(buffer);
648 +    
649 +    if (line.find("<SiteData>") == std::string::npos) {
650 +      // site data isn't required for a simulation, so skip
651 +      return;
652 +    }
653 +
654 +    while(inputStream.getline(buffer, bufferSize)) {
655 +      line = buffer;
656 +      
657 +      if(line.find("</SiteData>") != std::string::npos) {
658 +        break;
659 +      }
660 +
661 +      parseSiteLine(line);
662 +    }
663 +  
664 +  }
665 +
666    void DumpReader::readFrameProperties(std::istream& inputStream) {
667  
668      Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 521 | Line 708 | namespace OpenMD {
708          hmat(2, 2) = tokenizer.nextTokenAsDouble();
709          s->setHmat(hmat);      
710        } else if (propertyName == "Thermostat") {
711 <        RealType chi = tokenizer.nextTokenAsDouble();
712 <        RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
713 <        s->setChi(chi);
714 <        s->setIntegralOfChiDt(integralOfChiDt);        
711 >        pair<RealType, RealType> thermostat;
712 >        thermostat.first = tokenizer.nextTokenAsDouble();
713 >        thermostat.second = tokenizer.nextTokenAsDouble();
714 >        s->setThermostat(thermostat);
715       } else if (propertyName == "Barostat") {
716          Mat3x3d eta;
717          eta(0, 0) = tokenizer.nextTokenAsDouble();
# Line 536 | Line 723 | namespace OpenMD {
723          eta(2, 0) = tokenizer.nextTokenAsDouble();
724          eta(2, 1) = tokenizer.nextTokenAsDouble();
725          eta(2, 2) = tokenizer.nextTokenAsDouble();
726 <        s->setEta(eta);
726 >        s->setBarostat(eta);
727        } else {
728          sprintf(painCave.errMsg,
729                  "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());

Diff Legend

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