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

Comparing trunk/src/io/RestReader.cpp (file contents):
Revision 1241 by gezelter, Fri Apr 25 15:14:47 2008 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 1 | Line 1
1 < /*
2 < * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 < *
1 > /*
2 > * Copyright (c) 2005 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 < */
32 <
33 < #define _LARGEFILE_SOURCE64
34 < #define _FILE_OFFSET_BITS 64
35 <
36 < #include <sys/types.h>
37 < #include <sys/stat.h>
38 <
39 < #include <iostream>
40 < #include <math.h>
41 <
42 < #include <stdio.h>
52 < #include <stdlib.h>
53 < #include <string.h>
54 <
55 < #include "primitives/Molecule.hpp"
56 < #include "utils/MemoryUtils.hpp"
57 < #include "utils/StringTokenizer.hpp"
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 > #include "io/DumpReader.hpp"
43   #include "io/RestReader.hpp"
44 < #include "utils/simError.h"
44 > #include "primitives/Molecule.hpp"
45 > #include "restraints/ObjectRestraint.hpp"
46 > #include "restraints/MolecularRestraint.hpp"
47  
48 < #ifdef IS_MPI
49 < #include <mpi.h>
50 < #define TAKE_THIS_TAG_CHAR 0
64 < #define TAKE_THIS_TAG_INT 1
65 < #define TAKE_THIS_TAG_DOUBLE 2
66 < #endif // is_mpi
48 > namespace OpenMD {
49 >      
50 >  void RestReader::readReferenceStructure() {
51  
52 < namespace oopse {
53 <  
54 <  RestReader::RestReader( SimInfo* info ) : info_(info){
55 <    
56 <    idealName = "idealCrystal.in";
52 >    // some of this comes directly from DumpReader.
53 >
54 >    if (!isScanned_)
55 >      scanFile();
56 >        
57 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
58      
59 < #ifdef IS_MPI
60 <    if (worldRank == 0) {
61 < #endif
59 >    if (storageLayout & DataStorage::dslPosition) {
60 >      needPos_ = true;
61 >    } else {
62 >      needPos_ = false;
63 >    }
64 >    
65 >    needVel_ = false;
66 >    
67 >    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
68 >      needQuaternion_ = true;
69 >    } else {
70 >      needQuaternion_ = false;
71 >    }
72  
73 <      inIdealFile = new std::ifstream(idealName.c_str());
73 >    needAngMom_ = false;
74  
75 <      if(inIdealFile->fail()){
76 <        sprintf(painCave.errMsg,
77 <                "RestReader: Cannot open file: %s\n",
83 <                idealName.c_str());
84 <        painCave.isFatal = 1;
85 <        simError();
86 <      }
87 <      
88 < #ifdef IS_MPI
89 <    }
90 <    strcpy( checkPointMsg,
91 <            "File \"idealCrystal.in\" opened successfully for reading." );
92 <    errorCheckPoint();
93 < #endif
75 >    // We need temporary storage to keep track of all StuntDouble positions
76 >    // in case some of the restraints are molecular (i.e. if they use
77 >    // multiple SD positions to determine restrained orientations or positions:
78  
79 <    return;  
80 <  }
97 <  
98 <  RestReader :: ~RestReader( ){
99 < #ifdef IS_MPI
100 <    if (worldRank == 0) {
101 < #endif
79 >    all_pos_.clear();
80 >    all_pos_.resize(info_->getNGlobalIntegrableObjects() );
81  
103 –      delete inIdealFile;
104 –      delete inAngFile;
82  
83 < #ifdef IS_MPI
84 <    }
85 <    strcpy( checkPointMsg,
86 <            "File idealCrystal.in (and .zang0 if present) closed successfully." );
87 <    errorCheckPoint();
88 < #endif
83 >    // Restraint files are just standard dump files, but with the reference
84 >    // structure stored in the first frame (frame 0).
85 >    // RestReader overloads readSet and explicitly handles all of the
86 >    // ObjectRestraints in that method:
87 >
88 >    readSet(0);
89      
90 <    return;
91 <  }
92 <  
93 <  
94 <  void RestReader :: readIdealCrystal(){
118 <        
119 < #ifdef IS_MPI
120 <    int which_node;
121 <    int i, j;
122 < #endif //is_mpi
123 <    
124 <    const int BUFFERSIZE = 2000; // size of the read buffer
125 <    int nTotObjs; // the number of atoms
126 <    char read_buffer[BUFFERSIZE]; //the line buffer for reading
127 <    
128 <    char *parseErr;
129 <    
130 <    std::vector<StuntDouble*> integrableObjects;
131 <    
90 >    // all ObjectRestraints have been handled, now we have to worry about
91 >    // molecular restraints:
92 >
93 >    SimInfo::MoleculeIterator i;
94 >    Molecule::IntegrableObjectIterator j;
95      Molecule* mol;
96 <    StuntDouble* integrableObject;
134 <    SimInfo::MoleculeIterator mi;
135 <    Molecule::IntegrableObjectIterator ii;
136 <    
137 < #ifndef IS_MPI
138 <    
139 <    inIdealFile->getline(read_buffer, sizeof(read_buffer));
96 >    StuntDouble* sd;
97  
98 <    if( inIdealFile->eof() ){
99 <      sprintf( painCave.errMsg,
100 <               "RestraintReader error: error reading 1st line of \"%s\"\n",
144 <               idealName.c_str() );
145 <      painCave.isFatal = 1;
146 <      simError();
147 <    }
148 <    
149 <    nTotObjs = atoi( read_buffer );
150 <    
151 <    if( nTotObjs != info_->getNGlobalIntegrableObjects() ){
152 <      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 <    }
161 <    
162 <    // skip over the comment line
163 <    inIdealFile->getline(read_buffer, sizeof(read_buffer));
98 >    // no need to worry about parallel molecules, as molecules are not
99 >    // split across processor boundaries.  Just loop over all molecules
100 >    // we know about:
101  
102 <    if( inIdealFile->eof() ){
103 <      sprintf( painCave.errMsg,
104 <               "error in reading commment in %s\n",
105 <               idealName.c_str());
106 <      painCave.isFatal = 1;
170 <      simError();
171 <    }
172 <    
173 <    // parse the ideal crystal lines
174 <    /*
175 <     * Note: we assume that there is a one-to-one correspondence between
176 <     * 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)) {
102 >    for (mol = info_->beginMolecule(i); mol != NULL;
103 >         mol = info_->nextMolecule(i)) {          
104 >
105 >      // is this molecule restrained?    
106 >      GenericData* data = mol->getPropertyByName("Restraint");
107        
108 <      for (integrableObject = mol->beginIntegrableObject(ii);
183 <           integrableObject != NULL;
184 <           integrableObject = mol->nextIntegrableObject(ii)) {    
108 >      if (data != NULL) {
109  
110 <        inIdealFile->getline(read_buffer, sizeof(read_buffer));
110 >        // make sure we can reinterpret the generic data as restraint data:
111  
112 <        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;
112 >        RestraintData* restData= dynamic_cast<RestraintData*>(data);        
113  
114 <    nTotObjs = info_->getNGlobalIntegrableObjects();
222 <    haveError = 0;
114 >        if (restData != NULL) {
115  
116 <    if (worldRank == masterNode) {
117 <      inIdealFile->getline(read_buffer, sizeof(read_buffer));
116 >          // make sure we can reinterpet the restraint data as a
117 >          // pointer to a MolecularRestraint:
118  
119 <      if( inIdealFile->eof() ){
228 <        sprintf( painCave.errMsg,
229 <                 "Error reading 1st line of %s \n ",idealName.c_str());
230 <        painCave.isFatal = 1;
231 <        simError();
232 <      }
233 <      
234 <      nitems = atoi( read_buffer );
235 <      
236 <      // 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));
119 >          MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
120  
121 <      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 <      }
121 >          if (mRest != NULL) {          
122  
123 <      MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD);
123 >            // now we need to pack the stunt doubles for the reference
124 >            // structure:
125  
126 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
127 <        which_node = info_->getMolToProc(i);
128 <        
129 <        if(which_node == masterNode){
265 <          //molecules belong to master node
266 <          
267 <          mol = info_->getMoleculeByGlobalIndex(i);
268 <          
269 <          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 <          }
276 <          
277 <          for (integrableObject = mol->beginIntegrableObject(ii);
278 <               integrableObject != NULL;
279 <               integrableObject = mol->nextIntegrableObject(ii)){
126 >            std::vector<Vector3d> ref;
127 >            int count = 0;
128 >            RealType mass, mTot;
129 >            Vector3d COM(0.0);
130              
131 <            inIdealFile->getline(read_buffer, sizeof(read_buffer));
131 >            mTot = 0.0;
132              
133 <            if( inIdealFile->eof() ){
134 <              sprintf(painCave.errMsg,
135 <                      "RestReader Error: error in reading file %s\n"
136 <                      "natoms  = %d; index = %d\n"
137 <                      "error reading the line from the file.\n",
138 <                      idealName.c_str(), nTotObjs, i );
139 <              painCave.isFatal = 1;
140 <              simError();
133 >            // loop over the stunt doubles in this molecule in the order we
134 >            // will be looping them in the restraint code:
135 >            
136 >            for (sd = mol->beginIntegrableObject(j); sd != NULL;
137 >                 sd = mol->nextIntegrableObject(j)) {
138 >              
139 >              // push back the reference positions of the stunt
140 >              // doubles from the *globally* sorted array of
141 >              // positions:
142 >              
143 >              ref.push_back( all_pos_[sd->getGlobalIndex()] );
144 >              
145 >              mass = sd->getMass();
146 >              
147 >              COM = COM + mass * ref[count];
148 >              mTot = mTot + mass;
149 >              count = count + 1;
150              }
292 –        
293 –            parseIdealLine(read_buffer, integrableObject);
294 –        
295 –          }
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 –            }
151              
152 <            MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
317 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
318 <          }
319 <        }
320 <      }
321 <    } else {
322 <      //actions taken at slave nodes
323 <      MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode, MPI_COMM_WORLD);
152 >            COM /= mTot;
153  
154 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
326 <        int which_node = info_->getMolToProc(i);
154 >            mRest->setReferenceStructure(ref, COM);        
155  
328 –        if(which_node == worldRank){
329 –          //molecule with global index i belongs to this processor
330 –          
331 –          mol = info_->getMoleculeByGlobalIndex(i);
332 –          
333 –          if(mol == NULL) {
334 –            sprintf(painCave.errMsg,
335 –                    "RestReader Error: molecule not found on node %d\n",
336 –                    worldRank);
337 –            painCave.isFatal = 1;
338 –            simError();
156            }
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);
354 –            
355 –            if( parseErr != NULL ){
356 –              strcpy( painCave.errMsg, parseErr );
357 –              simError();
358 –            }
359 –          }
157          }
158        }
159      }
363 – #endif
160    }
365 –  
366 –  char* RestReader::parseIdealLine(char* readLine, StuntDouble* sd){
367 –    
368 –    RealType pos[3];        // position place holders
369 –    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;
376 –    
377 –    nTokens = tokenizer.countTokens();
161  
162 <    if (nTokens < 14) {
163 <      sprintf(painCave.errMsg,
164 <              "RestReader Error: Not enough Tokens.\n");
165 <      painCave.isFatal = 1;
166 <      simError();
167 <    }
168 <    
169 <    std::string name = tokenizer.nextToken();
162 >  
163 >  
164 >  void RestReader::parseDumpLine(const std::string& line) {        
165 >    StringTokenizer tokenizer(line);
166 >    int nTokens;
167 >    
168 >    nTokens = tokenizer.countTokens();
169 >    
170 >    if (nTokens < 2) {
171 >      sprintf(painCave.errMsg,
172 >              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
173 >      painCave.isFatal = 1;
174 >      simError();
175 >    }
176  
177 <    if (name != sd->getType()) {
178 <      
179 <      sprintf(painCave.errMsg,
391 <              "RestReader Error: Atom type [%s] in %s does not "
392 <              "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();
177 >    int index = tokenizer.nextTokenAsInt();
178 >
179 >    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
180  
181 <    // store the positions in the stuntdouble as generic data doubles
182 <    DoubleGenericData* refPosX = new DoubleGenericData();
183 <    refPosX->setID("refPosX");
406 <    refPosX->setData(pos[0]);
407 <    sd->addProperty(refPosX);
181 >    if (integrableObject == NULL) {
182 >      return;
183 >    }  
184  
185 <    DoubleGenericData* refPosY = new DoubleGenericData();
186 <    refPosY->setID("refPosY");
187 <    refPosY->setData(pos[1]);
188 <    sd->addProperty(refPosY);
189 <    
190 <    DoubleGenericData* refPosZ = new DoubleGenericData();
191 <    refPosZ->setID("refPosZ");
192 <    refPosZ->setData(pos[2]);
193 <    sd->addProperty(refPosZ);
194 <
195 <    // we don't need the velocities
196 <    uselessToken = tokenizer.nextTokenAsDouble();
197 <    uselessToken = tokenizer.nextTokenAsDouble();
198 <    uselessToken = tokenizer.nextTokenAsDouble();
199 <    
200 <    if (sd->isDirectional()) {
201 <      
202 <      q[0] = tokenizer.nextTokenAsDouble();
203 <      q[1] = tokenizer.nextTokenAsDouble();
204 <      q[2] = tokenizer.nextTokenAsDouble();
205 <      q[3] = tokenizer.nextTokenAsDouble();
206 <      
207 <      // now build the rotation matrix and find the unit vectors
208 <      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);
463 <    }
464 <    
465 <    // we don't need the angular velocities, so let's exit the line
466 <    return NULL;
467 <  }
468 <  
469 <  void RestReader::readZangle(){
470 <    
471 <    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 <      }
185 >    std::string type = tokenizer.nextToken();
186 >    int size = type.size();
187 >    Vector3d pos;
188 >    Vector3d vel;
189 >    Quat4d q;
190 >    Vector3d ji;
191 >    Vector3d force;
192 >    Vector3d torque;
193 >          
194 >    for(int i = 0; i < size; ++i) {
195 >      switch(type[i]) {
196 >        
197 >        case 'p': {
198 >            pos[0] = tokenizer.nextTokenAsDouble();
199 >            pos[1] = tokenizer.nextTokenAsDouble();
200 >            pos[2] = tokenizer.nextTokenAsDouble();
201 >            break;
202 >        }
203 >        case 'v' : {
204 >            vel[0] = tokenizer.nextTokenAsDouble();
205 >            vel[1] = tokenizer.nextTokenAsDouble();
206 >            vel[2] = tokenizer.nextTokenAsDouble();
207 >            break;
208 >        }
209  
210 <    }
211 <    
212 <    // listen to node 0 to see if we should exit this function
213 <    if (worldRank != 0) {
214 <      MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
215 <      if (!isPresent) {
216 <        zeroZangle();
217 <        return;
218 <      }
219 <    }
220 <    
221 <    strcpy( checkPointMsg, "zAngle file opened successfully for reading." );
222 <    errorCheckPoint();
223 < #endif
224 <    
225 < #ifndef IS_MPI
226 <    
227 <    // read the first line and die if there is a failure
228 <    inAngFile->getline(read_buffer, sizeof(read_buffer));
210 >        case 'q' : {
211 >           if (integrableObject->isDirectional()) {
212 >              
213 >             q[0] = tokenizer.nextTokenAsDouble();
214 >             q[1] = tokenizer.nextTokenAsDouble();
215 >             q[2] = tokenizer.nextTokenAsDouble();
216 >             q[3] = tokenizer.nextTokenAsDouble();
217 >              
218 >             RealType qlen = q.length();
219 >             if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
220 >                
221 >               sprintf(painCave.errMsg,
222 >                       "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
223 >               painCave.isFatal = 1;
224 >               simError();
225 >                
226 >             }  
227 >              
228 >             q.normalize();
229 >           }            
230 >           break;
231 >        }  
232 >        case 'j' : {
233 >          if (integrableObject->isDirectional()) {
234 >             ji[0] = tokenizer.nextTokenAsDouble();
235 >             ji[1] = tokenizer.nextTokenAsDouble();
236 >             ji[2] = tokenizer.nextTokenAsDouble();
237 >          }
238 >          break;
239 >        }  
240 >        case 'f': {
241 >          force[0] = tokenizer.nextTokenAsDouble();
242 >          force[1] = tokenizer.nextTokenAsDouble();
243 >          force[2] = tokenizer.nextTokenAsDouble();          
244  
245 <    if( inAngFile->eof() ){
246 <      sprintf( painCave.errMsg,
247 <               "RestraintReader error: error reading 1st line of \"%s\"\n",
248 <               angFile.c_str() );
249 <      painCave.isFatal = 1;
250 <      simError();
560 <    }
245 >          break;
246 >        }
247 >        case 't' : {
248 >           torque[0] = tokenizer.nextTokenAsDouble();
249 >           torque[1] = tokenizer.nextTokenAsDouble();
250 >           torque[2] = tokenizer.nextTokenAsDouble();          
251  
252 <    // read the file and load the values into a vector
253 <    inAngFile->getline(read_buffer, sizeof(read_buffer));
254 <    
255 <    while ( !inAngFile->eof() ) {
256 <      // check for and ignore blank lines
257 <      if ( read_buffer != NULL )
258 <        tempZangs.push_back( atof(read_buffer) );
259 <
260 <      inAngFile->getline(read_buffer, sizeof(read_buffer));
261 <    }
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)) {
590 <      
591 <      for (integrableObject = mol->beginIntegrableObject(ii);
592 <           integrableObject != NULL;
593 <           integrableObject = mol->nextIntegrableObject(ii)) {    
594 <        
595 <        integrableObject->setZangle(tempZangs[i]);
596 <        i++;
252 >           break;
253 >        }
254 >        default: {
255 >               sprintf(painCave.errMsg,
256 >                       "DumpReader Error: %s is an unrecognized type\n", type.c_str());
257 >               painCave.isFatal = 1;
258 >               simError();
259 >          break;  
260 >        }
261 >          
262        }
263      }
264 +    
265 +    // keep the position in case we need it for a molecular restraint:
266      
267 <    // MPI Section of code..........
601 < #else //IS_MPI
602 <    
603 <    // first thing first, suspend fatalities.
604 <    painCave.isEventLoop = 1;
267 >    all_pos_[index] = pos;
268  
269 <    int masterNode = 0;
269 >    // is this io restrained?    
270 >    GenericData* data = integrableObject->getPropertyByName("Restraint");
271 >    ObjectRestraint* oRest;
272  
273 <    int haveError;
274 <    int intObjIndex;    
275 <    int intObjIndexTransfer;    
276 <
277 <    int j;
278 <    int nCurObj;
279 <    RealType angleTransfer;
280 <    
281 <    nTotObjs = info_->getNGlobalIntegrableObjects();
282 <    haveError = 0;
283 <
284 <    if (worldRank == masterNode) {
285 <      
286 <      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();
273 >    if (data != NULL) {
274 >      // make sure we can reinterpret the generic data as restraint data:
275 >      RestraintData* restData= dynamic_cast<RestraintData*>(data);        
276 >      if (restData != NULL) {
277 >        // make sure we can reinterpet the restraint data as a pointer to
278 >        // an ObjectRestraint:
279 >        oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
280 >        if (oRest != NULL) {          
281 >          if (integrableObject->isDirectional()) {
282 >            oRest->setReferenceStructure(pos, q.toRotationMatrix3());
283 >          } else {                          
284 >            oRest->setReferenceStructure(pos);
285 >          }
286 >        }
287        }
288 <      
630 <      // let the master node read the file and load the temporary angle vector
631 <      inAngFile->getline(read_buffer, sizeof(read_buffer));
288 >    }
289  
290 <      while ( !inAngFile->eof() ) {
291 <        // check for and ignore blank lines
635 <        if ( read_buffer != NULL )
636 <          tempZangs.push_back( atof(read_buffer) );
290 >  }
291 >  
292  
638 –        inAngFile->getline(read_buffer, sizeof(read_buffer));
639 –      }
293  
294 <      // Check to see that the number of integrable objects in the
295 <      // intial configuration file is the same as derived from the
296 <      // 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);
294 >  void RestReader::readFrameProperties(std::istream& inputStream) {
295 >    inputStream.getline(buffer, bufferSize);
296 >    std::string line(buffer);
297  
298 <          if(mol == NULL) {
299 <            strcpy(painCave.errMsg, "Molecule not found on node 0!");
300 <            haveError = 1;
301 <            simError();
302 <          }
298 >    if (line.find("<FrameData>") == std::string::npos) {
299 >      sprintf(painCave.errMsg,
300 >              "DumpReader Error: Missing <FrameData>\n");
301 >      painCave.isFatal = 1;
302 >      simError();
303 >    }
304  
305 <          for (integrableObject = mol->beginIntegrableObject(ii);
306 <               integrableObject != NULL;
307 <               integrableObject = mol->nextIntegrableObject(ii)){
674 <            intObjIndex = integrableObject->getGlobalIntegrableObjectIndex();
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);
305 >    // restraints don't care about frame data (unless we need to wrap
306 >    // coordinates, but we'll worry about that later), so
307 >    // we'll just scan ahead until the end of the frame data:
308  
309 <        if (which_node == worldRank) {
310 <          
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->getGlobalIntegrableObjectIndex();
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
736 <  }
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))
309 >    while(inputStream.getline(buffer, bufferSize)) {
310 >      line = buffer;
311        
312 <      for (integrableObject = mol->beginIntegrableObject(ii);
313 <           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;
767 <
768 <    for (int i=0 ; i < info_->getNGlobalMolecules(); i++) {
769 <
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 <        }
312 >      if(line.find("</FrameData>") != std::string::npos) {
313 >        break;
314        }
315 +      
316      }
317 <
795 < #endif
317 >
318    }
319 <  
320 < } // end namespace oopse
319 >
320 >  
321 > }//end namespace OpenMD

Diff Legend

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