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 990 by chrisfen, Mon Jun 19 01:36:06 2006 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.
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";
74 > #ifdef IS_MPI
75      
76 <    isScanned = false;
77 <    
78 < #ifdef IS_MPI
79 <    if (worldRank == 0) {
80 < #endif
76 >    if (worldRank == 0) {
77 > #endif // is_mpi
78 >
79 >      inFile_->clear();
80 >      currPos = inFile_->tellg();
81 >      prevPos = currPos;
82        
83 <      inIdealFile = fopen(idealName, "r");
84 <      if(inIdealFile == NULL){
85 <        sprintf(painCave.errMsg,
86 <                "RestReader: Cannot open file: %s\n", idealName);
87 <        painCave.isFatal = 1;
88 <        simError();
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 <      inIdealFileName = idealName;
89 < #ifdef IS_MPI
97 > #ifdef IS_MPI
98      }
99 <    strcpy( checkPointMsg,
100 <            "File \"idealCrystal.in\" opened successfully for reading." );
101 <    MPIcheckPoint();
94 < #endif
99 >    MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD);
100 > #endif // is_mpi
101 >  }
102  
103 <    return;  
104 <  }
105 <  
106 <  RestReader :: ~RestReader( ){
107 < #ifdef IS_MPI
108 <    if (worldRank == 0) {
109 < #endif
110 <      int error;
111 <      error = fclose( inIdealFile );
103 >
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 <      if( error ){
121 <        sprintf( painCave.errMsg,
122 <                 "Error closing %s\n", inIdealFileName.c_str());
123 <        simError();
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 <      MemoryUtils::deletePointers(framePos);
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 < #ifdef IS_MPI
138 <    }
139 <    strcpy( checkPointMsg, "Restraint file closed successfully." );
140 <    MPIcheckPoint();
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 >    std::istream& inputStream = sstream;  
150   #endif
151      
152 <    return;
121 <  }
122 <  
123 <  
124 <  void RestReader :: readIdealCrystal(){
125 <        
126 <    int i;
127 <    unsigned int j;
152 >    inputStream.getline(buffer, bufferSize);
153  
154 < #ifdef IS_MPI
155 <    int done, which_node, which_atom; // loop counter
156 < #endif //is_mpi
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 <    const int BUFFERSIZE = 2000; // size of the read buffer
163 <    int nTotObjs; // the number of atoms
135 <    char read_buffer[BUFFERSIZE]; //the line buffer for reading
162 >    //read frameData
163 >    readFrameProperties(inputStream);
164      
165 <    char *eof_test; // ptr to see when we reach the end of the file
166 <    char *parseErr;
165 >    //read StuntDoubles
166 >    readStuntDoubles(inputStream);
167      
168 <    std::vector<StuntDouble*> integrableObjects;
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 +    // 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 +    // 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 +    scanFile();
193 +
194 +
195 +    readSet();
196 +
197 +
198 +    // all ObjectRestraints have been handled, now we have to worry about
199 +    // molecular restraints:
200 +
201 +    SimInfo::MoleculeIterator i;
202 +    Molecule::IntegrableObjectIterator j;
203      Molecule* mol;
204 <    StuntDouble* integrableObject;
144 <    SimInfo::MoleculeIterator mi;
145 <    Molecule::IntegrableObjectIterator ii;
204 >    StuntDouble* sd;
205      
206 < #ifndef IS_MPI
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      
210 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
211 <    if( eof_test == NULL ){
151 <      sprintf( painCave.errMsg,
152 <               "RestraintReader error: error reading 1st line of \"%s\"\n",
153 <               inIdealFileName.c_str() );
154 <      painCave.isFatal = 1;
155 <      simError();
156 <    }
157 <    
158 <    nTotObjs = atoi( read_buffer );
159 <    
160 <    if( nTotObjs != info_->getNGlobalIntegrableObjects() ){
161 <      sprintf( painCave.errMsg,
162 <               "RestraintReader error. %s nIntegrable, %d, "
163 <               "does not match the meta-data file's nIntegrable, %d.\n",
164 <               inIdealFileName.c_str(), nTotObjs,
165 <               info_->getNGlobalIntegrableObjects());
166 <      painCave.isFatal = 1;
167 <      simError();
168 <    }
169 <    
170 <    // skip over the comment line
171 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
172 <    if(eof_test == NULL){
173 <      sprintf( painCave.errMsg,
174 <               "error in reading commment in %s\n", inIdealFileName.c_str());
175 <      painCave.isFatal = 1;
176 <      simError();
177 <    }
178 <    
179 <    // parse the ideal crystal lines
180 <    /*
181 <     * Note: we assume that there is a one-to-one correspondence between
182 <     * integrable objects and lines in the idealCrystal.in file.  Thermodynamic
183 <     * integration is only supported for simple rigid bodies.
184 <     */
185 <    for (mol = info_->beginMolecule(mi); mol != NULL;
186 <         mol = info_->nextMolecule(mi)) {
210 >    for (mol = info_->beginMolecule(i); mol != NULL;
211 >         mol = info_->nextMolecule(i)) {          
212        
213 <      for (integrableObject = mol->beginIntegrableObject(ii);
214 <           integrableObject != NULL;
215 <           integrableObject = mol->nextIntegrableObject(ii)) {    
213 >      // is this molecule restrained?    
214 >      GenericData* data = mol->getPropertyByName("Restraint");
215 >      
216 >      if (data != NULL) {
217          
218 <        eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
193 <        if(eof_test == NULL){
194 <          sprintf(painCave.errMsg,
195 <                  "RestReader Error: error in reading file %s\n"
196 <                  "natoms  = %d; index = %d\n"
197 <                  "error reading the line from the file.\n",
198 <                  inIdealFileName.c_str(), nTotObjs,
199 <                  integrableObject->getGlobalIndex() );
200 <          painCave.isFatal = 1;
201 <          simError();
202 <        }
218 >        // make sure we can reinterpret the generic data as restraint data:
219          
220 <        parseErr = parseIdealLine( read_buffer, integrableObject);
205 <        if( parseErr != NULL ){
206 <          strcpy( painCave.errMsg, parseErr );
207 <          painCave.isFatal = 1;
208 <          simError();
209 <        }
210 <      }
211 <    }
212 <    
213 <    // MPI Section of code..........
214 < #else //IS_MPI
215 <    
216 <    // first thing first, suspend fatalities.
217 <    painCave.isEventLoop = 1;
218 <    
219 <    int masterNode = 0;
220 <    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
221 <    
222 <    MPI_Status istatus;
223 <    int nCurObj;
224 <    int nitems;
225 <    int haveError;
226 <
227 <    nTotObjs = info_->getNGlobalIntegrableObjects();
228 <    haveError = 0;
229 <
230 <    if (worldRank == masterNode) {
231 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
232 <      if( eof_test == NULL ){
233 <        sprintf( painCave.errMsg,
234 <                 "Error reading 1st line of %s \n ",inIdealFileName.c_str());
235 <        painCave.isFatal = 1;
236 <        simError();
237 <      }
238 <      
239 <      nitems = atoi( read_buffer );
240 <      
241 <      // Check to see that the number of integrable objects in the
242 <      // intial configuration file is the same as derived from the
243 <      // meta-data file.
244 <      if( nTotObjs != nitems){
245 <        sprintf( painCave.errMsg,
246 <                 "RestraintReader Error. %s nIntegrable, %d, "
247 <                 "does not match the meta-data file's nIntegrable, %d.\n",
248 <                 inIdealFileName.c_str(), nTotObjs,
249 <                 info_->getNGlobalIntegrableObjects());
250 <        painCave.isFatal = 1;
251 <        simError();
252 <      }
253 <      
254 <      // skip over the comment line
255 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
256 <      if(eof_test == NULL){
257 <        sprintf( painCave.errMsg,
258 <                 "error in reading commment in %s\n", inIdealFileName.c_str());
259 <        painCave.isFatal = 1;
260 <        simError();
261 <      }
262 <
263 <      MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD);
264 <
265 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
266 <        int which_node = info_->getMolToProc(i);
220 >        RestraintData* restData= dynamic_cast<RestraintData*>(data);        
221          
222 <        if(which_node == masterNode){
269 <          //molecules belong to master node
222 >        if (restData != NULL) {
223            
224 <          mol = info_->getMoleculeByGlobalIndex(i);
224 >          // make sure we can reinterpet the restraint data as a
225 >          // pointer to a MolecularRestraint:
226            
227 <          if(mol == NULL) {
274 <            sprintf(painCave.errMsg,
275 <                    "RestReader Error: Molecule not found on node %d!\n",
276 <                    worldRank);
277 <            painCave.isFatal = 1;
278 <            simError();
279 <          }
227 >          MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
228            
229 <          for (integrableObject = mol->beginIntegrableObject(ii);
282 <               integrableObject != NULL;
283 <               integrableObject = mol->nextIntegrableObject(ii)){
229 >          if (mRest != NULL) {          
230              
231 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
231 >            // now we need to pack the stunt doubles for the reference
232 >            // structure:
233              
234 <            if(eof_test == NULL){
235 <              sprintf(painCave.errMsg,
236 <                      "RestReader Error: error in reading file %s\n"
237 <                      "natoms  = %d; index = %d\n"
291 <                      "error reading the line from the file.\n",
292 <                      inIdealFileName.c_str(), nTotObjs, i );
293 <              painCave.isFatal = 1;
294 <              simError();
295 <            }
296 <        
297 <            parseIdealLine(read_buffer, integrableObject);
298 <        
299 <          }
300 <
301 <        } else {
302 <          //molecule belongs to slave nodes
303 <          
304 <          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
305 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
306 <          
307 <          for(j = 0; j < nCurObj; j++){
234 >            std::vector<Vector3d> ref;
235 >            int count = 0;
236 >            RealType mass, mTot;
237 >            Vector3d COM(0.0);
238              
239 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
310 <            if(eof_test == NULL){
311 <              sprintf(painCave.errMsg,
312 <                      "RestReader Error: error in reading file %s\n"
313 <                      "natoms  = %d; index = %d\n"
314 <                      "error reading the line from the file.\n",
315 <                      inIdealFileName.c_str(), nTotObjs, i );
316 <              painCave.isFatal = 1;
317 <              simError();
318 <            }
239 >            mTot = 0.0;
240              
241 <            MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
242 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
322 <          }
323 <        }
324 <      }
325 <    } else {
326 <      //actions taken at slave nodes
327 <      MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD);
328 <
329 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
330 <        int which_node = info_->getMolToProc(i);
331 <
332 <        if(which_node == worldRank){
333 <          //molecule with global index i belongs to this processor
334 <          
335 <          mol = info_->getMoleculeByGlobalIndex(i);
336 <          
337 <          if(mol == NULL) {
338 <            sprintf(painCave.errMsg,
339 <                    "RestReader Error: molecule not found on node %d\n",
340 <                    worldRank);
341 <            painCave.isFatal = 1;
342 <            simError();
343 <          }
344 <          
345 <          nCurObj = mol->getNIntegrableObjects();
346 <          
347 <          MPI_Send(&nCurObj, 1, MPI_INT, masterNode,
348 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
349 <          
350 <          for (integrableObject = mol->beginIntegrableObject(ii);
351 <               integrableObject != NULL;
352 <               integrableObject = mol->nextIntegrableObject(ii)){
241 >            // loop over the stunt doubles in this molecule in the order we
242 >            // will be looping them in the restraint code:
243              
244 <            MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode,
245 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
246 <            
247 <            parseErr = parseIdealLine(read_buffer, integrableObject);
248 <            
249 <            if( parseErr != NULL ){
250 <              strcpy( painCave.errMsg, parseErr );
251 <              simError();
244 >            for (sd = mol->beginIntegrableObject(j); sd != NULL;
245 >                 sd = mol->nextIntegrableObject(j)) {
246 >              
247 >              // push back the reference positions of the stunt
248 >              // doubles from the *globally* sorted array of
249 >              // positions:
250 >              
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              }
257 +            COM /= mTot;
258 +            mRest->setReferenceStructure(ref, COM);        
259            }
260          }
261        }
262      }
367 – #endif
263    }
369 –  
370 –  char* RestReader::parseIdealLine(char* readLine, StuntDouble* sd){
371 –    
372 –    RealType pos[3];        // position place holders
373 –    RealType q[4];          // the quaternions
374 –    RealType RfromQ[3][3];  // the rotation matrix
375 –    RealType normalize;     // to normalize the reference unit vector
376 –    RealType uX, uY, uZ;    // reference unit vector place holders
377 –    RealType uselessToken;
378 –    StringTokenizer tokenizer(readLine);
379 –    int nTokens;
380 –    
381 –    nTokens = tokenizer.countTokens();
264  
265 <    if (nTokens < 14) {
266 <      sprintf(painCave.errMsg,
267 <              "RestReader Error: Not enough Tokens.\n");
386 <      painCave.isFatal = 1;
387 <      simError();
388 <    }
389 <    
390 <    std::string name = tokenizer.nextToken();
265 >  
266 >  
267 >  void RestReader::parseDumpLine(const std::string& line) {        
268  
269 <    if (name != sd->getType()) {
269 >    StringTokenizer tokenizer(line);
270 >    int nTokens;
271 >    
272 >    nTokens = tokenizer.countTokens();
273 >    
274 >    if (nTokens < 2) {
275 >      sprintf(painCave.errMsg,
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();
282 >
283 >  
284 >    if (index < 1116){
285 >      std::string type = tokenizer.nextToken();
286        
287 <      sprintf(painCave.errMsg,
288 <              "RestReader Error: Atom type [%s] in %s does not "
289 <              "match Atom Type [%s] in .md file.\n",
290 <              name.c_str(), inIdealFileName.c_str(),
291 <              sd->getType().c_str());
292 <      painCave.isFatal = 1;
293 <      simError();        
294 <    }
402 <    
403 <    pos[0] = tokenizer.nextTokenAsDouble();
404 <    pos[1] = tokenizer.nextTokenAsDouble();
405 <    pos[2] = tokenizer.nextTokenAsDouble();
406 <
407 <    // store the positions in the stuntdouble as generic data doubles
408 <    DoubleGenericData* refPosX = new DoubleGenericData();
409 <    refPosX->setID("refPosX");
410 <    refPosX->setData(pos[0]);
411 <    sd->addProperty(refPosX);
412 <
413 <    DoubleGenericData* refPosY = new DoubleGenericData();
414 <    refPosY->setID("refPosY");
415 <    refPosY->setData(pos[1]);
416 <    sd->addProperty(refPosY);
417 <    
418 <    DoubleGenericData* refPosZ = new DoubleGenericData();
419 <    refPosZ->setID("refPosZ");
420 <    refPosZ->setData(pos[2]);
421 <    sd->addProperty(refPosZ);
422 <
423 <    // we don't need the velocities
424 <    uselessToken = tokenizer.nextTokenAsDouble();
425 <    uselessToken = tokenizer.nextTokenAsDouble();
426 <    uselessToken = tokenizer.nextTokenAsDouble();
427 <    
428 <    if (sd->isDirectional()) {
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 <      q[0] = tokenizer.nextTokenAsDouble();
297 <      q[1] = tokenizer.nextTokenAsDouble();
432 <      q[2] = tokenizer.nextTokenAsDouble();
433 <      q[3] = tokenizer.nextTokenAsDouble();
296 >      bool exist = false;
297 >      int indexCount = 0;
298        
299 <      // now build the rotation matrix and find the unit vectors
300 <      RfromQ[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
301 <      RfromQ[0][1] = 2*(q[1]*q[2] + q[0]*q[3]);
302 <      RfromQ[0][2] = 2*(q[1]*q[3] - q[0]*q[2]);
303 <      RfromQ[1][0] = 2*(q[1]*q[2] - q[0]*q[3]);
440 <      RfromQ[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
441 <      RfromQ[1][2] = 2*(q[2]*q[3] + q[0]*q[1]);
442 <      RfromQ[2][0] = 2*(q[1]*q[3] + q[0]*q[2]);
443 <      RfromQ[2][1] = 2*(q[2]*q[3] - q[0]*q[1]);
444 <      RfromQ[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
299 >      while ( (!exist) && (indexCount < stuntDoubleIndex_.size())){
300 >        if (index == stuntDoubleIndex_[indexCount])
301 >          exist = true;
302 >        indexCount++;
303 >      }
304        
305 <      normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1]
447 <                       + RfromQ[2][2]*RfromQ[2][2]);
448 <      uX = RfromQ[2][0]/normalize;
449 <      uY = RfromQ[2][1]/normalize;
450 <      uZ = RfromQ[2][2]/normalize;
305 >      StuntDouble* integrableObject;
306        
307 <      // store reference unit vectors as generic data in the stuntdouble
308 <      DoubleGenericData* refVectorX = new DoubleGenericData();
309 <      refVectorX->setID("refVectorX");
310 <      refVectorX->setData(uX);
311 <      sd->addProperty(refVectorX);
312 <      
313 <      DoubleGenericData* refVectorY = new DoubleGenericData();
314 <      refVectorY->setID("refVectorY");
315 <      refVectorY->setData(uY);
316 <      sd->addProperty(refVectorY);
317 <      
318 <      DoubleGenericData* refVectorZ = new DoubleGenericData();
319 <      refVectorZ->setID("refVectorZ");
320 <      refVectorZ->setData(uZ);
321 <      sd->addProperty(refVectorZ);
307 >      if (exist){
308 >        
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 >        // 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 >          }
362 >        }
363 >      }
364      }
468 –    
469 –    // we don't need the angular velocities, so let's exit the line
470 –    return NULL;
365    }
366    
367 <  void RestReader::readZangle(){
367 >  void  RestReader::readStuntDoubles(std::istream& inputStream) {
368      
369 <    int i;
370 <    unsigned int j;
477 <    int isPresent;
369 >    inputStream.getline(buffer, bufferSize);
370 >    std::string line(buffer);
371      
372 <    Molecule* mol;
373 <    StuntDouble* integrableObject;
374 <    SimInfo::MoleculeIterator mi;
375 <    Molecule::IntegrableObjectIterator ii;
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 < #ifdef IS_MPI
380 <    int done, which_node, which_atom; // loop counter
486 <    int nProc;
487 <    MPI_Status istatus;
488 < #endif //is_mpi
489 <    
490 <    const int BUFFERSIZE = 2000; // size of the read buffer
491 <    int nTotObjs; // the number of atoms
492 <    char read_buffer[BUFFERSIZE]; //the line buffer for reading
493 <    
494 <    char *eof_test; // ptr to see when we reach the end of the file
495 <    char *parseErr;
496 <    
497 <    std::vector<StuntDouble*> vecParticles;
498 <    std::vector<RealType> tempZangs;
379 >    while(inputStream.getline(buffer, bufferSize)) {
380 >      line = buffer;
381        
382 <    inAngFileName = info_->getRestFileName();
383 <    
502 <    inAngFileName += "0";
503 <    
504 <    // open the omega value file for reading
505 < #ifdef IS_MPI
506 <    if (worldRank == 0) {
507 < #endif
508 <      isPresent = 1;
509 <      inAngFile = fopen(inAngFileName.c_str(), "r");
510 <      if(!inAngFile){
511 <        sprintf(painCave.errMsg,
512 <                "Restraints Warning: %s file is not present\n"
513 <                "\tAll omega values will be initialized to zero. If the\n"
514 <                "\tsimulation is starting from the idealCrystal.in reference\n"
515 <                "\tconfiguration, this is the desired action. If this is not\n"
516 <                "\tthe case, the energy calculations will be incorrect.\n",
517 <                inAngFileName.c_str());
518 <        painCave.severity = OOPSE_WARNING;
519 <        painCave.isFatal = 0;
520 <        simError();  
521 <        isPresent = 0;
382 >      if(line.find("</StuntDoubles>") != std::string::npos) {
383 >        break;
384        }
385        
386 < #ifdef IS_MPI
525 <      // let the other nodes know the status of the file search
526 <      MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
527 < #endif // is_mpi
528 <      
529 <      if (!isPresent) {
530 <        zeroZangle();
531 <        return;
532 <      }
533 <      
534 < #ifdef IS_MPI
386 >      parseDumpLine(line);
387      }
388      
389 <    // listen to node 0 to see if we should exit this function
390 <    if (worldRank != 0) {
391 <      MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
392 <      if (!isPresent) {
393 <        zeroZangle();
394 <        return;
395 <      }
389 >  }
390 >
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 >              "RestReader Error: Missing <FrameData>\n");
399 >      painCave.isFatal = 1;
400 >      simError();
401      }
402 <    
403 <    strcpy( checkPointMsg, "zAngle file opened successfully for reading." );
404 <    MPIcheckPoint();
405 < #endif
406 <    
407 < #ifndef IS_MPI
408 <    
409 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
410 <    if( eof_test == NULL ){
411 <      sprintf( painCave.errMsg,
412 <               "RestraintReader error: error reading 1st line of \"%s\"\n",
556 <               inAngFileName.c_str() );
557 <      painCave.isFatal = 1;
558 <      simError();
559 <    }
560 <    
561 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
562 <    while ( eof_test != NULL ) {
563 <      // check for and ignore blank lines
564 <      if ( read_buffer != NULL )
565 <        tempZangs.push_back( atof(read_buffer) );
566 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
567 <    }
568 <    
569 <    nTotObjs = info_->getNGlobalIntegrableObjects();
570 <    
571 <    if( nTotObjs != tempZangs.size() ){
572 <      sprintf( painCave.errMsg,
573 <               "RestraintReader zAngle reading error. %s nIntegrable, %d, "
574 <               "does not match the meta-data file's nIntegrable, %d.\n",
575 <               inAngFileName.c_str(), tempZangs.size(), nTotObjs );
576 <      painCave.isFatal = 1;
577 <      simError();
578 <    }
579 <    
580 <    // load the zAngles into the integrable objects
581 <    i = 0;
582 <    for (mol = info_->beginMolecule(mi); mol != NULL;
583 <         mol = info_->nextMolecule(mi)) {
402 >
403 >    // restraints don't care about frame data (unless we need to wrap
404 >    // coordinates, but we'll worry about that later), so
405 >    // we'll just scan ahead until the end of the frame data:
406 >
407 >    while(inputStream.getline(buffer, bufferSize)) {
408 >      line = buffer;
409 >      
410 >      if(line.find("</FrameData>") != std::string::npos) {
411 >        break;
412 >      }
413        
585 –      for (integrableObject = mol->beginIntegrableObject(ii);
586 –           integrableObject != NULL;
587 –           integrableObject = mol->nextIntegrableObject(ii)) {    
588 –        
589 –        integrableObject->setZangle(tempZangs[i]);
590 –        i++;
591 –      }
414      }
415      
594 –    // MPI Section of code..........
595 – #else //IS_MPI
596 –    
597 –    // first thing first, suspend fatalities.
598 –    painCave.isEventLoop = 1;
599 –
600 –    int masterNode = 0;
601 –    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
602 –    int haveError;
603 –    int intObjIndex;    
604 –    int intObjIndexTransfer;    
605 –
606 –    int nCurObj;
607 –    RealType angleTranfer;
608 –    
609 –    nTotObjs = info_->getNGlobalIntegrableObjects();
610 –    haveError = 0;
611 –
612 –    if (worldRank == masterNode) {
613 –      
614 –      eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
615 –      if( eof_test == NULL ){
616 –        sprintf( painCave.errMsg,
617 –                 "Error reading 1st line of %s \n ",inAngFileName.c_str());
618 –        haveError = 1;
619 –        simError();
620 –      }
621 –      
622 –      // let node 0 load the temporary angle vector
623 –      eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
624 –      while ( eof_test != NULL ) {
625 –        // check for and ignore blank lines
626 –        if ( read_buffer != NULL )
627 –          tempZangs.push_back( atof(read_buffer) );
628 –        eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
629 –      }
630 –
631 –      // Check to see that the number of integrable objects in the
632 –      // intial configuration file is the same as derived from the
633 –      // meta-data file.
634 –      if( nTotObjs != tempZangs.size() ){
635 –        sprintf( painCave.errMsg,
636 –                 "RestraintReader zAngle reading Error. %s nIntegrable, %d, "
637 –                 "does not match the meta-data file's nIntegrable, %d.\n",
638 –                 inAngFileName.c_str(), tempZangs.size(), nTotObjs);
639 –        haveError= 1;
640 –        simError();
641 –      }
642 –      
643 –      // At this point, node 0 has a tempZangs vector completed, and
644 –      // everyone else has nada
645 –      
646 –      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
647 –        // Get the Node number which has this atom
648 –        which_node = info_->getMolToProc(i);
649 –        
650 –        if (which_node == masterNode) {
651 –          mol = info_->getMoleculeByGlobalIndex(i);
652 –
653 –          if(mol == NULL) {
654 –            strcpy(painCave.errMsg, "Molecule not found on node 0!");
655 –            haveError = 1;
656 –            simError();
657 –          }
658 –
659 –          for (integrableObject = mol->beginIntegrableObject(ii);
660 –               integrableObject != NULL;
661 –               integrableObject = mol->nextIntegrableObject(ii)){
662 –            intObjIndex = integrableObject->getGlobalIndex();
663 –            integrableObject->setZangle(tempZangs[intObjIndex]);
664 –          }    
665 –          
666 –        } else {
667 –          // I am MASTER OF THE UNIVERSE, but I don't own this molecule
668 –          // listen for the number of integrableObjects in the molecule
669 –          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
670 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
671 –          
672 –          for(j=0; j < nCurObj; j++){          
673 –            // listen for which integrableObject we need to send the value for
674 –            MPI_Recv(&intObjIndexTransfer, 1, MPI_INT, which_node,
675 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
676 –            angleTransfer = tempZangs[intObjIndexTransfer];
677 –            // send the value to the node so it can initialize the object
678 –            MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node,
679 –                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);              
680 –          }
681 –        }
682 –      }
683 –    } else {
684 –      // I am SLAVE TO THE MASTER
685 –      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
686 –        int which_node = info_->getMolToProc(i);
687 –
688 –        if (which_node == worldRank) {
689 –          
690 –          // BUT I OWN THIS MOLECULE!!!
691 –          
692 –          mol = info_->getMoleculeByGlobalIndex(i);
693 –
694 –          if(mol == NULL) {
695 –            sprintf(painCave.errMsg,
696 –                    "RestReader Error: molecule not found on node %d\n",
697 –                    worldRank);
698 –            painCave.isFatal = 1;
699 –            simError();
700 –          }
701 –
702 –          nCurObj = mol->getNIntegrableObjects();
703 –          // send the number of integrableObjects in the molecule
704 –          MPI_Send(&nCurObj, 1, MPI_INT, 0,
705 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
706 –          
707 –          for (integrableObject = mol->beginIntegrableObject(ii);
708 –               integrableObject != NULL;
709 –               integrableObject = mol->nextIntegrableObject(ii)){
710 –            intObjIndexTransfer = integrableObject->getGlobalIndex();
711 –            // send the global index of the integrableObject
712 –            MPI_Send(&intObjIndexTransfer, 1, MPI_INT, 0,
713 –                     TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
714 –            // listen for the value we want to set locally
715 –            MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0,
716 –                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
717 –
718 –            integrableObject->setZangle(angleTransfer);
719 –          }    
720 –        }
721 –      }
722 –    }
723 – #endif
416    }
725 –  
726 –  void RestReader :: zeroZangle(){
727 –    
728 –    int i;
729 –    unsigned int j;
730 –    int nTotObjs; // the number of atoms
731 –    
732 –    Molecule* mol;
733 –    StuntDouble* integrableObject;
734 –    SimInfo::MoleculeIterator mi;
735 –    Molecule::IntegrableObjectIterator ii;
736 –    
737 –    std::vector<StuntDouble*> vecParticles;
738 –    
739 – #ifndef IS_MPI
740 –    // set all zAngles to 0.0
741 –    for (mol = info_->beginMolecule(mi); mol != NULL;
742 –         mol = info_->nextMolecule(mi))
743 –      
744 –      for (integrableObject = mol->beginIntegrableObject(ii);
745 –           integrableObject != NULL;
746 –           integrableObject = mol->nextIntegrableObject(ii))    
747 –        integrableObject->setZangle( 0.0 );
748 –    
749 –    
750 –    // MPI Section of code..........
751 – #else //IS_MPI
752 –    
753 –    // first thing first, suspend fatalities.
754 –    painCave.isEventLoop = 1;
755 –    
756 –    int masterNode = 0;
757 –    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
758 –    int haveError;
759 –    int which_node;
760 –    
761 –    MPI_Status istatus;
762 –    
763 –    int nCurObj;
764 –    RealType angleTranfer;
765 –    
766 –    nTotObjs = info_->getNGlobalIntegrableObjects();
767 –    haveError = 0;
768 –    if (worldRank == masterNode) {
417  
418 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
419 <        // Get the Node number which has this atom
772 <        which_node = info_->getMolToProc(i);
773 <        
774 <        // let's let node 0 pass out constant values to all the processors
775 <        if (worldRank == masterNode) {
776 <          mol = info_->getMoleculeByGlobalIndex(i);
777 <          
778 <          if(mol == NULL) {
779 <            strcpy(painCave.errMsg, "Molecule not found on node 0!");
780 <            haveError = 1;
781 <            simError();
782 <          }
783 <          
784 <          for (integrableObject = mol->beginIntegrableObject(ii);
785 <               integrableObject != NULL;
786 <               integrableObject = mol->nextIntegrableObject(ii)){
787 <            
788 <            integrableObject->setZangle( 0.0 );
789 <            
790 <          }
791 <          
792 <        } else {
793 <          // I am MASTER OF THE UNIVERSE, but I don't own this molecule
794 <          
795 <          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
796 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
797 <          
798 <          for(j=0; j < nCurObj; j++){            
799 <            angleTransfer = 0.0;
800 <            MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node,
801 <                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);              
802 <            
803 <          }
804 <        }
805 <      }
806 <    } else {
807 <      // I am SLAVE TO THE MASTER
808 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
809 <        int which_node = info_->getMolToProc(i);
810 <        
811 <        if (which_node == worldRank) {
812 <          
813 <          // BUT I OWN THIS MOLECULE!!!
814 <          mol = info_->getMoleculeByGlobalIndex(i);
815 <          
816 <          if(mol == NULL) {
817 <            sprintf(painCave.errMsg,
818 <                    "RestReader Error: molecule not found on node %d\n",
819 <                    worldRank);
820 <            painCave.isFatal = 1;
821 <            simError();
822 <          }
823 <          
824 <          nCurObj = mol->getNIntegrableObjects();
825 <          
826 <          MPI_Send(&nCurObj, 1, MPI_INT, 0,
827 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
828 <          
829 <          for (integrableObject = mol->beginIntegrableObject(ii);
830 <               integrableObject != NULL;
831 <               integrableObject = mol->nextIntegrableObject(ii)){
832 <            
833 <            MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0,
834 <                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
835 <            vecParticles[j]->setZangle(angleTransfer);
836 <          }    
837 <        }
838 <      }
839 <    }
840 < #endif
841 <  }
842 <  
843 < } // end namespace oopse
418 >  
419 > }//end namespace OpenMD

Diff Legend

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