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 1871 by gezelter, Fri May 10 14:59:13 2013 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
# 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 +
145        while(inFile_->getline(buffer, bufferSize)) {
146          ++lineNo;
147          
# Line 202 | Line 202 | namespace OpenMD {
202   #ifdef IS_MPI
203      }
204      
205 <    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
205 >    MPI::COMM_WORLD.Bcast(&nframes_, 1, MPI::INT, 0);
206      
207   #endif // is_mpi
208      
# 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;
250 <      Vector3d comvel;
251 <      Vector3d comw;
252 <      if (needPos_ && needVel_){
253 <        info_->getComAll(com, comvel);
254 <        comw = info_->getAngularMomentum();
255 <      }else{
256 <        com = info_->getCom();
257 <        comvel = 0.0;
258 <        comw   = 0.0;
256 <      }
257 <      s->setCOMprops(com, comvel, comw);      
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      }
259 –
260    }
261    
262    void DumpReader::readSet(int whichFrame) {    
# Line 288 | Line 288 | namespace OpenMD {
288        }
289  
290        int sendBufferSize = sendBuffer.size();
291 <      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
292 <      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
291 >      MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode);    
292 >      MPI::COMM_WORLD.Bcast((void *)sendBuffer.c_str(), sendBufferSize,
293 >                            MPI::CHAR, masterNode);    
294        
295        sstream.str(sendBuffer);
296      } else {
297        int sendBufferSize;
298 <      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
298 >      MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode);
299        char * recvBuffer = new char[sendBufferSize+1];
300        assert(recvBuffer);
301        recvBuffer[sendBufferSize] = '\0';
302 <      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
302 >      MPI::COMM_WORLD.Bcast(recvBuffer, sendBufferSize, MPI::CHAR, masterNode);
303        sstream.str(recvBuffer);
304        delete [] recvBuffer;
305      }      
# Line 324 | Line 325 | namespace OpenMD {
325  
326      inputStream.getline(buffer, bufferSize);
327      line = buffer;
328 <    if (line.find("</Snapshot>") == std::string::npos) {
329 <      sprintf(painCave.errMsg,
330 <              "DumpReader Error: can not find </Snapshot>\n");
331 <      painCave.isFatal = 1;
332 <      simError();
333 <    }        
334 <  
328 >
329 >    if (line.find("<SiteData>") != std::string::npos) {
330 >      //read SiteData
331 >      readSiteData(inputStream);        
332 >    } else {
333 >      if (line.find("</Snapshot>") == std::string::npos) {
334 >        sprintf(painCave.errMsg,
335 >                "DumpReader Error: can not find </Snapshot>\n");
336 >        painCave.isFatal = 1;
337 >        simError();
338 >      }        
339 >    }
340    }
341    
342    void DumpReader::parseDumpLine(const std::string& line) {
# Line 350 | Line 356 | namespace OpenMD {
356  
357      int index = tokenizer.nextTokenAsInt();
358  
359 <    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
359 >    StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
360  
361 <    if (integrableObject == NULL) {
361 >    if (sd == NULL) {
362        return;
363      }
364      std::string type = tokenizer.nextToken();
# Line 372 | Line 378 | namespace OpenMD {
378        }
379      }
380      
381 <    if (integrableObject->isDirectional()) {
381 >    if (sd->isDirectional()) {
382        if (needQuaternion_) {
383          found = type.find("q");      
384          if (found == std::string::npos) {
# Line 395 | Line 401 | namespace OpenMD {
401              pos[1] = tokenizer.nextTokenAsDouble();
402              pos[2] = tokenizer.nextTokenAsDouble();
403              if (needPos_) {
404 <              integrableObject->setPos(pos);
404 >              sd->setPos(pos);
405              }            
406              break;
407          }
# Line 405 | Line 411 | namespace OpenMD {
411              vel[1] = tokenizer.nextTokenAsDouble();
412              vel[2] = tokenizer.nextTokenAsDouble();
413              if (needVel_) {
414 <              integrableObject->setVel(vel);
414 >              sd->setVel(vel);
415              }
416              break;
417          }
418  
419          case 'q' : {
420             Quat4d q;
421 <           if (integrableObject->isDirectional()) {
421 >           if (sd->isDirectional()) {
422                
423               q[0] = tokenizer.nextTokenAsDouble();
424               q[1] = tokenizer.nextTokenAsDouble();
# Line 431 | Line 437 | namespace OpenMD {
437                
438               q.normalize();
439               if (needQuaternion_) {            
440 <               integrableObject->setQ(q);
440 >               sd->setQ(q);
441               }              
442             }            
443             break;
444          }  
445          case 'j' : {
446            Vector3d ji;
447 <          if (integrableObject->isDirectional()) {
447 >          if (sd->isDirectional()) {
448               ji[0] = tokenizer.nextTokenAsDouble();
449               ji[1] = tokenizer.nextTokenAsDouble();
450               ji[2] = tokenizer.nextTokenAsDouble();
451               if (needAngMom_) {
452 <               integrableObject->setJ(ji);
452 >               sd->setJ(ji);
453               }
454            }
455            break;
# Line 454 | Line 460 | namespace OpenMD {
460            force[0] = tokenizer.nextTokenAsDouble();
461            force[1] = tokenizer.nextTokenAsDouble();
462            force[2] = tokenizer.nextTokenAsDouble();          
463 <          integrableObject->setFrc(force);
463 >          sd->setFrc(force);
464            break;
465          }
466          case 't' : {
# Line 463 | Line 469 | namespace OpenMD {
469             torque[0] = tokenizer.nextTokenAsDouble();
470             torque[1] = tokenizer.nextTokenAsDouble();
471             torque[2] = tokenizer.nextTokenAsDouble();          
472 <           integrableObject->setTrq(torque);          
472 >           sd->setTrq(torque);          
473             break;
474          }
475 +        case 'u' : {
476 +
477 +           RealType particlePot;
478 +           particlePot = tokenizer.nextTokenAsDouble();
479 +           sd->setParticlePot(particlePot);          
480 +           break;
481 +        }
482 +        case 'c' : {
483 +
484 +           RealType flucQPos;
485 +           flucQPos = tokenizer.nextTokenAsDouble();
486 +           sd->setFlucQPos(flucQPos);          
487 +           break;
488 +        }
489 +        case 'w' : {
490 +
491 +           RealType flucQVel;
492 +           flucQVel = tokenizer.nextTokenAsDouble();
493 +           sd->setFlucQVel(flucQVel);          
494 +           break;
495 +        }
496 +        case 'g' : {
497 +
498 +           RealType flucQFrc;
499 +           flucQFrc = tokenizer.nextTokenAsDouble();
500 +           sd->setFlucQFrc(flucQFrc);          
501 +           break;
502 +        }
503 +        case 'e' : {
504 +
505 +           Vector3d eField;
506 +           eField[0] = tokenizer.nextTokenAsDouble();
507 +           eField[1] = tokenizer.nextTokenAsDouble();
508 +           eField[2] = tokenizer.nextTokenAsDouble();          
509 +           sd->setElectricField(eField);          
510 +           break;
511 +        }
512          default: {
513                 sprintf(painCave.errMsg,
514                         "DumpReader Error: %s is an unrecognized type\n", type.c_str());
# Line 480 | Line 523 | namespace OpenMD {
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 +      RigidBody* rb = static_cast<RigidBody*>(sd);
563 +      sd = rb->getAtoms()[siteIndex];
564 +    }
565 +
566 +    /**
567 +     * The next token contains information on what follows.
568 +     */
569 +    std::string type = tokenizer.nextToken();
570 +    int size = type.size();
571 +    
572 +    for(int i = 0; i < size; ++i) {
573 +      switch(type[i]) {
574 +        
575 +      case 'u' : {
576 +        
577 +        RealType particlePot;
578 +        particlePot = tokenizer.nextTokenAsDouble();
579 +        sd->setParticlePot(particlePot);
580 +        break;
581 +      }
582 +      case 'c' : {
583 +        
584 +        RealType flucQPos;
585 +        flucQPos = tokenizer.nextTokenAsDouble();
586 +        sd->setFlucQPos(flucQPos);
587 +        break;
588 +      }
589 +      case 'w' : {
590 +        
591 +        RealType flucQVel;
592 +        flucQVel = tokenizer.nextTokenAsDouble();
593 +        sd->setFlucQVel(flucQVel);
594 +        break;
595 +      }
596 +      case 'g' : {
597 +        
598 +        RealType flucQFrc;
599 +        flucQFrc = tokenizer.nextTokenAsDouble();
600 +        sd->setFlucQFrc(flucQFrc);
601 +        break;
602 +      }
603 +      case 'e' : {
604 +        
605 +        Vector3d eField;
606 +        eField[0] = tokenizer.nextTokenAsDouble();
607 +        eField[1] = tokenizer.nextTokenAsDouble();
608 +        eField[2] = tokenizer.nextTokenAsDouble();  
609 +        sd->setElectricField(eField);          
610 +        break;
611 +      }
612 +      default: {
613 +        sprintf(painCave.errMsg,
614 +                "DumpReader Error: %s is an unrecognized type\n", type.c_str());
615 +        painCave.isFatal = 1;
616 +        simError();
617 +        break;  
618 +      }
619 +      }
620 +    }    
621 +  }
622 +  
623 +  
624 +  void  DumpReader::readStuntDoubles(std::istream& inputStream) {
625 +    
626      inputStream.getline(buffer, bufferSize);
627      std::string line(buffer);
628      
# Line 504 | Line 645 | namespace OpenMD {
645    
646    }
647  
648 +  void  DumpReader::readSiteData(std::istream& inputStream) {
649 +
650 +    std::string line(buffer);
651 +
652 +    // We already found the starting <SiteData> tag or we wouldn't be
653 +    // here, so just start parsing until we get to the ending
654 +    // </SiteData> tag:
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 549 | 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 564 | 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