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 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 1360 by cli2, Mon Sep 7 16:31:51 2009 UTC

# Line 1 | Line 1
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
19 < *    notice, this list of conditions and the following disclaimer.
20 < *
21 < * 3. Redistributions in binary form must reproduce the above copyright
22 < *    notice, this list of conditions and the following disclaimer in the
23 < *    documentation and/or other materials provided with the
24 < *    distribution.
25 < *
26 < * This software is provided "AS IS," without a warranty of any
27 < * kind. All express or implied conditions, representations and
28 < * warranties, including any implied warranty of merchantability,
29 < * fitness for a particular purpose or non-infringement, are hereby
30 < * excluded.  The University of Notre Dame and its licensors shall not
31 < * be liable for any damages suffered by licensee as a result of
32 < * using, modifying or distributing the software or its
33 < * derivatives. In no event will the University of Notre Dame or its
34 < * licensors be liable for any lost revenue, profit or data, or for
35 < * direct, indirect, special, consequential, incidental or punitive
36 < * damages, however caused and regardless of the theory of liability,
37 < * arising out of the use of or inability to use software, even if the
38 < * University of Notre Dame has been advised of the possibility of
39 < * such damages.
40 < */
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
19 > *    notice, this list of conditions and the following disclaimer.
20 > *
21 > * 3. Redistributions in binary form must reproduce the above copyright
22 > *    notice, this list of conditions and the following disclaimer in the
23 > *    documentation and/or other materials provided with the
24 > *    distribution.
25 > *
26 > * This software is provided "AS IS," without a warranty of any
27 > * kind. All express or implied conditions, representations and
28 > * warranties, including any implied warranty of merchantability,
29 > * fitness for a particular purpose or non-infringement, are hereby
30 > * excluded.  The University of Notre Dame and its licensors shall not
31 > * be liable for any damages suffered by licensee as a result of
32 > * using, modifying or distributing the software or its
33 > * derivatives. In no event will the University of Notre Dame or its
34 > * licensors be liable for any lost revenue, profit or data, or for
35 > * direct, indirect, special, consequential, incidental or punitive
36 > * damages, however caused and regardless of the theory of liability,
37 > * arising out of the use of or inability to use software, even if the
38 > * University of Notre Dame has been advised of the possibility of
39 > * such damages.
40 > */
41 >  
42 > #include "io/DumpReader.hpp"
43 > #include "io/RestReader.hpp"
44 > #include "primitives/Molecule.hpp"
45 > #include "restraints/ObjectRestraint.hpp"
46 > #include "restraints/MolecularRestraint.hpp"
47  
48 < #define _LARGEFILE_SOURCE64
49 < #define _FILE_OFFSET_BITS 64
48 > namespace oopse {
49 >      
50 >  void RestReader::readReferenceStructure() {
51  
52 < #include <sys/types.h>
46 < #include <sys/stat.h>
52 >    // some of this comes directly from DumpReader.
53  
54 < #include <iostream>
55 < #include <math.h>
54 >    if (!isScanned_)
55 >      scanFile();
56 >        
57 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
58 >    
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 < #include <stdio.h>
52 < #include <stdlib.h>
53 < #include <string.h>
73 >    needAngMom_ = false;
74  
75 < #include "primitives/Molecule.hpp"
76 < #include "utils/MemoryUtils.hpp"
77 < #include "utils/StringTokenizer.hpp"
58 < #include "io/RestReader.hpp"
59 < #include "utils/simError.h"
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 < #ifdef IS_MPI
80 < #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
79 >    all_pos_.clear();
80 >    all_pos_.resize(info_->getNGlobalIntegrableObjects() );
81  
68 – namespace oopse {
69 –  
70 –  RestReader::RestReader( SimInfo* info ) : info_(info){
71 –    
72 –    idealName = "idealCrystal.in";
73 –    
74 –    isScanned = false;
75 –    
76 – #ifdef IS_MPI
77 –    if (worldRank == 0) {
78 – #endif
79 –      
80 –      inIdealFile = fopen(idealName, "r");
81 –      if(inIdealFile == NULL){
82 –        sprintf(painCave.errMsg,
83 –                "RestReader: Cannot open file: %s\n", idealName);
84 –        painCave.isFatal = 1;
85 –        simError();
86 –      }
87 –      
88 –      inIdealFileName = idealName;
89 – #ifdef IS_MPI
90 –    }
91 –    strcpy( checkPointMsg,
92 –            "File \"idealCrystal.in\" opened successfully for reading." );
93 –    MPIcheckPoint();
94 – #endif
82  
83 <    return;  
84 <  }
85 <  
86 <  RestReader :: ~RestReader( ){
87 < #ifdef IS_MPI
88 <    if (worldRank == 0) {
102 < #endif
103 <      int error;
104 <      error = fclose( inIdealFile );
105 <      
106 <      if( error ){
107 <        sprintf( painCave.errMsg,
108 <                 "Error closing %s\n", inIdealFileName.c_str());
109 <        simError();
110 <      }
111 <      
112 <      MemoryUtils::deletePointers(framePos);
113 <      
114 < #ifdef IS_MPI
115 <    }
116 <    strcpy( checkPointMsg, "Restraint file closed successfully." );
117 <    MPIcheckPoint();
118 < #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(){
125 <    
126 <    int i;
127 <    unsigned int j;
128 <    
129 < #ifdef IS_MPI
130 <    int done, which_node, which_atom; // loop counter
131 < #endif //is_mpi
132 <    
133 <    const int BUFFERSIZE = 2000; // size of the read buffer
134 <    int nTotObjs; // the number of atoms
135 <    char read_buffer[BUFFERSIZE]; //the line buffer for reading
136 <    
137 <    char *eof_test; // ptr to see when we reach the end of the file
138 <    char *parseErr;
139 <    
140 <    std::vector<StuntDouble*> integrableObjects;
141 <    
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;
97 <    SimInfo::MoleculeIterator mi;
98 <    Molecule::IntegrableObjectIterator ii;
99 <    
100 < #ifndef IS_MPI
101 <    
102 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
103 <    if( eof_test == NULL ){
104 <      sprintf( painCave.errMsg,
105 <               "RestraintReader error: error reading 1st line of \"%s\"\n",
106 <               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)) {
96 >    StuntDouble* sd;
97 >
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 >    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);
189 <           integrableObject != NULL;
190 <           integrableObject = mol->nextIntegrableObject(ii)) {    
191 <        
192 <        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 <        }
203 <        
204 <        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;
108 >      if (data != NULL) {
109  
110 <    nTotObjs = info_->getNGlobalIntegrableObjects();
228 <    haveError = 0;
110 >        // make sure we can reinterpret the generic data as restraint data:
111  
112 <    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 <      }
112 >        RestraintData* restData= dynamic_cast<RestraintData*>(data);        
113  
114 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
115 <        int which_node = info_->getMolToProc(i);
116 <        
117 <        if(which_node == masterNode){
118 <          //molecules belong to master node
119 <          
120 <          mol = info_->getMoleculeByGlobalIndex(i);
121 <          
122 <          if(mol == NULL) {
123 <            sprintf(painCave.errMsg,
124 <                    "RestReader Error: Molecule not found on node %d!\n",
125 <                    worldRank);
126 <            painCave.isFatal = 1;
127 <            simError();
128 <          }
129 <          
279 <          for (integrableObject = mol->beginIntegrableObject(ii);
280 <               integrableObject != NULL;
281 <               integrableObject = mol->nextIntegrableObject(ii)){
114 >        if (restData != NULL) {
115 >
116 >          // make sure we can reinterpet the restraint data as a
117 >          // pointer to a MolecularRestraint:
118 >
119 >          MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
120 >
121 >          if (mRest != NULL) {          
122 >
123 >            // now we need to pack the stunt doubles for the reference
124 >            // structure:
125 >
126 >            std::vector<Vector3d> ref;
127 >            int count = 0;
128 >            RealType mass, mTot;
129 >            Vector3d COM(0.0);
130              
131 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
131 >            mTot = 0.0;
132              
133 <            if(eof_test == NULL){
134 <              sprintf(painCave.errMsg,
287 <                      "RestReader Error: error in reading file %s\n"
288 <                      "natoms  = %d; index = %d\n"
289 <                      "error reading the line from the file.\n",
290 <                      inIdealFileName.c_str(), nTotObjs, i );
291 <              painCave.isFatal = 1;
292 <              simError();
293 <            }
294 <        
295 <            parseIdealLine(read_buffer, integrableObject);
296 <        
297 <          }
298 <        } else {
299 <          //molecule belongs to slave nodes
300 <          
301 <          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
302 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
303 <          
304 <          for(j=0; j < nCurObj; j++){
133 >            // loop over the stunt doubles in this molecule in the order we
134 >            // will be looping them in the restraint code:
135              
136 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inIdealFile);
137 <            if(eof_test == NULL){
138 <              sprintf(painCave.errMsg,
139 <                      "RestReader Error: error in reading file %s\n"
140 <                      "natoms  = %d; index = %d\n"
141 <                      "error reading the line from the file.\n",
142 <                      inIdealFileName.c_str(), nTotObjs, i );
143 <              painCave.isFatal = 1;
144 <              simError();
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              }
151              
152 <            MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
153 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
152 >            COM /= mTot;
153 >
154 >            mRest->setReferenceStructure(ref, COM);        
155 >
156            }
157          }
158        }
322 –    } else {
323 –      //actions taken at slave nodes
324 –      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
325 –        int which_node = info_->getMolToProc(i);
326 –        
327 –        if(which_node == worldRank){
328 –          //molecule with global index i belongs to this processor
329 –          
330 –          mol = info_->getMoleculeByGlobalIndex(i);
331 –          
332 –          if(mol == NULL) {
333 –            sprintf(painCave.errMsg,
334 –                    "RestReader Error: molecule not found on node %d\n",
335 –                    worldRank);
336 –            painCave.isFatal = 1;
337 –            simError();
338 –          }
339 –          
340 –          nCurObj = mol->getNIntegrableObjects();
341 –          
342 –          MPI_Send(&nCurObj, 1, MPI_INT, masterNode,
343 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
344 –          
345 –          for (integrableObject = mol->beginIntegrableObject(ii);
346 –               integrableObject != NULL;
347 –               integrableObject = mol->nextIntegrableObject(ii)){
348 –            
349 –            MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, masterNode,
350 –                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
351 –            
352 –            parseErr = parseIdealLine(read_buffer, integrableObject);
353 –            
354 –            if( parseErr != NULL ){
355 –              strcpy( painCave.errMsg, parseErr );
356 –              simError();
357 –            }
358 –          }
359 –        }
360 –      }
159      }
362 – #endif
160    }
364 –  
365 –  char* RestReader::parseIdealLine(char* readLine, StuntDouble* sd){
366 –    
367 –    char *foo; // the pointer to the current string token
368 –    
369 –    RealType pos[3];        // position place holders
370 –    RealType q[4];          // the quaternions
371 –    RealType RfromQ[3][3];  // the rotation matrix
372 –    RealType normalize;     // to normalize the reference unit vector
373 –    RealType uX, uY, uZ;    // reference unit vector place holders
374 –    RealType uselessToken;
375 –    StringTokenizer tokenizer(readLine);
376 –    int nTokens;
377 –    
378 –    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,
392 <              "RestReader Error: Atom type [%s] in %s does not "
393 <              "match Atom Type [%s] in .md file.\n",
394 <              name.c_str(), inIdealFileName.c_str(),
395 <              sd->getType().c_str());
396 <      painCave.isFatal = 1;
397 <      simError();        
398 <    }
399 <    
400 <    pos[0] = tokenizer.nextTokenAsDouble();
401 <    pos[1] = tokenizer.nextTokenAsDouble();
402 <    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");
407 <    refPosX->setData(pos[0]);
408 <    sd->addProperty(refPosX);
409 <
410 <    DoubleGenericData* refPosY = new DoubleGenericData();
411 <    refPosY->setID("refPosY");
412 <    refPosY->setData(pos[1]);
413 <    sd->addProperty(refPosY);
414 <    
415 <    DoubleGenericData* refPosZ = new DoubleGenericData();
416 <    refPosZ->setID("refPosZ");
417 <    refPosZ->setData(pos[2]);
418 <    sd->addProperty(refPosZ);
181 >    if (integrableObject == NULL) {
182 >      return;
183 >    }  
184  
185 <    // we don't need the velocities
186 <    uselessToken = tokenizer.nextTokenAsDouble();
187 <    uselessToken = tokenizer.nextTokenAsDouble();
188 <    uselessToken = tokenizer.nextTokenAsDouble();
189 <    
190 <    if (sd->isDirectional()) {
191 <      
192 <      q[0] = tokenizer.nextTokenAsDouble();
193 <      q[1] = tokenizer.nextTokenAsDouble();
194 <      q[2] = tokenizer.nextTokenAsDouble();
195 <      q[3] = tokenizer.nextTokenAsDouble();
196 <      
197 <      // now build the rotation matrix and find the unit vectors
198 <      RfromQ[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
199 <      RfromQ[0][1] = 2*(q[1]*q[2] + q[0]*q[3]);
200 <      RfromQ[0][2] = 2*(q[1]*q[3] - q[0]*q[2]);
201 <      RfromQ[1][0] = 2*(q[1]*q[2] - q[0]*q[3]);
202 <      RfromQ[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
203 <      RfromQ[1][2] = 2*(q[2]*q[3] + q[0]*q[1]);
204 <      RfromQ[2][0] = 2*(q[1]*q[3] + q[0]*q[2]);
205 <      RfromQ[2][1] = 2*(q[2]*q[3] - q[0]*q[1]);
206 <      RfromQ[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
207 <      
208 <      normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1]
209 <                       + RfromQ[2][2]*RfromQ[2][2]);
210 <      uX = RfromQ[2][0]/normalize;
211 <      uY = RfromQ[2][1]/normalize;
212 <      uZ = RfromQ[2][2]/normalize;
213 <      
214 <      // store reference unit vectors as generic data in the stuntdouble
215 <      DoubleGenericData* refVectorX = new DoubleGenericData();
216 <      refVectorX->setID("refVectorX");
217 <      refVectorX->setData(uX);
218 <      sd->addProperty(refVectorX);
219 <      
220 <      DoubleGenericData* refVectorY = new DoubleGenericData();
221 <      refVectorY->setID("refVectorY");
222 <      refVectorY->setData(uY);
223 <      sd->addProperty(refVectorY);
224 <      
225 <      DoubleGenericData* refVectorZ = new DoubleGenericData();
226 <      refVectorZ->setID("refVectorZ");
227 <      refVectorZ->setData(uZ);
228 <      sd->addProperty(refVectorZ);
229 <    }
230 <    
231 <    // we don't need the angular velocities, so let's exit the line
232 <    return NULL;
233 <  }
234 <  
235 <  void RestReader::readZangle(){
236 <    
237 <    int i;
238 <    unsigned int j;
239 <    int isPresent;
240 <    
241 <    Molecule* mol;
242 <    StuntDouble* integrableObject;
243 <    SimInfo::MoleculeIterator mi;
244 <    Molecule::IntegrableObjectIterator ii;
245 <    
246 < #ifdef IS_MPI
247 <    int done, which_node, which_atom; // loop counter
248 <    int nProc;
249 <    MPI_Status istatus;
250 < #endif //is_mpi
251 <    
252 <    const int BUFFERSIZE = 2000; // size of the read buffer
253 <    int nTotObjs; // the number of atoms
254 <    char read_buffer[BUFFERSIZE]; //the line buffer for reading
255 <    
256 <    char *eof_test; // ptr to see when we reach the end of the file
257 <    char *parseErr;
258 <    
259 <    std::vector<StuntDouble*> vecParticles;
260 <    std::vector<RealType> tempZangs;
261 <      
497 <    inAngFileName = info_->getRestFileName();
498 <    
499 <    inAngFileName += "0";
500 <    
501 <    // open the omega value file for reading
502 < #ifdef IS_MPI
503 <    if (worldRank == 0) {
504 < #endif
505 <      isPresent = 1;
506 <      inAngFile = fopen(inAngFileName.c_str(), "r");
507 <      if(!inAngFile){
508 <        sprintf(painCave.errMsg,
509 <                "Restraints Warning: %s file is not present\n"
510 <                "\tAll omega values will be initialized to zero. If the\n"
511 <                "\tsimulation is starting from the idealCrystal.in reference\n"
512 <                "\tconfiguration, this is the desired action. If this is not\n"
513 <                "\tthe case, the energy calculations will be incorrect.\n",
514 <                inAngFileName.c_str());
515 <        painCave.severity = OOPSE_WARNING;
516 <        painCave.isFatal = 0;
517 <        simError();  
518 <        isPresent = 0;
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 >        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 < oopse::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 >          break;
246 >        }
247 >        case 't' : {
248 >           torque[0] = tokenizer.nextTokenAsDouble();
249 >           torque[1] = tokenizer.nextTokenAsDouble();
250 >           torque[2] = tokenizer.nextTokenAsDouble();          
251 >
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        }
520 –      
521 – #ifdef IS_MPI
522 –      // let the other nodes know the status of the file search
523 –      MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
524 – #endif // is_mpi
525 –      
526 –      if (!isPresent) {
527 –        zeroZangle();
528 –        return;
529 –      }
530 –      
531 – #ifdef IS_MPI
263      }
264 +    
265 +    // keep the position in case we need it for a molecular restraint:
266      
267 <    // listen to node 0 to see if we should exit this function
268 <    if (worldRank != 0) {
269 <      MPI_Bcast(&isPresent, 1, MPI_INT, 0, MPI_COMM_WORLD);
270 <      if (!isPresent) {
271 <        zeroZangle();
272 <        return;
267 >    all_pos_[index] = pos;
268 >
269 >    // is this io restrained?    
270 >    GenericData* data = integrableObject->getPropertyByName("Restraint");
271 >    ObjectRestraint* oRest;
272 >
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      }
289 <    
290 <    strcpy( checkPointMsg, "zAngle file opened successfully for reading." );
291 <    MPIcheckPoint();
292 < #endif
293 <    
294 < #ifndef IS_MPI
295 <    
296 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
297 <    if( eof_test == NULL ){
298 <      sprintf( painCave.errMsg,
299 <               "RestraintReader error: error reading 1st line of \"%s\"\n",
300 <               inAngFileName.c_str() );
301 <      painCave.isFatal = 1;
302 <      simError();
289 >
290 >  }
291 >  
292 >
293 >
294 >  void RestReader::readFrameProperties(std::istream& inputStream) {
295 >    inputStream.getline(buffer, bufferSize);
296 >    std::string line(buffer);
297 >
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 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
306 <    while ( eof_test != NULL ) {
307 <      // check for and ignore blank lines
308 <      if ( read_buffer != NULL )
309 <        tempZangs.push_back( atof(read_buffer) );
310 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
311 <    }
312 <    
313 <    nTotObjs = info_->getNGlobalIntegrableObjects();
314 <    
568 <    if( nTotObjs != tempZangs.size() ){
569 <      sprintf( painCave.errMsg,
570 <               "RestraintReader zAngle reading error. %s nIntegrable, %d, "
571 <               "does not match the meta-data file's nIntegrable, %d.\n",
572 <               inAngFileName.c_str(), tempZangs.size(), nTotObjs );
573 <      painCave.isFatal = 1;
574 <      simError();
575 <    }
576 <    
577 <    // load the zAngles into the integrable objects
578 <    i = 0;
579 <    for (mol = info_->beginMolecule(mi); mol != NULL;
580 <         mol = info_->nextMolecule(mi)) {
304 >
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 >    while(inputStream.getline(buffer, bufferSize)) {
310 >      line = buffer;
311 >      
312 >      if(line.find("</FrameData>") != std::string::npos) {
313 >        break;
314 >      }
315        
582 –      for (integrableObject = mol->beginIntegrableObject(ii);
583 –           integrableObject != NULL;
584 –           integrableObject = mol->nextIntegrableObject(ii)) {    
585 –        
586 –        integrableObject->setZangle(tempZangs[i]);
587 –        i++;
588 –      }
316      }
590 –    
591 –    // MPI Section of code..........
592 – #else //IS_MPI
593 –    
594 –    // first thing first, suspend fatalities.
595 –    painCave.isEventLoop = 1;
317  
597 –    int masterNode = 0;
598 –    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
599 –    int haveError;
600 –    int index;    
601 –
602 –    int nCurObj;
603 –    RealType angleTranfer;
604 –    
605 –    nTotObjs = info_->getNGlobalIntegrableObjects();
606 –    haveError = 0;
607 –
608 –    if (worldRank == masterNode) {
609 –      
610 –      eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
611 –      if( eof_test == NULL ){
612 –        sprintf( painCave.errMsg,
613 –                 "Error reading 1st line of %s \n ",inAngFileName.c_str());
614 –        haveError = 1;
615 –        simError();
616 –      }
617 –      
618 –      // let node 0 load the temporary angle vector
619 –      eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
620 –      while ( eof_test != NULL ) {
621 –        // check for and ignore blank lines
622 –        if ( read_buffer != NULL )
623 –          tempZangs.push_back( atof(read_buffer) );
624 –        eof_test = fgets(read_buffer, sizeof(read_buffer), inAngFile);
625 –      }
626 –      
627 –      // Check to see that the number of integrable objects in the
628 –      // intial configuration file is the same as derived from the
629 –      // meta-data file.
630 –      if( nTotObjs != tempZangs.size() ){
631 –        sprintf( painCave.errMsg,
632 –                 "RestraintReader zAngle reading Error. %s nIntegrable, %d, "
633 –                 "does not match the meta-data file's nIntegrable, %d.\n",
634 –                 inAngFileName.c_str(), tempZangs.size(), nTotObjs);
635 –        haveError= 1;
636 –        simError();
637 –      }
638 –      
639 –      // At this point, node 0 has a tempZangs vector completed, and
640 –      // everyone else has nada
641 –      index = 0;
642 –      
643 –      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
644 –        // Get the Node number which has this atom
645 –        which_node = info_->getMolToProc(i);
646 –        
647 –        if (worldRank == masterNode) {
648 –          mol = info_->getMoleculeByGlobalIndex(i);
649 –          
650 –          if(mol == NULL) {
651 –            strcpy(painCave.errMsg, "Molecule not found on node 0!");
652 –            haveError = 1;
653 –            simError();
654 –          }
655 –          
656 –          for (integrableObject = mol->beginIntegrableObject(ii);
657 –               integrableObject != NULL;
658 –               integrableObject = mol->nextIntegrableObject(ii)){
659 –            
660 –            integrableObject->setZangle(tempZangs[index]);
661 –            index++;
662 –          }    
663 –          
664 –        } else {
665 –          // I am MASTER OF THE UNIVERSE, but I don't own this molecule
666 –          
667 –          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
668 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
669 –          
670 –          for(j=0; j < nCurObj; j++){            
671 –            angleTransfer = tempZangs[index];
672 –            MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node,
673 –                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);              
674 –            index++;
675 –          }
676 –          
677 –        }
678 –      }
679 –    } else {
680 –      // I am SLAVE TO THE MASTER
681 –      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
682 –        int which_node = info_->getMolToProc(i);
683 –
684 –        if (which_node == worldRank) {
685 –          
686 –          // BUT I OWN THIS MOLECULE!!!
687 –          
688 –          mol = info_->getMoleculeByGlobalIndex(i);
689 –
690 –          if(mol == NULL) {
691 –            sprintf(painCave.errMsg,
692 –                    "RestReader Error: molecule not found on node %d\n",
693 –                    worldRank);
694 –            painCave.isFatal = 1;
695 –            simError();
696 –          }
697 –
698 –          nCurObj = mol->getNIntegrableObjects();
699 –        
700 –          MPI_Send(&nCurObj, 1, MPI_INT, 0,
701 –                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
702 –          
703 –          for (integrableObject = mol->beginIntegrableObject(ii);
704 –               integrableObject != NULL;
705 –               integrableObject = mol->nextIntegrableObject(ii)){
706 –            
707 –            MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0,
708 –                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
709 –
710 –            integrableObject->setZangle(angleTransfer);
711 –          }    
712 –        }
713 –      }
714 –    }
715 – #endif
318    }
717 –  
718 –  void RestReader :: zeroZangle(){
719 –    
720 –    int i;
721 –    unsigned int j;
722 –    int nTotObjs; // the number of atoms
723 –    
724 –    Molecule* mol;
725 –    StuntDouble* integrableObject;
726 –    SimInfo::MoleculeIterator mi;
727 –    Molecule::IntegrableObjectIterator ii;
728 –    
729 –    std::vector<StuntDouble*> vecParticles;
730 –    
731 – #ifndef IS_MPI
732 –    // set all zAngles to 0.0
733 –    for (mol = info_->beginMolecule(mi); mol != NULL;
734 –         mol = info_->nextMolecule(mi))
735 –      
736 –      for (integrableObject = mol->beginIntegrableObject(ii);
737 –           integrableObject != NULL;
738 –           integrableObject = mol->nextIntegrableObject(ii))    
739 –        integrableObject->setZangle( 0.0 );
740 –    
741 –    
742 –    // MPI Section of code..........
743 – #else //IS_MPI
744 –    
745 –    // first thing first, suspend fatalities.
746 –    painCave.isEventLoop = 1;
747 –    
748 –    int masterNode = 0;
749 –    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
750 –    int haveError;
751 –    int which_node;
752 –    
753 –    MPI_Status istatus;
754 –    
755 –    int nCurObj;
756 –    RealType angleTranfer;
757 –    
758 –    nTotObjs = info_->getNGlobalIntegrableObjects();
759 –    haveError = 0;
760 –    if (worldRank == masterNode) {
319  
320 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
321 <        // Get the Node number which has this atom
764 <        which_node = info_->getMolToProc(i);
765 <        
766 <        // let's let node 0 pass out constant values to all the processors
767 <        if (worldRank == masterNode) {
768 <          mol = info_->getMoleculeByGlobalIndex(i);
769 <          
770 <          if(mol == NULL) {
771 <            strcpy(painCave.errMsg, "Molecule not found on node 0!");
772 <            haveError = 1;
773 <            simError();
774 <          }
775 <          
776 <          for (integrableObject = mol->beginIntegrableObject(ii);
777 <               integrableObject != NULL;
778 <               integrableObject = mol->nextIntegrableObject(ii)){
779 <            
780 <            integrableObject->setZangle( 0.0 );
781 <            
782 <          }
783 <          
784 <        } else {
785 <          // I am MASTER OF THE UNIVERSE, but I don't own this molecule
786 <          
787 <          MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
788 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
789 <          
790 <          for(j=0; j < nCurObj; j++){            
791 <            angleTransfer = 0.0;
792 <            MPI_Send(&angleTransfer, 1, MPI_REALTYPE, which_node,
793 <                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD);              
794 <            
795 <          }
796 <        }
797 <      }
798 <    } else {
799 <      // I am SLAVE TO THE MASTER
800 <      for (i=0 ; i < info_->getNGlobalMolecules(); i++) {
801 <        int which_node = info_->getMolToProc(i);
802 <        
803 <        if (which_node == worldRank) {
804 <          
805 <          // BUT I OWN THIS MOLECULE!!!
806 <          mol = info_->getMoleculeByGlobalIndex(i);
807 <          
808 <          if(mol == NULL) {
809 <            sprintf(painCave.errMsg,
810 <                    "RestReader Error: molecule not found on node %d\n",
811 <                    worldRank);
812 <            painCave.isFatal = 1;
813 <            simError();
814 <          }
815 <          
816 <          nCurObj = mol->getNIntegrableObjects();
817 <          
818 <          MPI_Send(&nCurObj, 1, MPI_INT, 0,
819 <                   TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
820 <          
821 <          for (integrableObject = mol->beginIntegrableObject(ii);
822 <               integrableObject != NULL;
823 <               integrableObject = mol->nextIntegrableObject(ii)){
824 <            
825 <            MPI_Recv(&angleTransfer, 1, MPI_REALTYPE, 0,
826 <                     TAKE_THIS_TAG_DOUBLE, MPI_COMM_WORLD, &istatus);
827 <            vecParticles[j]->setZangle(angleTransfer);
828 <          }    
829 <        }
830 <      }
831 <    }
832 < #endif
833 <  }
834 <  
835 < } // end namespace oopse
320 >  
321 > }//end namespace oopse

Diff Legend

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