1 |
|
/* |
2 |
< |
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
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 |
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 "primitives/Molecule.hpp" |
41 |
> |
|
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 "restraints/ObjectRestraint.hpp" |
60 |
< |
#include "restraints/MolecularRestraint.hpp" |
60 |
> |
#include "restraints/MolecularRestraint.hpp" |
61 |
> |
|
62 |
> |
#ifdef IS_MPI |
63 |
> |
|
64 |
> |
#include <mpi.h> |
65 |
> |
#endif |
66 |
|
|
67 |
|
namespace OpenMD { |
68 |
+ |
|
69 |
+ |
void RestReader::scanFile(){ |
70 |
+ |
int lineNo = 0; |
71 |
+ |
std::streampos prevPos; |
72 |
+ |
std::streampos currPos; |
73 |
+ |
|
74 |
+ |
#ifdef IS_MPI |
75 |
+ |
|
76 |
+ |
if (worldRank == 0) { |
77 |
+ |
#endif // is_mpi |
78 |
+ |
|
79 |
+ |
inFile_->clear(); |
80 |
+ |
currPos = inFile_->tellg(); |
81 |
+ |
prevPos = currPos; |
82 |
|
|
83 |
< |
void RestReader::readReferenceStructure() { |
83 |
> |
bool foundOpenSnapshotTag = false; |
84 |
> |
|
85 |
> |
while(!foundOpenSnapshotTag && inFile_->getline(buffer, bufferSize)) { |
86 |
> |
++lineNo; |
87 |
> |
|
88 |
> |
std::string line = buffer; |
89 |
> |
currPos = inFile_->tellg(); |
90 |
> |
if (line.find("<Snapshot>")!= std::string::npos) { |
91 |
> |
foundOpenSnapshotTag = true; |
92 |
> |
framePos_ = prevPos; |
93 |
> |
} |
94 |
> |
prevPos = currPos; |
95 |
> |
} |
96 |
> |
|
97 |
> |
#ifdef IS_MPI |
98 |
> |
} |
99 |
> |
MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD); |
100 |
> |
#endif // is_mpi |
101 |
> |
} |
102 |
|
|
52 |
– |
// some of this comes directly from DumpReader. |
103 |
|
|
104 |
< |
if (!isScanned_) |
105 |
< |
scanFile(); |
106 |
< |
|
107 |
< |
int storageLayout = info_->getSnapshotManager()->getStorageLayout(); |
108 |
< |
|
109 |
< |
if (storageLayout & DataStorage::dslPosition) { |
110 |
< |
needPos_ = true; |
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 |
> |
inFile_->clear(); |
121 |
> |
inFile_->seekg(framePos_); |
122 |
> |
|
123 |
> |
while (inFile_->getline(buffer, bufferSize)) { |
124 |
> |
|
125 |
> |
line = buffer; |
126 |
> |
sendBuffer += line; |
127 |
> |
sendBuffer += '\n'; |
128 |
> |
if (line.find("</Snapshot>") != std::string::npos) { |
129 |
> |
break; |
130 |
> |
} |
131 |
> |
} |
132 |
> |
|
133 |
> |
int sendBufferSize = sendBuffer.size(); |
134 |
> |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
135 |
> |
MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
136 |
> |
|
137 |
> |
sstream.str(sendBuffer); |
138 |
|
} else { |
139 |
< |
needPos_ = false; |
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 |
> |
inputStream.getline(buffer, bufferSize); |
153 |
> |
|
154 |
> |
line = buffer; |
155 |
> |
if (line.find("<Snapshot>") == std::string::npos) { |
156 |
> |
sprintf(painCave.errMsg, |
157 |
> |
"RestReader Error: can not find <Snapshot>\n"); |
158 |
> |
painCave.isFatal = 1; |
159 |
> |
simError(); |
160 |
|
} |
161 |
< |
|
162 |
< |
needVel_ = false; |
163 |
< |
|
164 |
< |
if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) { |
165 |
< |
needQuaternion_ = true; |
166 |
< |
} else { |
167 |
< |
needQuaternion_ = false; |
168 |
< |
} |
161 |
> |
|
162 |
> |
//read frameData |
163 |
> |
readFrameProperties(inputStream); |
164 |
> |
|
165 |
> |
//read StuntDoubles |
166 |
> |
readStuntDoubles(inputStream); |
167 |
> |
|
168 |
> |
inputStream.getline(buffer, bufferSize); |
169 |
> |
line = buffer; |
170 |
> |
if (line.find("</Snapshot>") == std::string::npos) { |
171 |
> |
sprintf(painCave.errMsg, |
172 |
> |
"RestReader Error: can not find </Snapshot>\n"); |
173 |
> |
painCave.isFatal = 1; |
174 |
> |
simError(); |
175 |
> |
} |
176 |
> |
} |
177 |
> |
|
178 |
> |
void RestReader::readReferenceStructure() { |
179 |
|
|
73 |
– |
needAngMom_ = false; |
74 |
– |
|
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 |
< |
|
82 |
< |
|
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 |
< |
readSet(0); |
193 |
< |
|
192 |
> |
scanFile(); |
193 |
> |
|
194 |
> |
|
195 |
> |
readSet(); |
196 |
> |
|
197 |
> |
|
198 |
|
// all ObjectRestraints have been handled, now we have to worry about |
199 |
|
// molecular restraints: |
200 |
|
|
202 |
|
Molecule::IntegrableObjectIterator j; |
203 |
|
Molecule* mol; |
204 |
|
StuntDouble* sd; |
205 |
< |
|
205 |
> |
|
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 |
< |
|
209 |
> |
|
210 |
|
for (mol = info_->beginMolecule(i); mol != NULL; |
211 |
|
mol = info_->nextMolecule(i)) { |
212 |
< |
|
212 |
> |
|
213 |
|
// is this molecule restrained? |
214 |
|
GenericData* data = mol->getPropertyByName("Restraint"); |
215 |
|
|
216 |
|
if (data != NULL) { |
217 |
< |
|
217 |
> |
|
218 |
|
// make sure we can reinterpret the generic data as restraint data: |
219 |
< |
|
219 |
> |
|
220 |
|
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
221 |
< |
|
221 |
> |
|
222 |
|
if (restData != NULL) { |
223 |
< |
|
223 |
> |
|
224 |
|
// make sure we can reinterpet the restraint data as a |
225 |
|
// pointer to a MolecularRestraint: |
226 |
< |
|
226 |
> |
|
227 |
|
MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData()); |
228 |
< |
|
228 |
> |
|
229 |
|
if (mRest != NULL) { |
230 |
< |
|
230 |
> |
|
231 |
|
// now we need to pack the stunt doubles for the reference |
232 |
|
// structure: |
233 |
< |
|
233 |
> |
|
234 |
|
std::vector<Vector3d> ref; |
235 |
|
int count = 0; |
236 |
|
RealType mass, mTot; |
248 |
|
// doubles from the *globally* sorted array of |
249 |
|
// positions: |
250 |
|
|
251 |
< |
ref.push_back( all_pos_[sd->getGlobalIndex()] ); |
252 |
< |
|
145 |
< |
mass = sd->getMass(); |
146 |
< |
|
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 |
|
} |
151 |
– |
|
257 |
|
COM /= mTot; |
153 |
– |
|
258 |
|
mRest->setReferenceStructure(ref, COM); |
155 |
– |
|
259 |
|
} |
260 |
|
} |
261 |
|
} |
265 |
|
|
266 |
|
|
267 |
|
void RestReader::parseDumpLine(const std::string& line) { |
268 |
+ |
|
269 |
|
StringTokenizer tokenizer(line); |
270 |
|
int nTokens; |
271 |
|
|
273 |
|
|
274 |
|
if (nTokens < 2) { |
275 |
|
sprintf(painCave.errMsg, |
276 |
< |
"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); |
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(); |
178 |
– |
|
179 |
– |
StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index); |
282 |
|
|
283 |
< |
if (integrableObject == NULL) { |
284 |
< |
return; |
285 |
< |
} |
286 |
< |
|
287 |
< |
std::string type = tokenizer.nextToken(); |
288 |
< |
int size = type.size(); |
289 |
< |
Vector3d pos; |
290 |
< |
Vector3d vel; |
291 |
< |
Quat4d q; |
292 |
< |
Vector3d ji; |
293 |
< |
Vector3d force; |
294 |
< |
Vector3d torque; |
295 |
< |
|
296 |
< |
for(int i = 0; i < size; ++i) { |
297 |
< |
switch(type[i]) { |
283 |
> |
|
284 |
> |
if (index < 1116){ |
285 |
> |
std::string type = tokenizer.nextToken(); |
286 |
> |
|
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 |
> |
bool exist = false; |
297 |
> |
int indexCount = 0; |
298 |
> |
|
299 |
> |
while ( (!exist) && (indexCount < stuntDoubleIndex_.size())){ |
300 |
> |
if (index == stuntDoubleIndex_[indexCount]) |
301 |
> |
exist = true; |
302 |
> |
indexCount++; |
303 |
> |
} |
304 |
> |
|
305 |
> |
StuntDouble* integrableObject; |
306 |
> |
|
307 |
> |
if (exist){ |
308 |
|
|
309 |
< |
case 'p': { |
310 |
< |
pos[0] = tokenizer.nextTokenAsDouble(); |
311 |
< |
pos[1] = tokenizer.nextTokenAsDouble(); |
312 |
< |
pos[2] = tokenizer.nextTokenAsDouble(); |
313 |
< |
break; |
314 |
< |
} |
315 |
< |
case 'v' : { |
316 |
< |
vel[0] = tokenizer.nextTokenAsDouble(); |
317 |
< |
vel[1] = tokenizer.nextTokenAsDouble(); |
318 |
< |
vel[2] = tokenizer.nextTokenAsDouble(); |
319 |
< |
break; |
320 |
< |
} |
321 |
< |
|
322 |
< |
case 'q' : { |
323 |
< |
if (integrableObject->isDirectional()) { |
324 |
< |
|
325 |
< |
q[0] = tokenizer.nextTokenAsDouble(); |
326 |
< |
q[1] = tokenizer.nextTokenAsDouble(); |
327 |
< |
q[2] = tokenizer.nextTokenAsDouble(); |
328 |
< |
q[3] = tokenizer.nextTokenAsDouble(); |
329 |
< |
|
330 |
< |
RealType qlen = q.length(); |
331 |
< |
if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0 |
332 |
< |
|
333 |
< |
sprintf(painCave.errMsg, |
334 |
< |
"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n"); |
335 |
< |
painCave.isFatal = 1; |
336 |
< |
simError(); |
337 |
< |
|
226 |
< |
} |
227 |
< |
|
228 |
< |
q.normalize(); |
229 |
< |
} |
230 |
< |
break; |
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 |
< |
case 'j' : { |
340 |
< |
if (integrableObject->isDirectional()) { |
341 |
< |
ji[0] = tokenizer.nextTokenAsDouble(); |
342 |
< |
ji[1] = tokenizer.nextTokenAsDouble(); |
343 |
< |
ji[2] = tokenizer.nextTokenAsDouble(); |
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 |
|
} |
238 |
– |
break; |
239 |
– |
} |
240 |
– |
case 'f': { |
241 |
– |
force[0] = tokenizer.nextTokenAsDouble(); |
242 |
– |
force[1] = tokenizer.nextTokenAsDouble(); |
243 |
– |
force[2] = tokenizer.nextTokenAsDouble(); |
244 |
– |
|
245 |
– |
break; |
362 |
|
} |
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 |
– |
|
363 |
|
} |
364 |
|
} |
365 |
< |
|
366 |
< |
// keep the position in case we need it for a molecular restraint: |
365 |
> |
} |
366 |
> |
|
367 |
> |
void RestReader::readStuntDoubles(std::istream& inputStream) { |
368 |
|
|
369 |
< |
all_pos_[index] = pos; |
370 |
< |
|
371 |
< |
// is this io restrained? |
372 |
< |
GenericData* data = integrableObject->getPropertyByName("Restraint"); |
373 |
< |
ObjectRestraint* oRest; |
374 |
< |
|
375 |
< |
if (data != NULL) { |
376 |
< |
// make sure we can reinterpret the generic data as restraint data: |
377 |
< |
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
378 |
< |
if (restData != NULL) { |
379 |
< |
// make sure we can reinterpet the restraint data as a pointer to |
380 |
< |
// an ObjectRestraint: |
381 |
< |
oRest = dynamic_cast<ObjectRestraint*>(restData->getData()); |
382 |
< |
if (oRest != NULL) { |
383 |
< |
if (integrableObject->isDirectional()) { |
282 |
< |
oRest->setReferenceStructure(pos, q.toRotationMatrix3()); |
283 |
< |
} else { |
284 |
< |
oRest->setReferenceStructure(pos); |
285 |
< |
} |
286 |
< |
} |
369 |
> |
inputStream.getline(buffer, bufferSize); |
370 |
> |
std::string line(buffer); |
371 |
> |
|
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 |
> |
while(inputStream.getline(buffer, bufferSize)) { |
380 |
> |
line = buffer; |
381 |
> |
|
382 |
> |
if(line.find("</StuntDoubles>") != std::string::npos) { |
383 |
> |
break; |
384 |
|
} |
385 |
+ |
|
386 |
+ |
parseDumpLine(line); |
387 |
|
} |
388 |
< |
|
388 |
> |
|
389 |
|
} |
291 |
– |
|
390 |
|
|
391 |
< |
|
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 |
< |
"DumpReader Error: Missing <FrameData>\n"); |
398 |
> |
"RestReader Error: Missing <FrameData>\n"); |
399 |
|
painCave.isFatal = 1; |
400 |
|
simError(); |
401 |
|
} |
412 |
|
} |
413 |
|
|
414 |
|
} |
415 |
< |
|
415 |
> |
|
416 |
|
} |
417 |
|
|
418 |
|
|