ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/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 1408 by gezelter, Fri Jan 22 21:26:29 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"the
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 >    readSet();
195 >
196 >
197      // all ObjectRestraints have been handled, now we have to worry about
198      // molecular restraints:
199  
# Line 94 | Line 201 | namespace OpenMD {
201      Molecule::IntegrableObjectIterator j;
202      Molecule* mol;
203      StuntDouble* sd;
204 <
204 >    
205      // no need to worry about parallel molecules, as molecules are not
206      // split across processor boundaries.  Just loop over all molecules
207      // we know about:
208 <
208 >    
209      for (mol = info_->beginMolecule(i); mol != NULL;
210           mol = info_->nextMolecule(i)) {          
211 <
211 >      
212        // is this molecule restrained?    
213        GenericData* data = mol->getPropertyByName("Restraint");
214        
215        if (data != NULL) {
216 <
216 >        
217          // make sure we can reinterpret the generic data as restraint data:
218 <
218 >        
219          RestraintData* restData= dynamic_cast<RestraintData*>(data);        
220 <
220 >        
221          if (restData != NULL) {
222 <
222 >          
223            // make sure we can reinterpet the restraint data as a
224            // pointer to a MolecularRestraint:
225 <
225 >          
226            MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
227 <
227 >          
228            if (mRest != NULL) {          
229 <
229 >            
230              // now we need to pack the stunt doubles for the reference
231              // structure:
232 <
232 >            
233              std::vector<Vector3d> ref;
234              int count = 0;
235              RealType mass, mTot;
# Line 140 | Line 247 | namespace OpenMD {
247                // doubles from the *globally* sorted array of
248                // positions:
249                
250 <              ref.push_back( all_pos_[sd->getGlobalIndex()] );
251 <              
145 <              mass = sd->getMass();
146 <              
250 >              ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] );
251 >              mass = sd->getMass();              
252                COM = COM + mass * ref[count];
253                mTot = mTot + mass;
254                count = count + 1;
255              }
151 –            
256              COM /= mTot;
153 –
257              mRest->setReferenceStructure(ref, COM);        
155 –
258            }
259          }
260        }
# Line 162 | Line 264 | namespace OpenMD {
264    
265    
266    void RestReader::parseDumpLine(const std::string& line) {        
267 +
268      StringTokenizer tokenizer(line);
269      int nTokens;
270      
# Line 169 | Line 272 | namespace OpenMD {
272      
273      if (nTokens < 2) {
274        sprintf(painCave.errMsg,
275 <              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
275 >              "RestReader Error: Not enough Tokens.\n%s\n", line.c_str());
276        painCave.isFatal = 1;
277        simError();
278      }
279  
280      int index = tokenizer.nextTokenAsInt();
281 <
281 >
282      StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
283  
284      if (integrableObject == NULL) {
285        return;
286 <    }  
287 <
286 >    }
287 >  
288      std::string type = tokenizer.nextToken();
289      int size = type.size();
290 +    
291      Vector3d pos;
188 –    Vector3d vel;
292      Quat4d q;
293 <    Vector3d ji;
191 <    Vector3d force;
192 <    Vector3d torque;
193 <          
293 >
294      for(int i = 0; i < size; ++i) {
295        switch(type[i]) {
196 –        
197 –        case 'p': {
198 –            pos[0] = tokenizer.nextTokenAsDouble();
199 –            pos[1] = tokenizer.nextTokenAsDouble();
200 –            pos[2] = tokenizer.nextTokenAsDouble();
201 –            break;
202 –        }
203 –        case 'v' : {
204 –            vel[0] = tokenizer.nextTokenAsDouble();
205 –            vel[1] = tokenizer.nextTokenAsDouble();
206 –            vel[2] = tokenizer.nextTokenAsDouble();
207 –            break;
208 –        }
296  
297 <        case 'q' : {
298 <           if (integrableObject->isDirectional()) {
299 <              
300 <             q[0] = tokenizer.nextTokenAsDouble();
301 <             q[1] = tokenizer.nextTokenAsDouble();
302 <             q[2] = tokenizer.nextTokenAsDouble();
303 <             q[3] = tokenizer.nextTokenAsDouble();
304 <              
305 <             RealType qlen = q.length();
306 <             if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
307 <                
308 <               sprintf(painCave.errMsg,
309 <                       "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
223 <               painCave.isFatal = 1;
224 <               simError();
225 <                
226 <             }  
227 <              
228 <             q.normalize();
229 <           }            
230 <           break;
231 <        }  
232 <        case 'j' : {
233 <          if (integrableObject->isDirectional()) {
234 <             ji[0] = tokenizer.nextTokenAsDouble();
235 <             ji[1] = tokenizer.nextTokenAsDouble();
236 <             ji[2] = tokenizer.nextTokenAsDouble();
237 <          }
238 <          break;
239 <        }  
240 <        case 'f': {
241 <          force[0] = tokenizer.nextTokenAsDouble();
242 <          force[1] = tokenizer.nextTokenAsDouble();
243 <          force[2] = tokenizer.nextTokenAsDouble();          
297 >      case 'p': {
298 >        pos[0] = tokenizer.nextTokenAsDouble();
299 >        pos[1] = tokenizer.nextTokenAsDouble();
300 >        pos[2] = tokenizer.nextTokenAsDouble();
301 >        break;
302 >      }
303 >      case 'v' : {
304 >        Vector3d vel;
305 >        vel[0] = tokenizer.nextTokenAsDouble();
306 >        vel[1] = tokenizer.nextTokenAsDouble();
307 >        vel[2] = tokenizer.nextTokenAsDouble();
308 >        break;
309 >      }
310  
311 <          break;
311 >      case 'q' : {
312 >        if (integrableObject->isDirectional()) {
313 >          
314 >          q[0] = tokenizer.nextTokenAsDouble();
315 >          q[1] = tokenizer.nextTokenAsDouble();
316 >          q[2] = tokenizer.nextTokenAsDouble();
317 >          q[3] = tokenizer.nextTokenAsDouble();
318 >          
319 >          RealType qlen = q.length();
320 >          if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
321 >            
322 >            sprintf(painCave.errMsg,
323 >                    "RestReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
324 >            painCave.isFatal = 1;
325 >            simError();            
326 >          }  
327 >              
328 >          q.normalize();
329 >        }          
330 >        break;
331 >      }  
332 >      case 'j' : {
333 >        Vector3d ji;
334 >        if (integrableObject->isDirectional()) {
335 >          ji[0] = tokenizer.nextTokenAsDouble();
336 >          ji[1] = tokenizer.nextTokenAsDouble();
337 >          ji[2] = tokenizer.nextTokenAsDouble();
338          }
339 <        case 't' : {
340 <           torque[0] = tokenizer.nextTokenAsDouble();
341 <           torque[1] = tokenizer.nextTokenAsDouble();
342 <           torque[2] = tokenizer.nextTokenAsDouble();          
339 >        break;
340 >      }  
341 >      case 'f': {        
342 >        Vector3d force;
343 >        force[0] = tokenizer.nextTokenAsDouble();
344 >        force[1] = tokenizer.nextTokenAsDouble();
345 >        force[2] = tokenizer.nextTokenAsDouble();          
346 >        break;
347 >      }
348 >      case 't' : {        
349 >        Vector3d torque;
350 >        torque[0] = tokenizer.nextTokenAsDouble();
351 >        torque[1] = tokenizer.nextTokenAsDouble();
352 >        torque[2] = tokenizer.nextTokenAsDouble();          
353 >        break;
354 >      }
355 >      default: {
356 >        sprintf(painCave.errMsg,
357 >                "RestReader Error: %s is an unrecognized type\n", type.c_str());
358 >        painCave.isFatal = 1;
359 >        simError();    
360 >        break;  
361 >      }
362 >      }
363 >      // keep the position in case we need it for a molecular restraint:
364  
365 <           break;
365 >      all_pos_[index] = pos;      
366 >        
367 >      // is this io restrained?
368 >      GenericData* data = integrableObject->getPropertyByName("Restraint");
369 >      ObjectRestraint* oRest;
370 >      
371 >      if (data != NULL) {
372 >        // make sure we can reinterpret the generic data as restraint data:
373 >        RestraintData* restData= dynamic_cast<RestraintData*>(data);        
374 >        if (restData != NULL) {
375 >          // make sure we can reinterpet the restraint data as a pointer to
376 >            // an ObjectRestraint:
377 >          oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
378 >          if (oRest != NULL) {          
379 >            if (integrableObject->isDirectional()) {
380 >              oRest->setReferenceStructure(pos, q.toRotationMatrix3());
381 >            } else {                          
382 >              oRest->setReferenceStructure(pos);
383 >            }
384 >          }
385          }
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 –          
386        }
387      }
388 <    
389 <    // keep the position in case we need it for a molecular restraint:
388 >  }
389 >  
390 >  void  RestReader::readStuntDoubles(std::istream& inputStream) {
391      
392 <    all_pos_[index] = pos;
393 <
394 <    // is this io restrained?    
395 <    GenericData* data = integrableObject->getPropertyByName("Restraint");
396 <    ObjectRestraint* oRest;
397 <
398 <    if (data != NULL) {
399 <      // make sure we can reinterpret the generic data as restraint data:
400 <      RestraintData* restData= dynamic_cast<RestraintData*>(data);        
401 <      if (restData != NULL) {
402 <        // make sure we can reinterpet the restraint data as a pointer to
403 <        // an ObjectRestraint:
404 <        oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
405 <        if (oRest != NULL) {          
406 <          if (integrableObject->isDirectional()) {
282 <            oRest->setReferenceStructure(pos, q.toRotationMatrix3());
283 <          } else {                          
284 <            oRest->setReferenceStructure(pos);
285 <          }
286 <        }
392 >    inputStream.getline(buffer, bufferSize);
393 >    std::string line(buffer);
394 >    
395 >    if (line.find("<StuntDoubles>") == std::string::npos) {
396 >      sprintf(painCave.errMsg,
397 >              "RestReader Error: Missing <StuntDoubles>\n");
398 >      painCave.isFatal = 1;
399 >      simError();
400 >    }
401 >    
402 >    while(inputStream.getline(buffer, bufferSize)) {
403 >      line = buffer;
404 >      
405 >      if(line.find("</StuntDoubles>") != std::string::npos) {
406 >        break;
407        }
408 +      
409 +      parseDumpLine(line);
410      }
411 <
411 >    
412    }
291 –  
413  
414 <
414 >  
415    void RestReader::readFrameProperties(std::istream& inputStream) {
416      inputStream.getline(buffer, bufferSize);
417      std::string line(buffer);
418  
419      if (line.find("<FrameData>") == std::string::npos) {
420        sprintf(painCave.errMsg,
421 <              "DumpReader Error: Missing <FrameData>\n");
421 >              "RestReader Error: Missing <FrameData>\n");
422        painCave.isFatal = 1;
423        simError();
424      }
# Line 314 | Line 435 | namespace OpenMD {
435        }
436        
437      }
438 <
438 >    
439    }
440  
441    

Diff Legend

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