ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/io/RestReader.cpp

Comparing trunk/src/io/RestReader.cpp (file contents):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1796 by gezelter, Mon Sep 10 18:38:44 2012 UTC

#	Line 1 | Line 1
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
#	Line 36 | Line 36
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).
39	>	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	>	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42	<
43	<	#include "io/DumpReader.hpp"
44	<	#include "io/RestReader.hpp"
45	<	#include "primitives/Molecule.hpp"
42	>
43	>
44	>	#include <sys/types.h>
45	>	#include <sys/stat.h>
46	>
47	>	#include <math.h>
48	>	#include <string>
49	>	#include <sstream>
50	>	#include <iostream>
51	>	#include <stdio.h>
52	>	#include <stdlib.h>
53	>	#include <string.h>
54	>
55	>	#include "io/RestReader.hpp"
56	>	#include "primitives/Molecule.hpp"
57	>	#include "utils/simError.h"
58	>	#include "utils/MemoryUtils.hpp"
59	>	#include "utils/StringTokenizer.hpp"
60		#include "restraints/ObjectRestraint.hpp"
61	<	#include "restraints/MolecularRestraint.hpp"
61	>	#include "restraints/MolecularRestraint.hpp"
62	>
63	>	#ifdef IS_MPI
64	>
65	>	#include <mpi.h>
66	>	#endif
67
68		namespace OpenMD {
69	+
70	+	void RestReader::scanFile(){
71	+	int lineNo = 0;
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	<	void RestReader::readReferenceStructure() {
84	>	bool foundOpenSnapshotTag = false;
85	>
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_ = prevPos;
94	>	}
95	>	prevPos = currPos;
96	>	}
97	>
98	>	#ifdef IS_MPI
99	>	}
100	>	MPI::COMM_WORLD.Bcast(&framePos_, 1, MPI::INT, 0);
101	>	#endif // is_mpi
102	>	}
103
52	–	// some of this comes directly from DumpReader.
104
105	<	if (!isScanned_)
106	<	scanFile();
107	<
108	<	int storageLayout = info_->getSnapshotManager()->getStorageLayout();
109	<
110	<	if (storageLayout & DataStorage::dslPosition) {
111	<	needPos_ = true;
105	>	void RestReader::readSet(){
106	>	std::string line;
107	>
108	>	#ifndef IS_MPI
109	>
110	>	inFile_->clear();
111	>	inFile_->seekg(framePos_);
112	>
113	>	std::istream& inputStream = *inFile_;
114	>	#else
115	>
116	>	int masterNode = 0;
117	>	std::stringstream sstream;
118	>	if (worldRank == masterNode) {
119	>	std::string sendBuffer;
120	>
121	>	inFile_->clear();
122	>	inFile_->seekg(framePos_);
123	>
124	>	while (inFile_->getline(buffer, bufferSize)) {
125	>
126	>	line = buffer;
127	>	sendBuffer += line;
128	>	sendBuffer += '\n';
129	>	if (line.find("</Snapshot>") != std::string::npos) {
130	>	break;
131	>	}
132	>	}
133	>
134	>	int sendBufferSize = sendBuffer.size();
135	>	MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode);
136	>	MPI::COMM_WORLD.Bcast((void *)sendBuffer.c_str(), sendBufferSize,
137	>	MPI::CHAR, masterNode);
138	>
139	>	sstream.str(sendBuffer);
140		} else {
141	<	needPos_ = false;
141	>	int sendBufferSize;
142	>	MPI::COMM_WORLD.Bcast(&sendBufferSize, 1, MPI::INT, masterNode);
143	>	char * recvBuffer = new char[sendBufferSize+1];
144	>	assert(recvBuffer);
145	>	recvBuffer[sendBufferSize] = '\0';
146	>	MPI::COMM_WORLD.Bcast(recvBuffer, sendBufferSize, MPI::CHAR, masterNode);
147	>	sstream.str(recvBuffer);
148	>	delete [] recvBuffer;
149	>	}
150	>
151	>	std::istream& inputStream = sstream;
152	>	#endif
153	>
154	>	inputStream.getline(buffer, bufferSize);
155	>
156	>	line = buffer;
157	>	if (line.find("<Snapshot>") == std::string::npos) {
158	>	sprintf(painCave.errMsg,
159	>	"RestReader Error: can not find <Snapshot>\n");
160	>	painCave.isFatal = 1;
161	>	simError();
162		}
163	<
164	<	needVel_ = false;
165	<
166	<	if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
167	<	needQuaternion_ = true;
168	<	} else {
169	<	needQuaternion_ = false;
170	<	}
163	>
164	>	//read frameData
165	>	readFrameProperties(inputStream);
166	>
167	>	//read StuntDoubles
168	>	readStuntDoubles(inputStream);
169	>
170	>	inputStream.getline(buffer, bufferSize);
171	>	line = buffer;
172	>	if (line.find("</Snapshot>") == std::string::npos) {
173	>	sprintf(painCave.errMsg,
174	>	"RestReader Error: can not find </Snapshot>\n");
175	>	painCave.isFatal = 1;
176	>	simError();
177	>	}
178	>	}
179	>
180	>	void RestReader::readReferenceStructure() {
181
73	–	needAngMom_ = false;
74	–
182		// We need temporary storage to keep track of all StuntDouble positions
183		// in case some of the restraints are molecular (i.e. if they use
184		// multiple SD positions to determine restrained orientations or positions:
185
186		all_pos_.clear();
187	<	all_pos_.resize(info_->getNGlobalIntegrableObjects() );
188	<
82	<
187	>	all_pos_.resize(info_->getNGlobalIntegrableObjects()) ;
188	>
189		// Restraint files are just standard dump files, but with the reference
190		// structure stored in the first frame (frame 0).
191		// RestReader overloads readSet and explicitly handles all of the
192		// ObjectRestraints in that method:
193
194	<	readSet(0);
195	<
194	>	scanFile();
195	>
196	>	readSet();
197	>
198	>
199		// all ObjectRestraints have been handled, now we have to worry about
200		// molecular restraints:
201
#	Line 94 | Line 203 | namespace OpenMD {
203		Molecule::IntegrableObjectIterator j;
204		Molecule* mol;
205		StuntDouble* sd;
206	<
206	>
207		// no need to worry about parallel molecules, as molecules are not
208		// split across processor boundaries.  Just loop over all molecules
209		// we know about:
210	<
210	>
211		for (mol = info_->beginMolecule(i); mol != NULL;
212		mol = info_->nextMolecule(i)) {
213	<
213	>
214		// is this molecule restrained?
215		GenericData* data = mol->getPropertyByName("Restraint");
216
217		if (data != NULL) {
218	<
218	>
219		// make sure we can reinterpret the generic data as restraint data:
220	<
220	>
221		RestraintData* restData= dynamic_cast<RestraintData*>(data);
222	<
222	>
223		if (restData != NULL) {
224	<
224	>
225		// make sure we can reinterpet the restraint data as a
226		// pointer to a MolecularRestraint:
227	<
227	>
228		MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
229	<
229	>
230		if (mRest != NULL) {
231	<
231	>
232		// now we need to pack the stunt doubles for the reference
233		// structure:
234	<
234	>
235		std::vector<Vector3d> ref;
236		int count = 0;
237		RealType mass, mTot;
#	Line 140 | Line 249 | namespace OpenMD {
249		// doubles from the *globally* sorted array of
250		// positions:
251
252	<	ref.push_back( all_pos_[sd->getGlobalIndex()] );
253	<
145	<	mass = sd->getMass();
146	<
252	>	ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] );
253	>	mass = sd->getMass();
254		COM = COM + mass * ref[count];
255		mTot = mTot + mass;
256		count = count + 1;
257		}
151	–
258		COM /= mTot;
153	–
259		mRest->setReferenceStructure(ref, COM);
155	–
260		}
261		}
262		}
#	Line 162 | Line 266 | namespace OpenMD {
266
267
268		void RestReader::parseDumpLine(const std::string& line) {
269	+
270		StringTokenizer tokenizer(line);
271		int nTokens;
272
#	Line 169 | Line 274 | namespace OpenMD {
274
275		if (nTokens < 2) {
276		sprintf(painCave.errMsg,
277	<	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
277	>	"RestReader Error: Not enough Tokens.\n%s\n", line.c_str());
278		painCave.isFatal = 1;
279		simError();
280		}
281
282		int index = tokenizer.nextTokenAsInt();
178	–
179	–	StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
283
284	<	if (integrableObject == NULL) {
182	<	return;
183	<	}
284	>	StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
285
286	+	if (sd == NULL) {
287	+	return;
288	+	}
289	+
290		std::string type = tokenizer.nextToken();
291		int size = type.size();
292	+
293		Vector3d pos;
188	–	Vector3d vel;
294		Quat4d q;
295	<	Vector3d ji;
191	<	Vector3d force;
192	<	Vector3d torque;
193	<
295	>
296		for(int i = 0; i < size; ++i) {
297		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	–	}
298
299	<	case 'q' : {
300	<	if (integrableObject->isDirectional()) {
301	<
302	<	q[0] = tokenizer.nextTokenAsDouble();
303	<	q[1] = tokenizer.nextTokenAsDouble();
304	<	q[2] = tokenizer.nextTokenAsDouble();
305	<	q[3] = tokenizer.nextTokenAsDouble();
306	<
307	<	RealType qlen = q.length();
308	<	if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
309	<
310	<	sprintf(painCave.errMsg,
311	<	"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();
299	>	case 'p': {
300	>	pos[0] = tokenizer.nextTokenAsDouble();
301	>	pos[1] = tokenizer.nextTokenAsDouble();
302	>	pos[2] = tokenizer.nextTokenAsDouble();
303	>	break;
304	>	}
305	>	case 'v' : {
306	>	Vector3d vel;
307	>	vel[0] = tokenizer.nextTokenAsDouble();
308	>	vel[1] = tokenizer.nextTokenAsDouble();
309	>	vel[2] = tokenizer.nextTokenAsDouble();
310	>	break;
311	>	}
312
313	<	break;
313	>	case 'q' : {
314	>	if (sd->isDirectional()) {
315	>
316	>	q[0] = tokenizer.nextTokenAsDouble();
317	>	q[1] = tokenizer.nextTokenAsDouble();
318	>	q[2] = tokenizer.nextTokenAsDouble();
319	>	q[3] = tokenizer.nextTokenAsDouble();
320	>
321	>	RealType qlen = q.length();
322	>	if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
323	>
324	>	sprintf(painCave.errMsg,
325	>	"RestReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
326	>	painCave.isFatal = 1;
327	>	simError();
328	>	}
329	>
330	>	q.normalize();
331	>	}
332	>	break;
333	>	}
334	>	case 'j' : {
335	>	Vector3d ji;
336	>	if (sd->isDirectional()) {
337	>	ji[0] = tokenizer.nextTokenAsDouble();
338	>	ji[1] = tokenizer.nextTokenAsDouble();
339	>	ji[2] = tokenizer.nextTokenAsDouble();
340		}
341	<	case 't' : {
342	<	torque[0] = tokenizer.nextTokenAsDouble();
343	<	torque[1] = tokenizer.nextTokenAsDouble();
344	<	torque[2] = tokenizer.nextTokenAsDouble();
341	>	break;
342	>	}
343	>	case 'f': {
344	>	Vector3d force;
345	>	force[0] = tokenizer.nextTokenAsDouble();
346	>	force[1] = tokenizer.nextTokenAsDouble();
347	>	force[2] = tokenizer.nextTokenAsDouble();
348	>	break;
349	>	}
350	>	case 't' : {
351	>	Vector3d torque;
352	>	torque[0] = tokenizer.nextTokenAsDouble();
353	>	torque[1] = tokenizer.nextTokenAsDouble();
354	>	torque[2] = tokenizer.nextTokenAsDouble();
355	>	break;
356	>	}
357	>	default: {
358	>	sprintf(painCave.errMsg,
359	>	"RestReader Error: %s is an unrecognized type\n", type.c_str());
360	>	painCave.isFatal = 1;
361	>	simError();
362	>	break;
363	>	}
364	>	}
365	>	// keep the position in case we need it for a molecular restraint:
366
367	<	break;
367	>	all_pos_[index] = pos;
368	>
369	>	// is this io restrained?
370	>	GenericData* data = sd->getPropertyByName("Restraint");
371	>	ObjectRestraint* oRest;
372	>
373	>	if (data != NULL) {
374	>	// make sure we can reinterpret the generic data as restraint data:
375	>	RestraintData* restData= dynamic_cast<RestraintData*>(data);
376	>	if (restData != NULL) {
377	>	// make sure we can reinterpet the restraint data as a pointer to
378	>	// an ObjectRestraint:
379	>	oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
380	>	if (oRest != NULL) {
381	>	if (sd->isDirectional()) {
382	>	oRest->setReferenceStructure(pos, q.toRotationMatrix3());
383	>	} else {
384	>	oRest->setReferenceStructure(pos);
385	>	}
386	>	}
387		}
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	–
388		}
389		}
390	<
391	<	// keep the position in case we need it for a molecular restraint:
390	>	}
391	>
392	>	void  RestReader::readStuntDoubles(std::istream& inputStream) {
393
394	<	all_pos_[index] = pos;
395	<
396	<	// is this io restrained?
397	<	GenericData* data = integrableObject->getPropertyByName("Restraint");
398	<	ObjectRestraint* oRest;
399	<
400	<	if (data != NULL) {
401	<	// make sure we can reinterpret the generic data as restraint data:
402	<	RestraintData* restData= dynamic_cast<RestraintData*>(data);
403	<	if (restData != NULL) {
404	<	// make sure we can reinterpet the restraint data as a pointer to
405	<	// an ObjectRestraint:
406	<	oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
407	<	if (oRest != NULL) {
408	<	if (integrableObject->isDirectional()) {
282	<	oRest->setReferenceStructure(pos, q.toRotationMatrix3());
283	<	} else {
284	<	oRest->setReferenceStructure(pos);
285	<	}
286	<	}
394	>	inputStream.getline(buffer, bufferSize);
395	>	std::string line(buffer);
396	>
397	>	if (line.find("<StuntDoubles>") == std::string::npos) {
398	>	sprintf(painCave.errMsg,
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) {
408	>	break;
409		}
410	+
411	+	parseDumpLine(line);
412		}
413	<
413	>
414		}
291	–
415
416	<
416	>
417		void RestReader::readFrameProperties(std::istream& inputStream) {
418		inputStream.getline(buffer, bufferSize);
419		std::string line(buffer);
420
421		if (line.find("<FrameData>") == std::string::npos) {
422		sprintf(painCave.errMsg,
423	<	"DumpReader Error: Missing <FrameData>\n");
423	>	"RestReader Error: Missing <FrameData>\n");
424		painCave.isFatal = 1;
425		simError();
426		}
#	Line 314 | Line 437 | namespace OpenMD {
437		}
438
439		}
440	<
440	>
441		}
442
443

Comparing trunk/src/io/RestReader.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1796 by gezelter, Mon Sep 10 18:38:44 2012 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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