OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
RestReader.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "io/RestReader.hpp"
46
47#include <cmath>
48#include <cstdio>
49#include <cstdlib>
50#include <cstring>
51#include <iostream>
52#include <memory>
53#include <sstream>
54#include <string>
55
56#include <sys/stat.h>
57#include <sys/types.h>
58
59#ifdef IS_MPI
60#include <mpi.h>
61#endif
62
64#include "restraints/MolecularRestraint.hpp"
65#include "restraints/ObjectRestraint.hpp"
67#include "utils/simError.h"
68
69namespace OpenMD {
70
71 void RestReader::scanFile() {
72 std::streampos prevPos;
73 std::streampos currPos;
74
75#ifdef IS_MPI
76
77 if (worldRank == 0) {
78#endif // is_mpi
79
80 inFile_->clear();
81 currPos = inFile_->tellg();
82 prevPos = currPos;
83
84 bool foundOpenSnapshotTag = false;
85 int lineNo = 0;
86 while (!foundOpenSnapshotTag && inFile_->getline(buffer, bufferSize)) {
87 ++lineNo;
88
89 std::string line = buffer;
90 currPos = inFile_->tellg();
91 if (line.find("<Snapshot>") != std::string::npos) {
92 foundOpenSnapshotTag = true;
93 framePos_ = (long long)prevPos;
94 }
95 prevPos = currPos;
96 }
97
98#ifdef IS_MPI
99 }
100 MPI_Bcast(&framePos_, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
101#endif // is_mpi
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 primaryNode = 0;
116 std::stringstream sstream;
117 if (worldRank == primaryNode) {
118 std::string sendBuffer;
119
120 inFile_->clear();
121 inFile_->seekg(framePos_);
122
123 while (inFile_->getline(buffer, bufferSize)) {
124 line = buffer;
125 sendBuffer += line;
126 sendBuffer += '\n';
127 if (line.find("</Snapshot>") != std::string::npos) { break; }
128 }
129
130 int sendBufferSize = sendBuffer.size();
131 MPI_Bcast(&sendBufferSize, 1, MPI_INT, primaryNode, MPI_COMM_WORLD);
132 MPI_Bcast((void*)sendBuffer.c_str(), sendBufferSize, MPI_CHAR,
133 primaryNode, MPI_COMM_WORLD);
134
135 sstream.str(sendBuffer);
136 } else {
137 int sendBufferSize;
138 MPI_Bcast(&sendBufferSize, 1, MPI_INT, primaryNode, MPI_COMM_WORLD);
139 char* recvBuffer = new char[sendBufferSize + 1];
140 assert(recvBuffer);
141 recvBuffer[sendBufferSize] = '\0';
142 MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, primaryNode,
143 MPI_COMM_WORLD);
144 sstream.str(recvBuffer);
145 delete[] recvBuffer;
146 }
147
148 std::istream& inputStream = sstream;
149#endif
150
151 inputStream.getline(buffer, bufferSize);
152
153 line = buffer;
154 if (line.find("<Snapshot>") == std::string::npos) {
155 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
156 "RestReader Error: can not find <Snapshot>\n");
157 painCave.isFatal = 1;
158 simError();
159 }
160
161 // read frameData
162 readFrameProperties(inputStream);
163
164 // read StuntDoubles
165 readStuntDoubles(inputStream);
166
167 inputStream.getline(buffer, bufferSize);
168 line = buffer;
169 if (line.find("</Snapshot>") == std::string::npos) {
170 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
171 "RestReader Error: can not find </Snapshot>\n");
172 painCave.isFatal = 1;
173 simError();
174 }
175 }
176
177 void RestReader::readReferenceStructure() {
178 // We need temporary storage to keep track of all StuntDouble positions
179 // in case some of the restraints are molecular (i.e. if they use
180 // multiple SD positions to determine restrained orientations or positions:
181
182 all_pos_.clear();
183 all_pos_.resize(info_->getNGlobalIntegrableObjects());
184
185 // Restraint files are just standard dump files, but with the reference
186 // structure stored in the first frame (frame 0).
187 // RestReader overloads readSet and explicitly handles all of the
188 // ObjectRestraints in that method:
189
190 scanFile();
191
192 readSet();
193
194 // all ObjectRestraints have been handled, now we have to worry about
195 // molecular restraints:
196
197 SimInfo::MoleculeIterator i;
198 Molecule::IntegrableObjectIterator j;
199 Molecule* mol;
200 StuntDouble* sd;
201
202 // no need to worry about parallel molecules, as molecules are not
203 // split across processor boundaries. Just loop over all molecules
204 // we know about:
205
206 for (mol = info_->beginMolecule(i); mol != NULL;
207 mol = info_->nextMolecule(i)) {
208 // is this molecule restrained?
209 std::shared_ptr<GenericData> data = mol->getPropertyByName("Restraint");
210
211 if (data != nullptr) {
212 // make sure we can reinterpret the generic data as restraint data:
213
214 std::shared_ptr<RestraintData> restData =
215 std::dynamic_pointer_cast<RestraintData>(data);
216
217 if (restData != nullptr) {
218 // make sure we can reinterpet the restraint data as a
219 // pointer to a MolecularRestraint:
220
221 MolecularRestraint* mRest =
222 dynamic_cast<MolecularRestraint*>(restData->getData());
223
224 if (mRest == NULL) {
225 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
226 "Can not cast RestraintData to MolecularRestraint\n");
227 painCave.severity = OPENMD_ERROR;
228 painCave.isFatal = 1;
229 simError();
230 } else {
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 // loop over the stunt doubles in this molecule in the order we
241 // will be looping them in the restraint code:
242
243 for (sd = mol->beginIntegrableObject(j); sd != NULL;
244 sd = mol->nextIntegrableObject(j)) {
245 // push back the reference positions of the stunt
246 // doubles from the *globally* sorted array of
247 // positions:
248
249 ref.push_back(all_pos_[sd->getGlobalIntegrableObjectIndex()]);
250 mass = sd->getMass();
251 COM = COM + mass * ref[count];
252 mTot = mTot + mass;
253 count = count + 1;
254 }
255 COM /= mTot;
256 mRest->setReferenceStructure(ref, COM);
257 }
258 }
259 }
260 }
261 }
262
263 void RestReader::parseDumpLine(const std::string& line) {
264 StringTokenizer tokenizer(line);
265 int nTokens;
266
267 nTokens = tokenizer.countTokens();
268
269 if (nTokens < 2) {
270 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
271 "RestReader Error: Not enough Tokens.\n%s\n", line.c_str());
272 painCave.isFatal = 1;
273 simError();
274 }
275
276 int index = tokenizer.nextTokenAsInt();
277
278 StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
279
280 if (sd == NULL) { return; }
281
282 std::string type = tokenizer.nextToken();
283 int size = type.size();
284
285 Vector3d pos;
286 Quat4d q;
287
288 for (int i = 0; i < size; ++i) {
289 switch (type[i]) {
290 case 'p': {
291 pos[0] = tokenizer.nextTokenAsDouble();
292 pos[1] = tokenizer.nextTokenAsDouble();
293 pos[2] = tokenizer.nextTokenAsDouble();
294 break;
295 }
296 case 'v': {
297 Vector3d vel;
298 vel[0] = tokenizer.nextTokenAsDouble();
299 vel[1] = tokenizer.nextTokenAsDouble();
300 vel[2] = tokenizer.nextTokenAsDouble();
301 break;
302 }
303
304 case 'q': {
305 if (sd->isDirectional()) {
306 q[0] = tokenizer.nextTokenAsDouble();
307 q[1] = tokenizer.nextTokenAsDouble();
308 q[2] = tokenizer.nextTokenAsDouble();
309 q[3] = tokenizer.nextTokenAsDouble();
310
311 RealType qlen = q.length();
312 if (qlen < OpenMD::epsilon) { // check quaternion is not equal to 0
313
314 snprintf(
315 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
316 "RestReader Error: initial quaternion error (q0^2 + q1^2 + "
317 "q2^2 + "
318 "q3^2) ~ 0\n");
319 painCave.isFatal = 1;
320 simError();
321 }
322
323 q.normalize();
324 }
325 break;
326 }
327 case 'j': {
328 Vector3d ji;
329 if (sd->isDirectional()) {
330 ji[0] = tokenizer.nextTokenAsDouble();
331 ji[1] = tokenizer.nextTokenAsDouble();
332 ji[2] = tokenizer.nextTokenAsDouble();
333 }
334 break;
335 }
336 case 'f': {
337 Vector3d force;
338 force[0] = tokenizer.nextTokenAsDouble();
339 force[1] = tokenizer.nextTokenAsDouble();
340 force[2] = tokenizer.nextTokenAsDouble();
341 break;
342 }
343 case 't': {
344 Vector3d torque;
345 torque[0] = tokenizer.nextTokenAsDouble();
346 torque[1] = tokenizer.nextTokenAsDouble();
347 torque[2] = tokenizer.nextTokenAsDouble();
348 break;
349 }
350 default: {
351 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
352 "RestReader Error: %s is an unrecognized type\n",
353 type.c_str());
354 painCave.isFatal = 1;
355 simError();
356 break;
357 }
358 }
359 // keep the position in case we need it for a molecular restraint:
360
361 all_pos_[index] = pos;
362
363 // is this io restrained?
364 std::shared_ptr<GenericData> data = sd->getPropertyByName("Restraint");
365
366 if (data != nullptr) {
367 // make sure we can reinterpret the generic data as restraint data:
368 std::shared_ptr<RestraintData> restData =
369 std::dynamic_pointer_cast<RestraintData>(data);
370 if (restData != nullptr) {
371 // make sure we can reinterpet the restraint data as a pointer to
372 // an ObjectRestraint:
373 ObjectRestraint* oRest =
374 dynamic_cast<ObjectRestraint*>(restData->getData());
375 if (oRest == NULL) {
376 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
377 "Can not cast RestraintData to ObjectRestraint\n");
378 painCave.severity = OPENMD_ERROR;
379 painCave.isFatal = 1;
380 simError();
381 } else {
382 if (sd->isDirectional()) {
383 oRest->setReferenceStructure(pos, q.toRotationMatrix3());
384 } else {
385 oRest->setReferenceStructure(pos);
386 }
387 }
388 }
389 }
390 }
391 }
392
393 void RestReader::readStuntDoubles(std::istream& inputStream) {
394 inputStream.getline(buffer, bufferSize);
395 std::string line(buffer);
396
397 if (line.find("<StuntDoubles>") == std::string::npos) {
398 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
399 "RestReader Error: Missing <StuntDoubles>\n");
400 painCave.isFatal = 1;
401 simError();
402 }
403
404 while (inputStream.getline(buffer, bufferSize)) {
405 line = buffer;
406
407 if (line.find("</StuntDoubles>") != std::string::npos) { break; }
408
409 parseDumpLine(line);
410 }
411 }
412
413 void RestReader::readFrameProperties(std::istream& inputStream) {
414 inputStream.getline(buffer, bufferSize);
415 std::string line(buffer);
416
417 if (line.find("<FrameData>") == std::string::npos) {
418 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
419 "RestReader Error: Missing <FrameData>\n");
420 painCave.isFatal = 1;
421 simError();
422 }
423
424 // restraints don't care about frame data (unless we need to wrap
425 // coordinates, but we'll worry about that later), so
426 // we'll just scan ahead until the end of the frame data:
427
428 while (inputStream.getline(buffer, bufferSize)) {
429 line = buffer;
430
431 if (line.find("</FrameData>") != std::string::npos) { break; }
432 }
433 }
434
435} // namespace OpenMD
int getNGlobalIntegrableObjects()
Returns the total number of integrable objects (total number of rigid bodies plus the total number of...
Definition SimInfo.hpp:139
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
Definition SimInfo.cpp:240
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
Definition SimInfo.cpp:245
StuntDouble * getIOIndexToIntegrableObject(int index)
return an integral objects by its global index.
Definition SimInfo.cpp:1027
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.