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 996 by chrisfen, Wed Jun 28 14:35:14 2006 UTC vs.
Revision 1409 by gezelter, Fri Jan 22 21:31:12 2010 UTC

# Line 1 | Line 1
1 < /*
2 < * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 < *
1 > /*
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
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
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.
# Line 37 | Line 28
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 < */
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 > */
41  
42 < #define _LARGEFILE_SOURCE64
43 < #define _FILE_OFFSET_BITS 64
44 <
45 < #include <sys/types.h>
46 < #include <sys/stat.h>
47 <
48 < #include <iostream>
49 < #include <math.h>
50 <
51 < #include <stdio.h>
52 < #include <stdlib.h>
53 < #include <string.h>
54 <
55 < #include "primitives/Molecule.hpp"
56 < #include "utils/MemoryUtils.hpp"
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 "io/RestReader.hpp"
60 < #include "utils/simError.h"
59 > #include "restraints/ObjectRestraint.hpp"
60 > #include "restraints/MolecularRestraint.hpp"
61 >
62 > #ifdef IS_MPI
63 >
64 > #include <mpi.h>
65 > #endif
66  
67 < #ifdef IS_MPI
62 < #include <mpi.h>
63 < #define TAKE_THIS_TAG_CHAR 0
64 < #define TAKE_THIS_TAG_INT 1
65 < #define TAKE_THIS_TAG_DOUBLE 2
66 < #endif // is_mpi
67 > namespace OpenMD {
68  
69 < namespace oopse {
70 <  
71 <  RestReader::RestReader( SimInfo* info ) : info_(info){
69 >  void RestReader::scanFile(){
70 >    int lineNo = 0;
71 >    std::streampos prevPos;
72 >    std::streampos  currPos;
73      
74 <    idealName = "idealCrystal.in";
75 <    
76 < #ifdef IS_MPI
77 <    if (worldRank == 0) {
76 < #endif
74 > #ifdef IS_MPI
75 >    
76 >    if (worldRank == 0) {
77 > #endif // is_mpi
78  
79 <      inIdealFile = new std::ifstream(idealName.c_str());
80 <
81 <      if(inIdealFile->fail()){
82 <        sprintf(painCave.errMsg,
83 <                "RestReader: Cannot open file: %s\n",
84 <                idealName.c_str());
85 <        painCave.isFatal = 1;
86 <        simError();
79 >      inFile_->clear();
80 >      currPos = inFile_->tellg();
81 >      prevPos = currPos;
82 >      
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
97 > #ifdef IS_MPI
98      }
99 <    strcpy( checkPointMsg,
100 <            "File \"idealCrystal.in\" opened successfully for reading." );
101 <    MPIcheckPoint();
93 < #endif
99 >    MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD);
100 > #endif // is_mpi
101 >  }
102  
95 –    return;  
96 –  }
97 –  
98 –  RestReader :: ~RestReader( ){
99 – #ifdef IS_MPI
100 –    if (worldRank == 0) {
101 – #endif
103  
104 <      delete inIdealFile;
105 <      delete inAngFile;
105 <
106 < #ifdef IS_MPI
107 <    }
108 <    strcpy( checkPointMsg,
109 <            "File idealCrystal.in (and .zang0 if present) closed successfully." );
110 <    MPIcheckPoint();
111 < #endif
104 >  void RestReader::readSet(){
105 >    std::string line;
106      
107 <    return;
114 <  }
115 <  
116 <  
117 <  void RestReader :: readIdealCrystal(){
118 <        
119 < #ifdef IS_MPI
120 <    int which_node;
121 <    int i, j;
122 < #endif //is_mpi
107 > #ifndef IS_MPI
108      
109 <    const int BUFFERSIZE = 2000; // size of the read buffer
110 <    int nTotObjs; // the number of atoms
126 <    char read_buffer[BUFFERSIZE]; //the line buffer for reading
109 >    inFile_->clear();
110 >    inFile_->seekg(framePos_);
111      
112 <    char *parseErr;
112 >    std::istream& inputStream = *inFile_;
113 > #else
114      
115 <    std::vector<StuntDouble*> integrableObjects;
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 >      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 <    Molecule* mol;
150 <    StuntDouble* integrableObject;
134 <    SimInfo::MoleculeIterator mi;
135 <    Molecule::IntegrableObjectIterator ii;
149 >    std::istream& inputStream = sstream;  
150 > #endif
151      
152 < #ifndef IS_MPI
138 <    
139 <    inIdealFile->getline(read_buffer, sizeof(read_buffer));
152 >    inputStream.getline(buffer, bufferSize);
153  
154 <    if( inIdealFile->eof() ){
155 <      sprintf( painCave.errMsg,
156 <               "RestraintReader error: error reading 1st line of \"%s\"\n",
157 <               idealName.c_str() );
158 <      painCave.isFatal = 1;
159 <      simError();
160 <    }
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 <    nTotObjs = atoi( read_buffer );
162 >    //read frameData
163 >    readFrameProperties(inputStream);
164      
165 <    if( nTotObjs != info_->getNGlobalIntegrableObjects() ){
166 <      sprintf( painCave.errMsg,
153 <               "RestraintReader error. %s nIntegrable, %d, "
154 <               "does not match the meta-data file's nIntegrable, %d.\n",
155 <               idealName.c_str(),
156 <               nTotObjs,
157 <               info_->getNGlobalIntegrableObjects());
158 <      painCave.isFatal = 1;
159 <      simError();
160 <    }
165 >    //read StuntDoubles
166 >    readStuntDoubles(inputStream);
167      
168 <    // skip over the comment line
169 <    inIdealFile->getline(read_buffer, sizeof(read_buffer));
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  
180 <    if( inIdealFile->eof() ){
181 <      sprintf( painCave.errMsg,
182 <               "error in reading commment in %s\n",
183 <               idealName.c_str());
184 <      painCave.isFatal = 1;
185 <      simError();
171 <    }
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      
187 <    // parse the ideal crystal lines
188 <    /*
189 <     * Note: we assume that there is a one-to-one correspondence between
190 <     * integrable objects and lines in the idealCrystal.in file.  Thermodynamic
177 <     * integration is only supported for simple rigid bodies.
178 <     */
179 <    for (mol = info_->beginMolecule(mi); mol != NULL;
180 <         mol = info_->nextMolecule(mi)) {
181 <      
182 <      for (integrableObject = mol->beginIntegrableObject(ii);
183 <           integrableObject != NULL;
184 <           integrableObject = mol->nextIntegrableObject(ii)) {    
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 <        inIdealFile->getline(read_buffer, sizeof(read_buffer));
192 >    scanFile();
193  
194 <        if( inIdealFile->eof() ){
189 <          sprintf(painCave.errMsg,
190 <                  "RestReader Error: error in reading file %s\n"
191 <                  "natoms  = %d; index = %d\n"
192 <                  "error reading the line from the file.\n",
193 <                  idealName.c_str(), nTotObjs,
194 <                  integrableObject->getGlobalIndex() );
195 <          painCave.isFatal = 1;
196 <          simError();
197 <        }
198 <        
199 <        parseErr = parseIdealLine( read_buffer, integrableObject);
200 <        if( parseErr != NULL ){
201 <          strcpy( painCave.errMsg, parseErr );
202 <          painCave.isFatal = 1;
203 <          simError();
204 <        }
205 <      }
206 <    }
207 <    
208 <    // MPI Section of code..........
209 < #else //IS_MPI
210 <    
211 <    // first thing first, suspend fatalities.
212 <    painCave.isEventLoop = 1;
213 <    
214 <    int masterNode = 0;
215 <    
216 <    MPI_Status istatus;
217 <    int nCurObj;
218 <    int nitems;
219 <    int haveError;
194 >    readSet();
195  
221 –    nTotObjs = info_->getNGlobalIntegrableObjects();
222 –    haveError = 0;
196  
197 <    if (worldRank == masterNode) {
198 <      inIdealFile->getline(read_buffer, sizeof(read_buffer));
197 >    // all ObjectRestraints have been handled, now we have to worry about
198 >    // molecular restraints:
199  
200 <      if( inIdealFile->eof() ){
201 <        sprintf( painCave.errMsg,
202 <                 "Error reading 1st line of %s \n ",idealName.c_str());
203 <        painCave.isFatal = 1;
204 <        simError();
205 <      }
200 >    SimInfo::MoleculeIterator i;
201 >    Molecule::IntegrableObjectIterator j;
202 >    Molecule* mol;
203 >    StuntDouble* sd;
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 >    
209 >    for (mol = info_->beginMolecule(i); mol != NULL;
210 >         mol = info_->nextMolecule(i)) {          
211        
212 <      nitems = atoi( read_buffer );
212 >      // is this molecule restrained?    
213 >      GenericData* data = mol->getPropertyByName("Restraint");
214        
215 <      // Check to see that the number of integrable objects in the
237 <      // intial configuration file is the same as derived from the
238 <      // meta-data file.
239 <      if( nTotObjs != nitems){
240 <        sprintf( painCave.errMsg,
241 <                 "RestraintReader Error. %s nIntegrable, %d, "
242 <                 "does not match the meta-data file's nIntegrable, %d.\n",
243 <                 idealName.c_str(), nTotObjs,
244 <                 info_->getNGlobalIntegrableObjects());
245 <        painCave.isFatal = 1;
246 <        simError();
247 <      }
248 <      
249 <      // skip over the comment line
250 <      inIdealFile->getline(read_buffer, sizeof(read_buffer));
251 <
252 <      if( inIdealFile->eof() ){
253 <        sprintf( painCave.errMsg,
254 <                 "error in reading commment in %s\n", idealName.c_str());
255 <        painCave.isFatal = 1;
256 <        simError();
257 <      }
258 <
259 <      MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD);
260 <
261 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
262 <        which_node = info_->getMolToProc(i);
215 >      if (data != NULL) {
216          
217 <        if(which_node == masterNode){
218 <          //molecules belong to master node
217 >        // make sure we can reinterpret the generic data as restraint data:
218 >        
219 >        RestraintData* restData= dynamic_cast<RestraintData*>(data);        
220 >        
221 >        if (restData != NULL) {
222            
223 <          mol = info_->getMoleculeByGlobalIndex(i);
223 >          // make sure we can reinterpet the restraint data as a
224 >          // pointer to a MolecularRestraint:
225            
226 <          if(mol == NULL) {
270 <            sprintf(painCave.errMsg,
271 <                    "RestReader Error: Molecule not found on node %d!\n",
272 <                    worldRank);
273 <            painCave.isFatal = 1;
274 <            simError();
275 <          }
226 >          MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
227            
228 <          for (integrableObject = mol->beginIntegrableObject(ii);
278 <               integrableObject != NULL;
279 <               integrableObject = mol->nextIntegrableObject(ii)){
228 >          if (mRest != NULL) {          
229              
230 <            inIdealFile->getline(read_buffer, sizeof(read_buffer));
230 >            // now we need to pack the stunt doubles for the reference
231 >            // structure:
232              
233 <            if( inIdealFile->eof() ){
234 <              sprintf(painCave.errMsg,
235 <                      "RestReader Error: error in reading file %s\n"
236 <                      "natoms  = %d; index = %d\n"
237 <                      "error reading the line from the file.\n",
238 <                      idealName.c_str(), nTotObjs, i );
239 <              painCave.isFatal = 1;
240 <              simError();
233 >            std::vector<Vector3d> ref;
234 >            int count = 0;
235 >            RealType mass, mTot;
236 >            Vector3d COM(0.0);
237 >            
238 >            mTot = 0.0;
239 >            
240 >            // 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 >              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              }
256 <        
257 <            parseIdealLine(read_buffer, integrableObject);
294 <        
256 >            COM /= mTot;
257 >            mRest->setReferenceStructure(ref, COM);        
258            }
296 –
297 –        } else {
298 –          //molecule belongs to slave nodes
299 –          
300 –          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
301 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
302 –          
303 –          for(j = 0; j < nCurObj; j++){
304 –            inIdealFile->getline(read_buffer, sizeof(read_buffer));
305 –
306 –            if( inIdealFile->eof() ){
307 –              sprintf(painCave.errMsg,
308 –                      "RestReader Error: error in reading file %s\n"
309 –                      "natoms  = %d; index = %d\n"
310 –                      "error reading the line from the file.\n",
311 –                      idealName.c_str(), nTotObjs, i );
312 –              painCave.isFatal = 1;
313 –              simError();
314 –            }
315 –            
316 –            MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
317 –                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
318 –          }
259          }
260        }
261 <    } else {
262 <      //actions taken at slave nodes
323 <      MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD);
261 >    }
262 >  }
263  
264 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
265 <        int which_node = info_->getMolToProc(i);
264 >  
265 >  
266 >  void RestReader::parseDumpLine(const std::string& line) {        
267  
268 <        if(which_node == worldRank){
269 <          //molecule with global index i belongs to this processor
268 >    StringTokenizer tokenizer(line);
269 >    int nTokens;
270 >    
271 >    nTokens = tokenizer.countTokens();
272 >    
273 >    if (nTokens < 2) {
274 >      sprintf(painCave.errMsg,
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 >
282 >    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
283 >
284 >    if (integrableObject == NULL) {
285 >      return;
286 >    }
287 >  
288 >    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 >        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 >      case 'q' : {
312 >        if (integrableObject->isDirectional()) {
313            
314 <          mol = info_->getMoleculeByGlobalIndex(i);
314 >          q[0] = tokenizer.nextTokenAsDouble();
315 >          q[1] = tokenizer.nextTokenAsDouble();
316 >          q[2] = tokenizer.nextTokenAsDouble();
317 >          q[3] = tokenizer.nextTokenAsDouble();
318            
319 <          if(mol == NULL) {
320 <            sprintf(painCave.errMsg,
335 <                    "RestReader Error: molecule not found on node %d\n",
336 <                    worldRank);
337 <            painCave.isFatal = 1;
338 <            simError();
339 <          }
340 <          
341 <          nCurObj = mol->getNIntegrableObjects();
342 <          
343 <          MPI_Send(&nCurObj, 1, MPI_INT, masterNode,
344 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
345 <          
346 <          for (integrableObject = mol->beginIntegrableObject(ii);
347 <               integrableObject != NULL;
348 <               integrableObject = mol->nextIntegrableObject(ii)){
349 <            
350 <            MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode,
351 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
352 <            
353 <            parseErr = parseIdealLine(read_buffer, integrableObject);
319 >          RealType qlen = q.length();
320 >          if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
321              
322 <            if( parseErr != NULL ){
323 <              strcpy( painCave.errMsg, parseErr );
324 <              simError();
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 >        
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          }
386        }
387      }
363 – #endif
388    }
389    
390 <  char* RestReader::parseIdealLine(char* readLine, StuntDouble* sd){
390 >  void  RestReader::readStuntDoubles(std::istream& inputStream) {
391      
392 <    RealType pos[3];        // position place holders
393 <    RealType q[4];          // the quaternions
370 <    RealType RfromQ[3][3];  // the rotation matrix
371 <    RealType normalize;     // to normalize the reference unit vector
372 <    RealType uX, uY, uZ;    // reference unit vector place holders
373 <    RealType uselessToken;
374 <    StringTokenizer tokenizer(readLine);
375 <    int nTokens;
392 >    inputStream.getline(buffer, bufferSize);
393 >    std::string line(buffer);
394      
395 <    nTokens = tokenizer.countTokens();
396 <
397 <    if (nTokens < 14) {
398 <      sprintf(painCave.errMsg,
399 <              "RestReader Error: Not enough Tokens.\n");
382 <      painCave.isFatal = 1;
383 <      simError();
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 <    std::string name = tokenizer.nextToken();
403 <
388 <    if (name != sd->getType()) {
402 >    while(inputStream.getline(buffer, bufferSize)) {
403 >      line = buffer;
404        
405 <      sprintf(painCave.errMsg,
406 <              "RestReader Error: Atom type [%s] in %s does not "
407 <              "match Atom Type [%s] in .md file.\n",
393 <              name.c_str(), idealName.c_str(),
394 <              sd->getType().c_str());
395 <      painCave.isFatal = 1;
396 <      simError();        
397 <    }
398 <    
399 <    pos[0] = tokenizer.nextTokenAsDouble();
400 <    pos[1] = tokenizer.nextTokenAsDouble();
401 <    pos[2] = tokenizer.nextTokenAsDouble();
402 <
403 <    // store the positions in the stuntdouble as generic data doubles
404 <    DoubleGenericData* refPosX = new DoubleGenericData();
405 <    refPosX->setID("refPosX");
406 <    refPosX->setData(pos[0]);
407 <    sd->addProperty(refPosX);
408 <
409 <    DoubleGenericData* refPosY = new DoubleGenericData();
410 <    refPosY->setID("refPosY");
411 <    refPosY->setData(pos[1]);
412 <    sd->addProperty(refPosY);
413 <    
414 <    DoubleGenericData* refPosZ = new DoubleGenericData();
415 <    refPosZ->setID("refPosZ");
416 <    refPosZ->setData(pos[2]);
417 <    sd->addProperty(refPosZ);
418 <
419 <    // we don't need the velocities
420 <    uselessToken = tokenizer.nextTokenAsDouble();
421 <    uselessToken = tokenizer.nextTokenAsDouble();
422 <    uselessToken = tokenizer.nextTokenAsDouble();
423 <    
424 <    if (sd->isDirectional()) {
405 >      if(line.find("</StuntDoubles>") != std::string::npos) {
406 >        break;
407 >      }
408        
409 <      q[0] = tokenizer.nextTokenAsDouble();
427 <      q[1] = tokenizer.nextTokenAsDouble();
428 <      q[2] = tokenizer.nextTokenAsDouble();
429 <      q[3] = tokenizer.nextTokenAsDouble();
430 <      
431 <      // now build the rotation matrix and find the unit vectors
432 <      RfromQ[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
433 <      RfromQ[0][1] = 2*(q[1]*q[2] + q[0]*q[3]);
434 <      RfromQ[0][2] = 2*(q[1]*q[3] - q[0]*q[2]);
435 <      RfromQ[1][0] = 2*(q[1]*q[2] - q[0]*q[3]);
436 <      RfromQ[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
437 <      RfromQ[1][2] = 2*(q[2]*q[3] + q[0]*q[1]);
438 <      RfromQ[2][0] = 2*(q[1]*q[3] + q[0]*q[2]);
439 <      RfromQ[2][1] = 2*(q[2]*q[3] - q[0]*q[1]);
440 <      RfromQ[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
441 <      
442 <      normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1]
443 <                       + RfromQ[2][2]*RfromQ[2][2]);
444 <      uX = RfromQ[2][0]/normalize;
445 <      uY = RfromQ[2][1]/normalize;
446 <      uZ = RfromQ[2][2]/normalize;
447 <      
448 <      // store reference unit vectors as generic data in the stuntdouble
449 <      DoubleGenericData* refVectorX = new DoubleGenericData();
450 <      refVectorX->setID("refVectorX");
451 <      refVectorX->setData(uX);
452 <      sd->addProperty(refVectorX);
453 <      
454 <      DoubleGenericData* refVectorY = new DoubleGenericData();
455 <      refVectorY->setID("refVectorY");
456 <      refVectorY->setData(uY);
457 <      sd->addProperty(refVectorY);
458 <      
459 <      DoubleGenericData* refVectorZ = new DoubleGenericData();
460 <      refVectorZ->setID("refVectorZ");
461 <      refVectorZ->setData(uZ);
462 <      sd->addProperty(refVectorZ);
409 >      parseDumpLine(line);
410      }
411      
465 –    // we don't need the angular velocities, so let's exit the line
466 –    return NULL;
412    }
413 +
414    
415 <  void RestReader::readZangle(){
416 <    
417 <    int i;
472 <    int isPresent;
473 <    
474 <    Molecule* mol;
475 <    StuntDouble* integrableObject;
476 <    SimInfo::MoleculeIterator mi;
477 <    Molecule::IntegrableObjectIterator ii;
478 <    
479 < #ifdef IS_MPI
480 <    int which_node;
481 <    MPI_Status istatus;
482 < #endif //is_mpi
483 <    
484 <    const int BUFFERSIZE = 2000; // size of the read buffer
485 <    unsigned int nTotObjs; // the number of atoms
486 <    char read_buffer[BUFFERSIZE]; //the line buffer for reading
487 <    
488 <    std::vector<StuntDouble*> vecParticles;
489 <    std::vector<RealType> tempZangs;
490 <      
491 <    angFile = info_->getRestFileName();
492 <    
493 <    angFile += "0";
494 <    
495 <    // open the omega value file for reading
496 < #ifdef IS_MPI
497 <    if (worldRank == 0) {
498 < #endif
499 <      isPresent = 1;
500 <      
501 <      inAngFile = new std::ifstream(angFile.c_str());
502 <      
503 <      if(inAngFile->fail()){
504 <        sprintf(painCave.errMsg,
505 <                "Restraints Warning: %s file is not present\n"
506 <                "\tAll omega values will be initialized to zero. If the\n"
507 <                "\tsimulation is starting from the idealCrystal.in reference\n"
508 <                "\tconfiguration, this is the desired action. If this is not\n"
509 <                "\tthe case, the energy calculations will be incorrect.\n",
510 <                angFile.c_str());
511 <        painCave.severity = OOPSE_WARNING;
512 <        painCave.isFatal = 0;
513 <        simError();  
514 <        isPresent = 0;
515 <      }
516 <      
517 < #ifdef IS_MPI
518 <      // let the other nodes know the status of the file search
519 <      MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
520 < #endif // is_mpi
521 <      
522 <      if (!isPresent) {
523 <        zeroZangle();
524 <        return;
525 <      }
526 <      
527 < #ifdef IS_MPI
528 <      if (!isPresent) {
529 <        // master node zeroes out its zAngles if .zang0 isn't present
530 <        zeroZangle();
531 <        return;
532 <      }
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 +              "RestReader Error: Missing <FrameData>\n");
422 +      painCave.isFatal = 1;
423 +      simError();
424      }
425 <    
426 <    // listen to node 0 to see if we should exit this function
427 <    if (worldRank != 0) {
428 <      MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
429 <      if (!isPresent) {
430 <        zeroZangle();
431 <        return;
425 >
426 >    // 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 >      
433 >      if(line.find("</FrameData>") != std::string::npos) {
434 >        break;
435        }
543 –    }
544 –    
545 –    strcpy( checkPointMsg, "zAngle file opened successfully for reading." );
546 –    MPIcheckPoint();
547 – #endif
548 –    
549 – #ifndef IS_MPI
550 –    
551 –    // read the first line and die if there is a failure
552 –    inAngFile->getline(read_buffer, sizeof(read_buffer));
553 –
554 –    if( inAngFile->eof() ){
555 –      sprintf( painCave.errMsg,
556 –               "RestraintReader error: error reading 1st line of \"%s\"\n",
557 –               angFile.c_str() );
558 –      painCave.isFatal = 1;
559 –      simError();
560 –    }
561 –
562 –    // read the file and load the values into a vector
563 –    inAngFile->getline(read_buffer, sizeof(read_buffer));
564 –    
565 –    while ( !inAngFile->eof() ) {
566 –      // check for and ignore blank lines
567 –      if ( read_buffer != NULL )
568 –        tempZangs.push_back( atof(read_buffer) );
569 –
570 –      inAngFile->getline(read_buffer, sizeof(read_buffer));
571 –    }
572 –    
573 –    nTotObjs = info_->getNGlobalIntegrableObjects();
574 –    
575 –    if( nTotObjs != tempZangs.size() ){
576 –      sprintf( painCave.errMsg,
577 –               "RestraintReader zAngle reading error. %s nIntegrable, %d, "
578 –               "does not match the meta-data file's nIntegrable, %i.\n",
579 –               angFile.c_str(),
580 –               tempZangs.size(),
581 –               nTotObjs );
582 –      painCave.isFatal = 1;
583 –      simError();
584 –    }
585 –    
586 –    // load the zAngles into the integrable objects
587 –    i = 0;
588 –    for (mol = info_->beginMolecule(mi); mol != NULL;
589 –         mol = info_->nextMolecule(mi)) {
436        
591 –      for (integrableObject = mol->beginIntegrableObject(ii);
592 –           integrableObject != NULL;
593 –           integrableObject = mol->nextIntegrableObject(ii)) {    
594 –        
595 –        integrableObject->setZangle(tempZangs[i]);
596 –        i++;
597 –      }
437      }
438      
600 –    // MPI Section of code..........
601 – #else //IS_MPI
602 –    
603 –    // first thing first, suspend fatalities.
604 –    painCave.isEventLoop = 1;
605 –
606 –    int masterNode = 0;
607 –
608 –    int haveError;
609 –    int intObjIndex;    
610 –    int intObjIndexTransfer;    
611 –
612 –    int j;
613 –    int nCurObj;
614 –    RealType angleTransfer;
615 –    
616 –    nTotObjs = info_->getNGlobalIntegrableObjects();
617 –    haveError = 0;
618 –
619 –    if (worldRank == masterNode) {
620 –      
621 –      inAngFile->getline(read_buffer, sizeof(read_buffer));
622 –
623 –      if( inAngFile->eof() ){
624 –        sprintf( painCave.errMsg,
625 –                 "Error reading 1st line of %s \n ",angFile.c_str());
626 –        haveError = 1;
627 –        simError();
628 –      }
629 –      
630 –      // let the master node read the file and load the temporary angle vector
631 –      inAngFile->getline(read_buffer, sizeof(read_buffer));
632 –
633 –      while ( !inAngFile->eof() ) {
634 –        // check for and ignore blank lines
635 –        if ( read_buffer != NULL )
636 –          tempZangs.push_back( atof(read_buffer) );
637 –
638 –        inAngFile->getline(read_buffer, sizeof(read_buffer));
639 –      }
640 –
641 –      // Check to see that the number of integrable objects in the
642 –      // intial configuration file is the same as derived from the
643 –      // meta-data file.
644 –      if( nTotObjs != tempZangs.size() ){
645 –        sprintf( painCave.errMsg,
646 –                 "RestraintReader zAngle reading Error. %s nIntegrable, %d, "
647 –                 "does not match the meta-data file's nIntegrable, %d.\n",
648 –                 angFile.c_str(),
649 –                 tempZangs.size(),
650 –                 nTotObjs);
651 –        haveError= 1;
652 –        simError();
653 –      }
654 –      
655 –      // At this point, node 0 has a tempZangs vector completed, and
656 –      // everyone else has nada
657 –      
658 –      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
659 –        // Get the Node number which has this atom
660 –        which_node = info_->getMolToProc(i);
661 –        
662 –        if (which_node == masterNode) {
663 –          mol = info_->getMoleculeByGlobalIndex(i);
664 –
665 –          if(mol == NULL) {
666 –            strcpy(painCave.errMsg, "Molecule not found on node 0!");
667 –            haveError = 1;
668 –            simError();
669 –          }
670 –
671 –          for (integrableObject = mol->beginIntegrableObject(ii);
672 –               integrableObject != NULL;
673 –               integrableObject = mol->nextIntegrableObject(ii)){
674 –            intObjIndex = integrableObject->getGlobalIndex();
675 –            integrableObject->setZangle(tempZangs[intObjIndex]);
676 –          }    
677 –          
678 –        } else {
679 –          // I am MASTER OF THE UNIVERSE, but I don't own this molecule
680 –          // listen for the number of integrableObjects in the molecule
681 –          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
682 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
683 –          
684 –          for(j=0; j < nCurObj; j++){          
685 –            // listen for which integrableObject we need to send the value for
686 –            MPI_Recv(&intObjIndexTransfer, 1, MPI_INT, which_node,
687 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
688 –            angleTransfer = tempZangs[intObjIndexTransfer];
689 –            // send the value to the node so it can initialize the object
690 –            MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node,
691 –                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);              
692 –          }
693 –        }
694 –      }
695 –    } else {
696 –      // I am SLAVE TO THE MASTER
697 –      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
698 –        which_node = info_->getMolToProc(i);
699 –
700 –        if (which_node == worldRank) {
701 –          
702 –          // BUT I OWN THIS MOLECULE!!!
703 –          
704 –          mol = info_->getMoleculeByGlobalIndex(i);
705 –
706 –          if(mol == NULL) {
707 –            sprintf(painCave.errMsg,
708 –                    "RestReader Error: molecule not found on node %d\n",
709 –                    worldRank);
710 –            painCave.isFatal = 1;
711 –            simError();
712 –          }
713 –
714 –          nCurObj = mol->getNIntegrableObjects();
715 –          // send the number of integrableObjects in the molecule
716 –          MPI_Send(&nCurObj, 1, MPI_INT, 0,
717 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
718 –          
719 –          for (integrableObject = mol->beginIntegrableObject(ii);
720 –               integrableObject != NULL;
721 –               integrableObject = mol->nextIntegrableObject(ii)){
722 –            intObjIndexTransfer = integrableObject->getGlobalIndex();
723 –            // send the global index of the integrableObject
724 –            MPI_Send(&intObjIndexTransfer, 1, MPI_INT, 0,
725 –                     TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
726 –            // listen for the value we want to set locally
727 –            MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0,
728 –                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
729 –
730 –            integrableObject->setZangle(angleTransfer);
731 –          }    
732 –        }
733 –      }
734 –    }
735 – #endif
439    }
737 –  
738 –  void RestReader :: zeroZangle(){
739 –    
740 –    Molecule* mol;
741 –    StuntDouble* integrableObject;
742 –    SimInfo::MoleculeIterator mi;
743 –    Molecule::IntegrableObjectIterator ii;
744 –    
745 – #ifndef IS_MPI
746 –    // set all zAngles to 0.0
747 –    for (mol = info_->beginMolecule(mi); mol != NULL;
748 –         mol = info_->nextMolecule(mi))
749 –      
750 –      for (integrableObject = mol->beginIntegrableObject(ii);
751 –           integrableObject != NULL;
752 –           integrableObject = mol->nextIntegrableObject(ii))    
753 –        integrableObject->setZangle( 0.0 );
754 –    
755 –    
756 –    // MPI Section of code..........
757 – #else //IS_MPI
758 –    
759 –    // first thing first, suspend fatalities.
760 –    painCave.isEventLoop = 1;
761 –    
762 –    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
763 –    int haveError;
764 –    int which_node;
765 –    
766 –    haveError = 0;
440  
441 <    for (int i=0 ; i < info_->getNGlobalMolecules(); i++) {
442 <
770 <      // Get the Node number which has this atom
771 <      which_node = info_->getMolToProc(i);
772 <        
773 <      // each processor zeroes its own integrable objects
774 <      if (which_node == worldRank) {
775 <        mol = info_->getMoleculeByGlobalIndex(i);
776 <        
777 <        if(mol == NULL) {
778 <          sprintf( painCave.errMsg,
779 <                  "Molecule not found on node %i!",
780 <                  which_node );
781 <          haveError = 1;
782 <          simError();
783 <        }
784 <        
785 <        for (integrableObject = mol->beginIntegrableObject(ii);
786 <             integrableObject != NULL;
787 <             integrableObject = mol->nextIntegrableObject(ii)){
788 <          
789 <          integrableObject->setZangle( 0.0 );
790 <          
791 <        }
792 <      }
793 <    }
794 <
795 < #endif
796 <  }
797 <  
798 < } // end namespace oopse
441 >  
442 > }//end namespace OpenMD

Diff Legend

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