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

Comparing trunk/src/io/RestReader.cpp (file contents):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1407 by cli2, Wed Jan 20 16:04:40 2010 UTC

# Line 1 | Line 1
1   /*
2 < * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
2 > * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
5   * non-exclusive, royalty free, license to use, modify and
# Line 38 | Line 38
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39   * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41 <  
42 < #include "io/DumpReader.hpp"
43 < #include "io/RestReader.hpp"
44 < #include "primitives/Molecule.hpp"
41 >
42 >
43 > #include <sys/types.h>
44 > #include <sys/stat.h>
45 >
46 > #include <math.h>
47 > #include <string>
48 > #include <sstream>
49 > #include <iostream>
50 > #include <stdio.h>
51 > #include <stdlib.h>
52 > #include <string.h>
53 >
54 > #include "io/RestReader.hpp"
55 > #include "primitives/Molecule.hpp"
56 > #include "utils/simError.h"
57 > #include "utils/MemoryUtils.hpp"
58 > #include "utils/StringTokenizer.hpp"
59   #include "restraints/ObjectRestraint.hpp"
60 < #include "restraints/MolecularRestraint.hpp"
60 > #include "restraints/MolecularRestraint.hpp"
61 >
62 > #ifdef IS_MPI
63 >
64 > #include <mpi.h>
65 > #endif
66  
67   namespace OpenMD {
68 +
69 +  void RestReader::scanFile(){
70 +    int lineNo = 0;
71 +    std::streampos prevPos;
72 +    std::streampos  currPos;
73 +    
74 + #ifdef IS_MPI
75 +    
76 +    if (worldRank == 0) {
77 + #endif // is_mpi
78 +
79 +      inFile_->clear();
80 +      currPos = inFile_->tellg();
81 +      prevPos = currPos;
82        
83 <  void RestReader::readReferenceStructure() {
83 >      bool foundOpenSnapshotTag = false;
84 >      
85 >      while(!foundOpenSnapshotTag && inFile_->getline(buffer, bufferSize)) {
86 >        ++lineNo;
87 >        
88 >        std::string line = buffer;
89 >        currPos = inFile_->tellg();
90 >        if (line.find("<Snapshot>")!= std::string::npos) {
91 >          foundOpenSnapshotTag = true;
92 >          framePos_ = prevPos;
93 >        }
94 >        prevPos = currPos;
95 >      }
96 >      
97 > #ifdef IS_MPI
98 >    }
99 >    MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD);
100 > #endif // is_mpi
101 >  }
102  
52 –    // some of this comes directly from DumpReader.
103  
104 <    if (!isScanned_)
105 <      scanFile();
106 <        
107 <    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
108 <    
109 <    if (storageLayout & DataStorage::dslPosition) {
110 <      needPos_ = true;
104 >  void RestReader::readSet(){
105 >    std::string line;
106 >    
107 > #ifndef IS_MPI
108 >    
109 >    inFile_->clear();
110 >    inFile_->seekg(framePos_);
111 >    
112 >    std::istream& inputStream = *inFile_;
113 > #else
114 >    
115 >    int masterNode = 0;
116 >    std::stringstream sstream;
117 >    if (worldRank == masterNode) {
118 >      std::string sendBuffer;
119 >      
120 >      inFile_->clear();  
121 >      inFile_->seekg(framePos_);
122 >      
123 >      while (inFile_->getline(buffer, bufferSize)) {
124 >        
125 >        line = buffer;
126 >        sendBuffer += line;
127 >        sendBuffer += '\n';
128 >        if (line.find("</Snapshot>") != std::string::npos) {
129 >          break;
130 >        }        
131 >      }
132 >      
133 >      int sendBufferSize = sendBuffer.size();
134 >      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
135 >      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
136 >      
137 >      sstream.str(sendBuffer);
138      } else {
139 <      needPos_ = false;
139 >      int sendBufferSize;
140 >      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);    
141 >      char * recvBuffer = new char[sendBufferSize+1];
142 >      assert(recvBuffer);
143 >      recvBuffer[sendBufferSize] = '\0';
144 >      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);    
145 >      sstream.str(recvBuffer);
146 >      delete [] recvBuffer;
147 >    }      
148 >    
149 >    std::istream& inputStream = sstream;  
150 > #endif
151 >    
152 >    inputStream.getline(buffer, bufferSize);
153 >
154 >    line = buffer;
155 >    if (line.find("<Snapshot>") == std::string::npos) {
156 >      sprintf(painCave.errMsg,
157 >              "RestReader Error: can not find <Snapshot>\n");
158 >      painCave.isFatal = 1;
159 >      simError();
160      }
161 <    
162 <    needVel_ = false;
163 <    
164 <    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
165 <      needQuaternion_ = true;
166 <    } else {
167 <      needQuaternion_ = false;
168 <    }
161 >    
162 >    //read frameData
163 >    readFrameProperties(inputStream);
164 >    
165 >    //read StuntDoubles
166 >    readStuntDoubles(inputStream);
167 >    
168 >    inputStream.getline(buffer, bufferSize);
169 >    line = buffer;
170 >    if (line.find("</Snapshot>") == std::string::npos) {
171 >      sprintf(painCave.errMsg,
172 >              "RestReader Error: can not find </Snapshot>\n");
173 >      painCave.isFatal = 1;
174 >      simError();
175 >    }            
176 >  }
177 >    
178 >  void RestReader::readReferenceStructure() {
179  
73 –    needAngMom_ = false;
74 –
180      // We need temporary storage to keep track of all StuntDouble positions
181      // in case some of the restraints are molecular (i.e. if they use
182      // multiple SD positions to determine restrained orientations or positions:
183  
184      all_pos_.clear();
185 <    all_pos_.resize(info_->getNGlobalIntegrableObjects() );
186 <
82 <
185 >    all_pos_.resize(info_->getNGlobalIntegrableObjects()) ;
186 >    
187      // Restraint files are just standard dump files, but with the reference
188      // structure stored in the first frame (frame 0).
189      // RestReader overloads readSet and explicitly handles all of the
190      // ObjectRestraints in that method:
191  
192 <    readSet(0);
193 <    
192 >    scanFile();
193 >
194 >
195 >    readSet();
196 >
197 >
198      // all ObjectRestraints have been handled, now we have to worry about
199      // molecular restraints:
200  
# Line 94 | Line 202 | namespace OpenMD {
202      Molecule::IntegrableObjectIterator j;
203      Molecule* mol;
204      StuntDouble* sd;
205 <
205 >    
206      // no need to worry about parallel molecules, as molecules are not
207      // split across processor boundaries.  Just loop over all molecules
208      // we know about:
209 <
209 >    
210      for (mol = info_->beginMolecule(i); mol != NULL;
211           mol = info_->nextMolecule(i)) {          
212 <
212 >      
213        // is this molecule restrained?    
214        GenericData* data = mol->getPropertyByName("Restraint");
215        
216        if (data != NULL) {
217 <
217 >        
218          // make sure we can reinterpret the generic data as restraint data:
219 <
219 >        
220          RestraintData* restData= dynamic_cast<RestraintData*>(data);        
221 <
221 >        
222          if (restData != NULL) {
223 <
223 >          
224            // make sure we can reinterpet the restraint data as a
225            // pointer to a MolecularRestraint:
226 <
226 >          
227            MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
228 <
228 >          
229            if (mRest != NULL) {          
230 <
230 >            
231              // now we need to pack the stunt doubles for the reference
232              // structure:
233 <
233 >            
234              std::vector<Vector3d> ref;
235              int count = 0;
236              RealType mass, mTot;
# Line 140 | Line 248 | namespace OpenMD {
248                // doubles from the *globally* sorted array of
249                // positions:
250                
251 <              ref.push_back( all_pos_[sd->getGlobalIndex()] );
252 <              
145 <              mass = sd->getMass();
146 <              
251 >              ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] );
252 >              mass = sd->getMass();              
253                COM = COM + mass * ref[count];
254                mTot = mTot + mass;
255                count = count + 1;
256              }
151 –            
257              COM /= mTot;
153 –
258              mRest->setReferenceStructure(ref, COM);        
155 –
259            }
260          }
261        }
# Line 162 | Line 265 | namespace OpenMD {
265    
266    
267    void RestReader::parseDumpLine(const std::string& line) {        
268 +
269      StringTokenizer tokenizer(line);
270      int nTokens;
271      
# Line 169 | Line 273 | namespace OpenMD {
273      
274      if (nTokens < 2) {
275        sprintf(painCave.errMsg,
276 <              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
276 >              "RestReader Error: Not enough Tokens.\n%s\n", line.c_str());
277        painCave.isFatal = 1;
278        simError();
279      }
280  
281      int index = tokenizer.nextTokenAsInt();
178 –
179 –    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
282  
283 <    if (integrableObject == NULL) {
284 <      return;
285 <    }  
286 <
287 <    std::string type = tokenizer.nextToken();
288 <    int size = type.size();
289 <    Vector3d pos;
290 <    Vector3d vel;
291 <    Quat4d q;
292 <    Vector3d ji;
293 <    Vector3d force;
294 <    Vector3d torque;
295 <          
296 <    for(int i = 0; i < size; ++i) {
297 <      switch(type[i]) {
283 >  
284 >    if (index < 1116){
285 >      std::string type = tokenizer.nextToken();
286 >      
287 >      Vector3d pos;
288 >      
289 >      pos[0] = tokenizer.nextTokenAsDouble();
290 >      pos[1] = tokenizer.nextTokenAsDouble();
291 >      pos[2] = tokenizer.nextTokenAsDouble();
292 >      
293 >      all_pos_[index] = pos;
294 >    }else{
295 >      
296 >      bool exist = false;
297 >      int indexCount = 0;
298 >      
299 >      while ( (!exist) && (indexCount < stuntDoubleIndex_.size())){
300 >        if (index == stuntDoubleIndex_[indexCount])
301 >          exist = true;
302 >        indexCount++;
303 >      }
304 >      
305 >      StuntDouble* integrableObject;
306 >      
307 >      if (exist){
308          
309 <        case 'p': {
310 <            pos[0] = tokenizer.nextTokenAsDouble();
311 <            pos[1] = tokenizer.nextTokenAsDouble();
312 <            pos[2] = tokenizer.nextTokenAsDouble();
313 <            break;
314 <        }
315 <        case 'v' : {
316 <            vel[0] = tokenizer.nextTokenAsDouble();
317 <            vel[1] = tokenizer.nextTokenAsDouble();
318 <            vel[2] = tokenizer.nextTokenAsDouble();
319 <            break;
320 <        }
321 <
322 <        case 'q' : {
323 <           if (integrableObject->isDirectional()) {
324 <              
325 <             q[0] = tokenizer.nextTokenAsDouble();
326 <             q[1] = tokenizer.nextTokenAsDouble();
327 <             q[2] = tokenizer.nextTokenAsDouble();
328 <             q[3] = tokenizer.nextTokenAsDouble();
329 <              
330 <             RealType qlen = q.length();
331 <             if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
332 <                
333 <               sprintf(painCave.errMsg,
334 <                       "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
335 <               painCave.isFatal = 1;
336 <               simError();
337 <                
226 <             }  
227 <              
228 <             q.normalize();
229 <           }            
230 <           break;
309 >        integrableObject = info_->getIOIndexToIntegrableObject(index);
310 >        
311 >        int compareindex = integrableObject->getGlobalIntegrableObjectIndex();
312 >        
313 >        if (integrableObject == NULL) {
314 >          return;
315 >        }  
316 >        
317 >        std::string type = tokenizer.nextToken();
318 >        
319 >        Vector3d pos;
320 >        
321 >        pos[0] = tokenizer.nextTokenAsDouble();
322 >        pos[1] = tokenizer.nextTokenAsDouble();
323 >        pos[2] = tokenizer.nextTokenAsDouble();
324 >        
325 >        Vector3d vel;
326 >        vel[0] = tokenizer.nextTokenAsDouble();
327 >        vel[1] = tokenizer.nextTokenAsDouble();
328 >        vel[2] = tokenizer.nextTokenAsDouble();
329 >        
330 >        
331 >        Quat4d q;
332 >        if (integrableObject->isDirectional()) {
333 >          
334 >          q[0] = tokenizer.nextTokenAsDouble();
335 >          q[1] = tokenizer.nextTokenAsDouble();
336 >          q[2] = tokenizer.nextTokenAsDouble();
337 >          q[3] = tokenizer.nextTokenAsDouble();
338          }  
339 <        case 'j' : {
340 <          if (integrableObject->isDirectional()) {
341 <             ji[0] = tokenizer.nextTokenAsDouble();
342 <             ji[1] = tokenizer.nextTokenAsDouble();
343 <             ji[2] = tokenizer.nextTokenAsDouble();
339 >        // keep the position in case we need it for a molecular restraint:
340 >        
341 >        all_pos_[index] = pos;
342 >        
343 >        // is this io restrained?
344 >        GenericData* data = integrableObject->getPropertyByName("Restraint");
345 >        ObjectRestraint* oRest;
346 >        
347 >        if (data != NULL) {
348 >          // make sure we can reinterpret the generic data as restraint data:
349 >          RestraintData* restData= dynamic_cast<RestraintData*>(data);        
350 >          if (restData != NULL) {
351 >            // make sure we can reinterpet the restraint data as a pointer to
352 >            // an ObjectRestraint:
353 >            oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
354 >            if (oRest != NULL) {          
355 >              if (integrableObject->isDirectional()) {
356 >                oRest->setReferenceStructure(pos, q.toRotationMatrix3());
357 >              } else {                          
358 >                oRest->setReferenceStructure(pos);
359 >              }
360 >            }
361            }
238 –          break;
239 –        }  
240 –        case 'f': {
241 –          force[0] = tokenizer.nextTokenAsDouble();
242 –          force[1] = tokenizer.nextTokenAsDouble();
243 –          force[2] = tokenizer.nextTokenAsDouble();          
244 –
245 –          break;
362          }
247 –        case 't' : {
248 –           torque[0] = tokenizer.nextTokenAsDouble();
249 –           torque[1] = tokenizer.nextTokenAsDouble();
250 –           torque[2] = tokenizer.nextTokenAsDouble();          
251 –
252 –           break;
253 –        }
254 –        default: {
255 –               sprintf(painCave.errMsg,
256 –                       "DumpReader Error: %s is an unrecognized type\n", type.c_str());
257 –               painCave.isFatal = 1;
258 –               simError();
259 –          break;  
260 –        }
261 –          
363        }
364      }
365 <    
366 <    // keep the position in case we need it for a molecular restraint:
365 >  }
366 >  
367 >  void  RestReader::readStuntDoubles(std::istream& inputStream) {
368      
369 <    all_pos_[index] = pos;
370 <
371 <    // is this io restrained?    
372 <    GenericData* data = integrableObject->getPropertyByName("Restraint");
373 <    ObjectRestraint* oRest;
374 <
375 <    if (data != NULL) {
376 <      // make sure we can reinterpret the generic data as restraint data:
377 <      RestraintData* restData= dynamic_cast<RestraintData*>(data);        
378 <      if (restData != NULL) {
379 <        // make sure we can reinterpet the restraint data as a pointer to
380 <        // an ObjectRestraint:
381 <        oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
382 <        if (oRest != NULL) {          
383 <          if (integrableObject->isDirectional()) {
282 <            oRest->setReferenceStructure(pos, q.toRotationMatrix3());
283 <          } else {                          
284 <            oRest->setReferenceStructure(pos);
285 <          }
286 <        }
369 >    inputStream.getline(buffer, bufferSize);
370 >    std::string line(buffer);
371 >    
372 >    if (line.find("<StuntDoubles>") == std::string::npos) {
373 >      sprintf(painCave.errMsg,
374 >              "RestReader Error: Missing <StuntDoubles>\n");
375 >      painCave.isFatal = 1;
376 >      simError();
377 >    }
378 >    
379 >    while(inputStream.getline(buffer, bufferSize)) {
380 >      line = buffer;
381 >      
382 >      if(line.find("</StuntDoubles>") != std::string::npos) {
383 >        break;
384        }
385 +      
386 +      parseDumpLine(line);
387      }
388 <
388 >    
389    }
291 –  
390  
391 <
391 >  
392    void RestReader::readFrameProperties(std::istream& inputStream) {
393      inputStream.getline(buffer, bufferSize);
394      std::string line(buffer);
395  
396      if (line.find("<FrameData>") == std::string::npos) {
397        sprintf(painCave.errMsg,
398 <              "DumpReader Error: Missing <FrameData>\n");
398 >              "RestReader Error: Missing <FrameData>\n");
399        painCave.isFatal = 1;
400        simError();
401      }
# Line 314 | Line 412 | namespace OpenMD {
412        }
413        
414      }
415 <
415 >    
416    }
417  
418    

Diff Legend

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