1 |
/* |
2 |
* Copyright (c) 2009 The University of Notre Dame. All Rights Reserved. |
3 |
* |
4 |
* The University of Notre Dame grants you ("Licensee") a |
5 |
* non-exclusive, royalty free, license to use, modify and |
6 |
* redistribute this software in source and binary code form, provided |
7 |
* that the following conditions are met: |
8 |
* |
9 |
* 1. Redistributions of source code must retain the above copyright |
10 |
* notice, this list of conditions and the following disclaimer. |
11 |
* |
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. |
16 |
* |
17 |
* This software is provided "AS IS," without a warranty of any |
18 |
* kind. All express or implied conditions, representations and |
19 |
* warranties, including any implied warranty of merchantability, |
20 |
* fitness for a particular purpose or non-infringement, are hereby |
21 |
* excluded. The University of Notre Dame and its licensors shall not |
22 |
* be liable for any damages suffered by licensee as a result of |
23 |
* using, modifying or distributing the software or its |
24 |
* derivatives. In no event will the University of Notre Dame or its |
25 |
* licensors be liable for any lost revenue, profit or data, or for |
26 |
* direct, indirect, special, consequential, incidental or punitive |
27 |
* damages, however caused and regardless of the theory of liability, |
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 |
* 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 |
|
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" |
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 |
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 |
|
103 |
|
104 |
void RestReader::readSet(){ |
105 |
std::string line; |
106 |
|
107 |
#ifndef IS_MPI |
108 |
|
109 |
inFile_->clear(); |
110 |
inFile_->seekg(framePos_); |
111 |
|
112 |
std::istream& inputStream = *inFile_; |
113 |
#else |
114 |
|
115 |
int masterNode = 0; |
116 |
std::stringstream sstream; |
117 |
if (worldRank == masterNode) { |
118 |
std::string sendBuffer; |
119 |
|
120 |
inFile_->clear(); |
121 |
inFile_->seekg(framePos_); |
122 |
|
123 |
while (inFile_->getline(buffer, bufferSize)) { |
124 |
|
125 |
line = buffer; |
126 |
sendBuffer += line; |
127 |
sendBuffer += '\n'; |
128 |
if (line.find("</Snapshot>") != std::string::npos) { |
129 |
break; |
130 |
} |
131 |
} |
132 |
|
133 |
int sendBufferSize = sendBuffer.size(); |
134 |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
135 |
MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
136 |
|
137 |
sstream.str(sendBuffer); |
138 |
} else { |
139 |
int sendBufferSize; |
140 |
MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
141 |
char * recvBuffer = new char[sendBufferSize+1]; |
142 |
assert(recvBuffer); |
143 |
recvBuffer[sendBufferSize] = '\0'; |
144 |
MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
145 |
sstream.str(recvBuffer); |
146 |
delete [] recvBuffer; |
147 |
} |
148 |
|
149 |
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 |
//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 |
|
180 |
// We need temporary storage to keep track of all StuntDouble positions |
181 |
// in case some of the restraints are molecular (i.e. if they use |
182 |
// multiple SD positions to determine restrained orientations or positions: |
183 |
|
184 |
all_pos_.clear(); |
185 |
all_pos_.resize(info_->getNGlobalIntegrableObjects()) ; |
186 |
|
187 |
// Restraint files are just standard dump files, but with the reference |
188 |
// structure stored in the first frame (frame 0). |
189 |
// RestReader overloads readSet and explicitly handles all of the |
190 |
// ObjectRestraints in that method: |
191 |
|
192 |
scanFile(); |
193 |
|
194 |
|
195 |
readSet(); |
196 |
|
197 |
|
198 |
// all ObjectRestraints have been handled, now we have to worry about |
199 |
// molecular restraints: |
200 |
|
201 |
SimInfo::MoleculeIterator i; |
202 |
Molecule::IntegrableObjectIterator j; |
203 |
Molecule* mol; |
204 |
StuntDouble* sd; |
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 |
|
210 |
for (mol = info_->beginMolecule(i); mol != NULL; |
211 |
mol = info_->nextMolecule(i)) { |
212 |
|
213 |
// is this molecule restrained? |
214 |
GenericData* data = mol->getPropertyByName("Restraint"); |
215 |
|
216 |
if (data != NULL) { |
217 |
|
218 |
// make sure we can reinterpret the generic data as restraint data: |
219 |
|
220 |
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
221 |
|
222 |
if (restData != NULL) { |
223 |
|
224 |
// make sure we can reinterpet the restraint data as a |
225 |
// pointer to a MolecularRestraint: |
226 |
|
227 |
MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData()); |
228 |
|
229 |
if (mRest != NULL) { |
230 |
|
231 |
// now we need to pack the stunt doubles for the reference |
232 |
// structure: |
233 |
|
234 |
std::vector<Vector3d> ref; |
235 |
int count = 0; |
236 |
RealType mass, mTot; |
237 |
Vector3d COM(0.0); |
238 |
|
239 |
mTot = 0.0; |
240 |
|
241 |
// loop over the stunt doubles in this molecule in the order we |
242 |
// will be looping them in the restraint code: |
243 |
|
244 |
for (sd = mol->beginIntegrableObject(j); sd != NULL; |
245 |
sd = mol->nextIntegrableObject(j)) { |
246 |
|
247 |
// push back the reference positions of the stunt |
248 |
// doubles from the *globally* sorted array of |
249 |
// positions: |
250 |
|
251 |
ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] ); |
252 |
mass = sd->getMass(); |
253 |
COM = COM + mass * ref[count]; |
254 |
mTot = mTot + mass; |
255 |
count = count + 1; |
256 |
} |
257 |
COM /= mTot; |
258 |
mRest->setReferenceStructure(ref, COM); |
259 |
} |
260 |
} |
261 |
} |
262 |
} |
263 |
} |
264 |
|
265 |
|
266 |
|
267 |
void RestReader::parseDumpLine(const std::string& line) { |
268 |
|
269 |
StringTokenizer tokenizer(line); |
270 |
int nTokens; |
271 |
|
272 |
nTokens = tokenizer.countTokens(); |
273 |
|
274 |
if (nTokens < 2) { |
275 |
sprintf(painCave.errMsg, |
276 |
"RestReader Error: Not enough Tokens.\n%s\n", line.c_str()); |
277 |
painCave.isFatal = 1; |
278 |
simError(); |
279 |
} |
280 |
|
281 |
int index = tokenizer.nextTokenAsInt(); |
282 |
|
283 |
|
284 |
if (index < 1116){ |
285 |
std::string type = tokenizer.nextToken(); |
286 |
|
287 |
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 |
integrableObject = info_->getIOIndexToIntegrableObject(index); |
310 |
|
311 |
int compareindex = integrableObject->getGlobalIntegrableObjectIndex(); |
312 |
|
313 |
if (integrableObject == NULL) { |
314 |
return; |
315 |
} |
316 |
|
317 |
std::string type = tokenizer.nextToken(); |
318 |
|
319 |
Vector3d pos; |
320 |
|
321 |
pos[0] = tokenizer.nextTokenAsDouble(); |
322 |
pos[1] = tokenizer.nextTokenAsDouble(); |
323 |
pos[2] = tokenizer.nextTokenAsDouble(); |
324 |
|
325 |
Vector3d vel; |
326 |
vel[0] = tokenizer.nextTokenAsDouble(); |
327 |
vel[1] = tokenizer.nextTokenAsDouble(); |
328 |
vel[2] = tokenizer.nextTokenAsDouble(); |
329 |
|
330 |
|
331 |
Quat4d q; |
332 |
if (integrableObject->isDirectional()) { |
333 |
|
334 |
q[0] = tokenizer.nextTokenAsDouble(); |
335 |
q[1] = tokenizer.nextTokenAsDouble(); |
336 |
q[2] = tokenizer.nextTokenAsDouble(); |
337 |
q[3] = tokenizer.nextTokenAsDouble(); |
338 |
} |
339 |
// keep the position in case we need it for a molecular restraint: |
340 |
|
341 |
all_pos_[index] = pos; |
342 |
|
343 |
// is this io restrained? |
344 |
GenericData* data = integrableObject->getPropertyByName("Restraint"); |
345 |
ObjectRestraint* oRest; |
346 |
|
347 |
if (data != NULL) { |
348 |
// make sure we can reinterpret the generic data as restraint data: |
349 |
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
350 |
if (restData != NULL) { |
351 |
// make sure we can reinterpet the restraint data as a pointer to |
352 |
// an ObjectRestraint: |
353 |
oRest = dynamic_cast<ObjectRestraint*>(restData->getData()); |
354 |
if (oRest != NULL) { |
355 |
if (integrableObject->isDirectional()) { |
356 |
oRest->setReferenceStructure(pos, q.toRotationMatrix3()); |
357 |
} else { |
358 |
oRest->setReferenceStructure(pos); |
359 |
} |
360 |
} |
361 |
} |
362 |
} |
363 |
} |
364 |
} |
365 |
} |
366 |
|
367 |
void RestReader::readStuntDoubles(std::istream& inputStream) { |
368 |
|
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 |
|
389 |
} |
390 |
|
391 |
|
392 |
void RestReader::readFrameProperties(std::istream& inputStream) { |
393 |
inputStream.getline(buffer, bufferSize); |
394 |
std::string line(buffer); |
395 |
|
396 |
if (line.find("<FrameData>") == std::string::npos) { |
397 |
sprintf(painCave.errMsg, |
398 |
"RestReader Error: Missing <FrameData>\n"); |
399 |
painCave.isFatal = 1; |
400 |
simError(); |
401 |
} |
402 |
|
403 |
// restraints don't care about frame data (unless we need to wrap |
404 |
// coordinates, but we'll worry about that later), so |
405 |
// we'll just scan ahead until the end of the frame data: |
406 |
|
407 |
while(inputStream.getline(buffer, bufferSize)) { |
408 |
line = buffer; |
409 |
|
410 |
if(line.find("</FrameData>") != std::string::npos) { |
411 |
break; |
412 |
} |
413 |
|
414 |
} |
415 |
|
416 |
} |
417 |
|
418 |
|
419 |
}//end namespace OpenMD |