src/io/RestReader.cpp

/* 
 * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved. 
 * 
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Vardeman & Gezelter, in progress (2009).                        
 */ 

 
#include <sys/types.h> 
#include <sys/stat.h> 
 
#include <math.h> 
#include <string>
#include <sstream>
#include <iostream> 
#include <stdio.h> 
#include <stdlib.h> 
#include <string.h> 
 
#include "io/RestReader.hpp" 
#include "primitives/Molecule.hpp" 
#include "utils/simError.h" 
#include "utils/MemoryUtils.hpp" 
#include "utils/StringTokenizer.hpp"
#include "restraints/ObjectRestraint.hpp"
#include "restraints/MolecularRestraint.hpp" 
 
#ifdef IS_MPI 
 
#include <mpi.h> 
#endif

namespace OpenMD { 

  void RestReader::scanFile(){
    int lineNo = 0; 
    std::streampos prevPos;
    std::streampos  currPos;
    
#ifdef IS_MPI 
    
    if (worldRank == 0) { 
#endif // is_mpi 

      inFile_->clear();
      currPos = inFile_->tellg();
      prevPos = currPos;
      
      bool foundOpenSnapshotTag = false;
      
      while(!foundOpenSnapshotTag && inFile_->getline(buffer, bufferSize)) {
        ++lineNo;
        
        std::string line = buffer;
        currPos = inFile_->tellg(); 
        if (line.find("<Snapshot>")!= std::string::npos) {
          foundOpenSnapshotTag = true;
          framePos_ = prevPos;
        }
        prevPos = currPos;
      }
      
#ifdef IS_MPI 
    }
    MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD); 
#endif // is_mpi 
  } 


  void RestReader::readSet(){
    std::string line;
    
#ifndef IS_MPI
    
    inFile_->clear();
    inFile_->seekg(framePos_);
    
    std::istream& inputStream = *inFile_;
#else
    
    int masterNode = 0;
    std::stringstream sstream;
    if (worldRank == masterNode) {
      std::string sendBuffer;
      
      inFile_->clear();  
      inFile_->seekg(framePos_); 
      
      while (inFile_->getline(buffer, bufferSize)) {
        
        line = buffer;
        sendBuffer += line;
        sendBuffer += '\n';
        if (line.find("</Snapshot>") != std::string::npos) {
          break;
        }        
      }
      
      int sendBufferSize = sendBuffer.size();
      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);     
      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);     
      
      sstream.str(sendBuffer);
    } else {
      int sendBufferSize;
      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);     
      char * recvBuffer = new char[sendBufferSize+1];
      assert(recvBuffer);
      recvBuffer[sendBufferSize] = '\0';
      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);     
      sstream.str(recvBuffer);
      delete [] recvBuffer;
    }      
    
    std::istream& inputStream = sstream;  
