ViewVC Help

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

root/OpenMD/branches/development/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.
branches/development/src/io/RestReader.cpp (file contents), Revision 1769 by gezelter, Mon Jul 9 14:15:52 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_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD);
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_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
136	>	MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
137	>
138	>	sstream.str(sendBuffer);
139		} else {
140	<	needPos_ = false;
140	>	int sendBufferSize;
141	>	MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
142	>	char * recvBuffer = new char[sendBufferSize+1];
143	>	assert(recvBuffer);
144	>	recvBuffer[sendBufferSize] = '\0';
145	>	MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
146	>	sstream.str(recvBuffer);
147	>	delete [] recvBuffer;
148	>	}
149	>
150	>	std::istream& inputStream = sstream;
151	>	#endif
152	>
153	>	inputStream.getline(buffer, bufferSize);
154	>
155	>	line = buffer;
156	>	if (line.find("<Snapshot>") == std::string::npos) {
157	>	sprintf(painCave.errMsg,
158	>	"RestReader Error: can not find <Snapshot>\n");
159	>	painCave.isFatal = 1;
160	>	simError();
161		}
162	<
163	<	needVel_ = false;
164	<
165	<	if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
166	<	needQuaternion_ = true;
167	<	} else {
168	<	needQuaternion_ = false;
169	<	}
162	>
163	>	//read frameData
164	>	readFrameProperties(inputStream);
165	>
166	>	//read StuntDoubles
167	>	readStuntDoubles(inputStream);
168	>
169	>	inputStream.getline(buffer, bufferSize);
170	>	line = buffer;
171	>	if (line.find("</Snapshot>") == std::string::npos) {
172	>	sprintf(painCave.errMsg,
173	>	"RestReader Error: can not find </Snapshot>\n");
174	>	painCave.isFatal = 1;
175	>	simError();
176	>	}
177	>	}
178	>
179	>	void RestReader::readReferenceStructure() {
180
73	–	needAngMom_ = false;
74	–
181		// We need temporary storage to keep track of all StuntDouble positions
182		// in case some of the restraints are molecular (i.e. if they use
183		// multiple SD positions to determine restrained orientations or positions:
184
185		all_pos_.clear();
186	<	all_pos_.resize(info_->getNGlobalIntegrableObjects() );
187	<
82	<
186	>	all_pos_.resize(info_->getNGlobalIntegrableObjects()) ;
187	>
188		// Restraint files are just standard dump files, but with the reference
189		// structure stored in the first frame (frame 0).
190		// RestReader overloads readSet and explicitly handles all of the
191		// ObjectRestraints in that method:
192
193	<	readSet(0);
194	<
193	>	scanFile();
194	>
195	>	readSet();
196	>
197	>
198		// all ObjectRestraints have been handled, now we have to worry about
199		// molecular restraints:
200
#	Line 94 | Line 202 | namespace OpenMD {
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;
#	Line 140 | Line 248 | namespace OpenMD {
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		}
#	Line 162 | Line 265 | namespace OpenMD {
265
266
267		void RestReader::parseDumpLine(const std::string& line) {
268	+
269		StringTokenizer tokenizer(line);
270		int nTokens;
271
#	Line 169 | Line 273 | namespace OpenMD {
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) {
182	<	return;
183	<	}
283	>	StuntDouble* sd = info_->getIOIndexToIntegrableObject(index);
284
285	+	if (sd == NULL) {
286	+	return;
287	+	}
288	+
289		std::string type = tokenizer.nextToken();
290		int size = type.size();
291	+
292		Vector3d pos;
188	–	Vector3d vel;
293		Quat4d q;
294	<	Vector3d ji;
191	<	Vector3d force;
192	<	Vector3d torque;
193	<
294	>
295		for(int i = 0; i < size; ++i) {
296		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	–	}
297
298	<	case 'q' : {
299	<	if (integrableObject->isDirectional()) {
300	<
301	<	q[0] = tokenizer.nextTokenAsDouble();
302	<	q[1] = tokenizer.nextTokenAsDouble();
303	<	q[2] = tokenizer.nextTokenAsDouble();
304	<	q[3] = tokenizer.nextTokenAsDouble();
305	<
306	<	RealType qlen = q.length();
307	<	if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
308	<
309	<	sprintf(painCave.errMsg,
310	<	"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();
298	>	case 'p': {
299	>	pos[0] = tokenizer.nextTokenAsDouble();
300	>	pos[1] = tokenizer.nextTokenAsDouble();
301	>	pos[2] = tokenizer.nextTokenAsDouble();
302	>	break;
303	>	}
304	>	case 'v' : {
305	>	Vector3d vel;
306	>	vel[0] = tokenizer.nextTokenAsDouble();
307	>	vel[1] = tokenizer.nextTokenAsDouble();
308	>	vel[2] = tokenizer.nextTokenAsDouble();
309	>	break;
310	>	}
311
312	<	break;
312	>	case 'q' : {
313	>	if (sd->isDirectional()) {
314	>
315	>	q[0] = tokenizer.nextTokenAsDouble();
316	>	q[1] = tokenizer.nextTokenAsDouble();
317	>	q[2] = tokenizer.nextTokenAsDouble();
318	>	q[3] = tokenizer.nextTokenAsDouble();
319	>
320	>	RealType qlen = q.length();
321	>	if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
322	>
323	>	sprintf(painCave.errMsg,
324	>	"RestReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
325	>	painCave.isFatal = 1;
326	>	simError();
327	>	}
328	>
329	>	q.normalize();
330	>	}
331	>	break;
332	>	}
333	>	case 'j' : {
334	>	Vector3d ji;
335	>	if (sd->isDirectional()) {
336	>	ji[0] = tokenizer.nextTokenAsDouble();
337	>	ji[1] = tokenizer.nextTokenAsDouble();
338	>	ji[2] = tokenizer.nextTokenAsDouble();
339		}
340	<	case 't' : {
341	<	torque[0] = tokenizer.nextTokenAsDouble();
342	<	torque[1] = tokenizer.nextTokenAsDouble();
343	<	torque[2] = tokenizer.nextTokenAsDouble();
340	>	break;
341	>	}
342	>	case 'f': {
343	>	Vector3d force;
344	>	force[0] = tokenizer.nextTokenAsDouble();
345	>	force[1] = tokenizer.nextTokenAsDouble();
346	>	force[2] = tokenizer.nextTokenAsDouble();
347	>	break;
348	>	}
349	>	case 't' : {
350	>	Vector3d torque;
351	>	torque[0] = tokenizer.nextTokenAsDouble();
352	>	torque[1] = tokenizer.nextTokenAsDouble();
353	>	torque[2] = tokenizer.nextTokenAsDouble();
354	>	break;
355	>	}
356	>	default: {
357	>	sprintf(painCave.errMsg,
358	>	"RestReader Error: %s is an unrecognized type\n", type.c_str());
359	>	painCave.isFatal = 1;
360	>	simError();
361	>	break;
362	>	}
363	>	}
364	>	// keep the position in case we need it for a molecular restraint:
365
366	<	break;
366	>	all_pos_[index] = pos;
367	>
368	>	// is this io restrained?
369	>	GenericData* data = sd->getPropertyByName("Restraint");
370	>	ObjectRestraint* oRest;
371	>
372	>	if (data != NULL) {
373	>	// make sure we can reinterpret the generic data as restraint data:
374	>	RestraintData* restData= dynamic_cast<RestraintData*>(data);
375	>	if (restData != NULL) {
376	>	// make sure we can reinterpet the restraint data as a pointer to
377	>	// an ObjectRestraint:
378	>	oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
379	>	if (oRest != NULL) {
380	>	if (sd->isDirectional()) {
381	>	oRest->setReferenceStructure(pos, q.toRotationMatrix3());
382	>	} else {
383	>	oRest->setReferenceStructure(pos);
384	>	}
385	>	}
386		}
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	–
387		}
388		}
389	<
390	<	// keep the position in case we need it for a molecular restraint:
389	>	}
390	>
391	>	void  RestReader::readStuntDoubles(std::istream& inputStream) {
392
393	<	all_pos_[index] = pos;
394	<
395	<	// is this io restrained?
396	<	GenericData* data = integrableObject->getPropertyByName("Restraint");
397	<	ObjectRestraint* oRest;
398	<
399	<	if (data != NULL) {
400	<	// make sure we can reinterpret the generic data as restraint data:
401	<	RestraintData* restData= dynamic_cast<RestraintData*>(data);
402	<	if (restData != NULL) {
403	<	// make sure we can reinterpet the restraint data as a pointer to
404	<	// an ObjectRestraint:
405	<	oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
406	<	if (oRest != NULL) {
407	<	if (integrableObject->isDirectional()) {
282	<	oRest->setReferenceStructure(pos, q.toRotationMatrix3());
283	<	} else {
284	<	oRest->setReferenceStructure(pos);
285	<	}
286	<	}
393	>	inputStream.getline(buffer, bufferSize);
394	>	std::string line(buffer);
395	>
396	>	if (line.find("<StuntDoubles>") == std::string::npos) {
397	>	sprintf(painCave.errMsg,
398	>	"RestReader Error: Missing <StuntDoubles>\n");
399	>	painCave.isFatal = 1;
400	>	simError();
401	>	}
402	>
403	>	while(inputStream.getline(buffer, bufferSize)) {
404	>	line = buffer;
405	>
406	>	if(line.find("</StuntDoubles>") != std::string::npos) {
407	>	break;
408		}
409	+
410	+	parseDumpLine(line);
411		}
412	<
412	>
413		}
291	–
414
415	<
415	>
416		void RestReader::readFrameProperties(std::istream& inputStream) {
417		inputStream.getline(buffer, bufferSize);
418		std::string line(buffer);
419
420		if (line.find("<FrameData>") == std::string::npos) {
421		sprintf(painCave.errMsg,
422	<	"DumpReader Error: Missing <FrameData>\n");
422	>	"RestReader Error: Missing <FrameData>\n");
423		painCave.isFatal = 1;
424		simError();
425		}
#	Line 314 | Line 436 | namespace OpenMD {
436		}
437
438		}
439	<
439	>
440		}
441
442

Comparing:
trunk/src/io/RestReader.cpp (property svn:keywords), Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
branches/development/src/io/RestReader.cpp (property svn:keywords), Revision 1769 by gezelter, Mon Jul 9 14:15:52 2012 UTC

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

Diff Legend

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