#endif
    
    inputStream.getline(buffer, bufferSize);

    line = buffer;
    if (line.find("<Snapshot>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "RestReader Error: can not find <Snapshot>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    } 
    
    //read frameData
    readFrameProperties(inputStream);
    
    //read StuntDoubles
    readStuntDoubles(inputStream);
    
    inputStream.getline(buffer, bufferSize);
    line = buffer;
    if (line.find("</Snapshot>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "RestReader Error: can not find </Snapshot>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }            
  }
    
  void RestReader::readReferenceStructure() {

    // We need temporary storage to keep track of all StuntDouble positions
    // in case some of the restraints are molecular (i.e. if they use
    // multiple SD positions to determine restrained orientations or positions:

    all_pos_.clear();
    all_pos_.resize(info_->getNGlobalIntegrableObjects()) ;
    
    // Restraint files are just standard dump files, but with the reference
    // structure stored in the first frame (frame 0).
    // RestReader overloads readSet and explicitly handles all of the
    // ObjectRestraints in that method:

    scanFile();

    readSet();


    // all ObjectRestraints have been handled, now we have to worry about
    // molecular restraints:

    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator j;
    Molecule* mol;
    StuntDouble* sd;
    
    // no need to worry about parallel molecules, as molecules are not
    // split across processor boundaries.  Just loop over all molecules
    // we know about:
    
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {          
      
      // is this molecule restrained?    
      GenericData* data = mol->getPropertyByName("Restraint");
      
      if (data != NULL) {
        
        // make sure we can reinterpret the generic data as restraint data:
        
        RestraintData* restData= dynamic_cast<RestraintData*>(data);        
        
        if (restData != NULL) {
          
          // make sure we can reinterpet the restraint data as a
          // pointer to a MolecularRestraint:
          
          MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
          
          if (mRest != NULL) {          
            
            // now we need to pack the stunt doubles for the reference
            // structure:
            
            std::vector<Vector3d> ref;
            int count = 0;
            RealType mass, mTot;
            Vector3d COM(0.0);
            
            mTot = 0.0;
            
            // loop over the stunt doubles in this molecule in the order we
            // will be looping them in the restraint code:
            
            for (sd = mol->beginIntegrableObject(j); sd != NULL; 
                 sd = mol->nextIntegrableObject(j)) {
              
              // push back the reference positions of the stunt
              // doubles from the *globally* sorted array of
              // positions:
              
              ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] );
              mass = sd->getMass();              
              COM = COM + mass * ref[count];
              mTot = mTot + mass;
              count = count + 1;
            }
            COM /= mTot;
            mRest->setReferenceStructure(ref, COM);         
          }
        }
      }
    }
  }

   
  void RestReader::parseDumpLine(const std::string& line) {        

    StringTokenizer tokenizer(line); 
    int nTokens; 
     
    nTokens = tokenizer.countTokens(); 
     
    if (nTokens < 2) { 
      sprintf(painCave.errMsg, 
              "RestReader Error: Not enough Tokens.\n%s\n", line.c_str()); 
      painCave.isFatal = 1; 
      simError(); 
    } 

    int index = tokenizer.nextTokenAsInt();

    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);

    if (integrableObject == NULL) {
      return;
    }
  
    std::string type = tokenizer.nextToken(); 
    int size = type.size();
    
    Vector3d pos;
    Quat4d q;

    for(int i = 0; i < size; ++i) {
      switch(type[i]) {

      case 'p': {
        pos[0] = tokenizer.nextTokenAsDouble(); 
        pos[1] = tokenizer.nextTokenAsDouble(); 
        pos[2] = tokenizer.nextTokenAsDouble(); 
        break;
      }
      case 'v' : {
        Vector3d vel;
        vel[0] = tokenizer.nextTokenAsDouble(); 
        vel[1] = tokenizer.nextTokenAsDouble(); 
        vel[2] = tokenizer.nextTokenAsDouble(); 
        break;
      }

      case 'q' : {
        if (integrableObject->isDirectional()) { 
          
          q[0] = tokenizer.nextTokenAsDouble(); 
          q[1] = tokenizer.nextTokenAsDouble(); 
          q[2] = tokenizer.nextTokenAsDouble(); 
          q[3] = tokenizer.nextTokenAsDouble(); 
          
          RealType qlen = q.length(); 
          if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0 
            
            sprintf(painCave.errMsg, 
                    "RestReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n"); 
            painCave.isFatal = 1; 
            simError();             
          }  
              
          q.normalize(); 
        }          
        break;
      }  
      case 'j' : {
        Vector3d ji;
        if (integrableObject->isDirectional()) {
          ji[0] = tokenizer.nextTokenAsDouble(); 
          ji[1] = tokenizer.nextTokenAsDouble(); 
          ji[2] = tokenizer.nextTokenAsDouble(); 
        }
        break;
      }  
      case 'f': {        
        Vector3d force;
        force[0] = tokenizer.nextTokenAsDouble(); 
        force[1] = tokenizer.nextTokenAsDouble(); 
        force[2] = tokenizer.nextTokenAsDouble();           
        break;
      }
      case 't' : {        
        Vector3d torque;
        torque[0] = tokenizer.nextTokenAsDouble(); 
        torque[1] = tokenizer.nextTokenAsDouble(); 
        torque[2] = tokenizer.nextTokenAsDouble();           
        break;
      }
      default: {
        sprintf(painCave.errMsg, 
                "RestReader Error: %s is an unrecognized type\n", type.c_str()); 
        painCave.isFatal = 1; 
        simError();     
        break;   
      }
      }
      // keep the position in case we need it for a molecular restraint:

      all_pos_[index] = pos;      
        
      // is this io restrained?
      GenericData* data = integrableObject->getPropertyByName("Restraint");
      ObjectRestraint* oRest;
      
      if (data != NULL) {
        // make sure we can reinterpret the generic data as restraint data:
        RestraintData* restData= dynamic_cast<RestraintData*>(data);        
        if (restData != NULL) {
          // make sure we can reinterpet the restraint data as a pointer to
            // an ObjectRestraint:
          oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
          if (oRest != NULL) {          
            if (integrableObject->isDirectional()) {
              oRest->setReferenceStructure(pos, q.toRotationMatrix3());
            } else {                           
              oRest->setReferenceStructure(pos);
            }
          }
        }
      }
    }
  }
  
  void  RestReader::readStuntDoubles(std::istream& inputStream) {
    
    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);
    
    if (line.find("<StuntDoubles>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "RestReader Error: Missing <StuntDoubles>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }
    
    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</StuntDoubles>") != std::string::npos) {
        break;
      }
      
      parseDumpLine(line);
    }
    
  }

  
  void RestReader::readFrameProperties(std::istream& inputStream) {
    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);

    if (line.find("<FrameData>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "RestReader Error: Missing <FrameData>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }

    // restraints don't care about frame data (unless we need to wrap
    // coordinates, but we'll worry about that later), so 
    // we'll just scan ahead until the end of the frame data:

    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</FrameData>") != std::string::npos) {
        break;
      }
      
    }
    
  }

   
}//end namespace OpenMD 
Revision:	1465
Committed:	Fri Jul 9 23:08:25 2010 UTC (15 years, 3 months ago) by chuckv
File size:	13522 byte(s)
Log Message:	Creating busticated version of OpenMD
#	User	Rev	Content
1	cli2	1360	/*
2	cli2	1407	* Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
3	cli2	1360	*
4	gezelter	1390	* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Redistributions of source code must retain the above copyright
10			* notice, this list of conditions and the following disclaimer.
11			*
12			* 2. Redistributions in binary form must reproduce the above copyright
13			* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31			*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
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).
40	cli2	1360	*/
41	cli2	1407
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	gezelter	1409	#include "utils/StringTokenizer.hpp"
59	cli2	1360	#include "restraints/ObjectRestraint.hpp"
60	cli2	1407	#include "restraints/MolecularRestraint.hpp"
61
62			#ifdef IS_MPI
63
64			#include <mpi.h>
65			#endif
66	chrisfen	417
67	gezelter	1390	namespace OpenMD {
68	cli2	1407
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	cli2	1360
83	cli2	1407	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	chrisfen	417
103
104	cli2	1407	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	cli2	1360	} else {
139	cli2	1407	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	cli2	1360	}
161	cli2	1407
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	chrisfen	417
180	cli2	1360	// 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	chrisfen	417
184	cli2	1360	all_pos_.clear();
185	cli2	1407	all_pos_.resize(info_->getNGlobalIntegrableObjects()) ;
186
187	cli2	1360	// 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	chrisfen	996
192	cli2	1407	scanFile();
193
194			readSet();
195
196
197	cli2	1360	// all ObjectRestraints have been handled, now we have to worry about
198			// molecular restraints:
199	chrisfen	431
200	cli2	1360	SimInfo::MoleculeIterator i;
201			Molecule::IntegrableObjectIterator j;
202			Molecule* mol;
203			StuntDouble* sd;
204	cli2	1407
205	cli2	1360	// 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	cli2	1407
209	cli2	1360	for (mol = info_->beginMolecule(i); mol != NULL;
210			mol = info_->nextMolecule(i)) {
211	cli2	1407
212	cli2	1360	// is this molecule restrained?
213			GenericData* data = mol->getPropertyByName("Restraint");
214	chrisfen	417
215	cli2	1360	if (data != NULL) {
216	cli2	1407
217	cli2	1360	// make sure we can reinterpret the generic data as restraint data:
218	cli2	1407
219	cli2	1360	RestraintData* restData= dynamic_cast<RestraintData*>(data);
220	cli2	1407
221	cli2	1360	if (restData != NULL) {
222	cli2	1407
223	cli2	1360	// make sure we can reinterpet the restraint data as a
224			// pointer to a MolecularRestraint:
225	cli2	1407
226	cli2	1360	MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
227	cli2	1407
228	cli2	1360	if (mRest != NULL) {
229	cli2	1407
230	cli2	1360	// now we need to pack the stunt doubles for the reference
231			// structure:
232	cli2	1407
233	cli2	1360	std::vector<Vector3d> ref;
234			int count = 0;
235			RealType mass, mTot;
236			Vector3d COM(0.0);
237	chrisfen	417
238	cli2	1360	mTot = 0.0;
239	chrisfen	417
240	cli2	1360	// loop over the stunt doubles in this molecule in the order we
241			// will be looping them in the restraint code:
242
243			for (sd = mol->beginIntegrableObject(j); sd != NULL;
244			sd = mol->nextIntegrableObject(j)) {
245
246			// push back the reference positions of the stunt
247			// doubles from the globally sorted array of
248			// positions:
249
250	cli2	1407	ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] );
251			mass = sd->getMass();
252	cli2	1360	COM = COM + mass * ref[count];
253			mTot = mTot + mass;
254			count = count + 1;
255	chrisfen	417	}
256	cli2	1360	COM /= mTot;
257			mRest->setReferenceStructure(ref, COM);
258	chrisfen	417	}
259			}
260			}
261			}
262			}
263	chrisfen	431
264	cli2	1360
265
266			void RestReader::parseDumpLine(const std::string& line) {
267	cli2	1407
268	cli2	1360	StringTokenizer tokenizer(line);
269			int nTokens;
270
271			nTokens = tokenizer.countTokens();
272
273			if (nTokens < 2) {
274			sprintf(painCave.errMsg,
275	cli2	1407	"RestReader Error: Not enough Tokens.\n%s\n", line.c_str());
276	cli2	1360	painCave.isFatal = 1;
277			simError();
278			}
279	chrisfen	431
280	cli2	1360	int index = tokenizer.nextTokenAsInt();
281	chrisfen	431
282	gezelter	1408	StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
283
284			if (integrableObject == NULL) {
285			return;
286			}
287	cli2	1407
288	gezelter	1408	std::string type = tokenizer.nextToken();
289			int size = type.size();
290
291			Vector3d pos;
292			Quat4d q;
293
294			for(int i = 0; i < size; ++i) {
295			switch(type[i]) {
296
297			case 'p': {
298	cli2	1407	pos[0] = tokenizer.nextTokenAsDouble();
299			pos[1] = tokenizer.nextTokenAsDouble();
300			pos[2] = tokenizer.nextTokenAsDouble();
301	gezelter	1408	break;
302			}
303			case 'v' : {
304	cli2	1407	Vector3d vel;
305			vel[0] = tokenizer.nextTokenAsDouble();
306			vel[1] = tokenizer.nextTokenAsDouble();
307			vel[2] = tokenizer.nextTokenAsDouble();
308	gezelter	1408	break;
309			}
310
311			case 'q' : {
312	cli2	1407	if (integrableObject->isDirectional()) {
313	cli2	1360
314	cli2	1407	q[0] = tokenizer.nextTokenAsDouble();
315			q[1] = tokenizer.nextTokenAsDouble();
316			q[2] = tokenizer.nextTokenAsDouble();
317			q[3] = tokenizer.nextTokenAsDouble();
318	gezelter	1408
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			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			all_pos_[index] = pos;
366	cli2	1360
367	gezelter	1408	// 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	cli2	1407	// an ObjectRestraint:
377	gezelter	1408	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	cli2	1407	}
384	cli2	1360	}
385			}
386	chrisfen	417	}
387			}
388	cli2	1407	}
389
390			void RestReader::readStuntDoubles(std::istream& inputStream) {
391	chrisfen	417
392	cli2	1407	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	chrisfen	417	}
408	cli2	1407
409			parseDumpLine(line);
410	chrisfen	417	}
411	cli2	1407
412	cli2	1360	}
413	chrisfen	996
414	cli2	1407
415	cli2	1360	void RestReader::readFrameProperties(std::istream& inputStream) {
416			inputStream.getline(buffer, bufferSize);
417			std::string line(buffer);
418	chrisfen	423
419	cli2	1360	if (line.find("<FrameData>") == std::string::npos) {
420			sprintf(painCave.errMsg,
421	cli2	1407	"RestReader Error: Missing <FrameData>\n");
422	cli2	1360	painCave.isFatal = 1;
423			simError();
424			}
425	chrisfen	996
426	cli2	1360	// restraints don't care about frame data (unless we need to wrap
427			// coordinates, but we'll worry about that later), so
428			// we'll just scan ahead until the end of the frame data:
429
430			while(inputStream.getline(buffer, bufferSize)) {
431			line = buffer;
432	chrisfen	417
433	cli2	1360	if(line.find("</FrameData>") != std::string::npos) {
434			break;
435	chrisfen	417	}
436
437	cli2	1360	}
438	cli2	1407
439	chrisfen	417	}
440	chrisfen	1030
441	cli2	1360
442	gezelter	1390	}//end namespace OpenMD