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 1971 by gezelter, Fri Feb 28 13:25:13 2014 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 35 | Line 35
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).
38	>	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
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	>	#ifdef IS_MPI
44	>	#include <mpi.h>
45	>	#endif
46	>
47	>	#include <sys/types.h>
48	>	#include <sys/stat.h>
49	>
50	>	#include <math.h>
51	>	#include <string>
52	>	#include <sstream>
53	>	#include <iostream>
54	>	#include <stdio.h>
55	>	#include <stdlib.h>
56	>	#include <string.h>
57	>
58	>	#include "io/RestReader.hpp"
59	>	#include "primitives/Molecule.hpp"
60	>	#include "utils/simError.h"
61	>	#include "utils/MemoryUtils.hpp"
62	>	#include "utils/StringTokenizer.hpp"
63		#include "restraints/ObjectRestraint.hpp"
64	<	#include "restraints/MolecularRestraint.hpp"
64	>	#include "restraints/MolecularRestraint.hpp"
65
66		namespace OpenMD {
67	+
68	+	void RestReader::scanFile(){
69	+
70	+	std::streampos prevPos;
71	+	std::streampos currPos;
72	+
73	+	#ifdef IS_MPI
74	+
75	+	if (worldRank == 0) {
76	+	#endif // is_mpi
77	+
78	+	inFile_->clear();
79	+	currPos = inFile_->tellg();
80	+	prevPos = currPos;
81
82	<	void RestReader::readReferenceStructure() {
82	>	bool foundOpenSnapshotTag = false;
83	>	int lineNo = 0;
84	>	while(!foundOpenSnapshotTag && inFile_->getline(buffer, bufferSize)) {
85	>	++lineNo;
86	>
87	>	std::string line = buffer;
88	>	currPos = inFile_->tellg();
89	>	if (line.find("<Snapshot>")!= std::string::npos) {
90	>	foundOpenSnapshotTag = true;
91	>	framePos_ = prevPos;
92	>	}
93	>	prevPos = currPos;
94	>	}
95	>
96	>	#ifdef IS_MPI
97	>	}
98	>	MPI_Bcast(&framePos_, 1, MPI_INT, 0, MPI_COMM_WORLD);
99	>	#endif // is_mpi
100	>	}
101
52	–	// some of this comes directly from DumpReader.
102
103	<	if (!isScanned_)
104	<	scanFile();
105	<
106	<	int storageLayout = info_->getSnapshotManager()->getStorageLayout();
107	<
108	<	if (storageLayout & DataStorage::dslPosition) {
109	<	needPos_ = true;
103	>	void RestReader::readSet(){
104	>	std::string line;
105	>
106	>	#ifndef IS_MPI
107	>
108	>	inFile_->clear();
109	>	inFile_->seekg(framePos_);
110	>
111	>	std::istream& inputStream = *inFile_;
112	>	#else
113	>
114	>	int masterNode = 0;
115	>	std::stringstream sstream;
116	>	if (worldRank == masterNode) {
117	>	std::string sendBuffer;
118	>
119	>	inFile_->clear();
120	>	inFile_->seekg(framePos_);
121	>
122	>	while (inFile_->getline(buffer, bufferSize)) {
123	>
124	>	line = buffer;
125	>	sendBuffer += line;
126	>	sendBuffer += '\n';
127	>	if (line.find("</Snapshot>") != std::string::npos) {
128	>	break;
129	>	}
130	>	}
131	>
132	>	int sendBufferSize = sendBuffer.size();
133	>	MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
134	>	MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize,
135	>	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,
145	>	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	>
371	>	if (data != NULL) {
372	>	// make sure we can reinterpret the generic data as restraint data:
373	>	RestraintData* restData= dynamic_cast<RestraintData*>(data);
374	>	if (restData != NULL) {
375	>	// make sure we can reinterpet the restraint data as a pointer to
376	>	// an ObjectRestraint:
377	>	ObjectRestraint* oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
378	>	if (oRest != NULL) {
379	>	if (sd->isDirectional()) {
380	>	oRest->setReferenceStructure(pos, q.toRotationMatrix3());
381	>	} else {
382	>	oRest->setReferenceStructure(pos);
383	>	}
384	>	}
385		}
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	–
386		}
387		}
388	<
389	<	// keep the position in case we need it for a molecular restraint:
388	>	}
389	>
390	>	void  RestReader::readStuntDoubles(std::istream& inputStream) {
391
392	<	all_pos_[index] = pos;
393	<
394	<	// is this io restrained?
395	<	GenericData* data = integrableObject->getPropertyByName("Restraint");
396	<	ObjectRestraint* oRest;
397	<
398	<	if (data != NULL) {
399	<	// make sure we can reinterpret the generic data as restraint data:
400	<	RestraintData* restData= dynamic_cast<RestraintData*>(data);
401	<	if (restData != NULL) {
402	<	// make sure we can reinterpet the restraint data as a pointer to
403	<	// an ObjectRestraint:
404	<	oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
405	<	if (oRest != NULL) {
406	<	if (integrableObject->isDirectional()) {
282	<	oRest->setReferenceStructure(pos, q.toRotationMatrix3());
283	<	} else {
284	<	oRest->setReferenceStructure(pos);
285	<	}
286	<	}
392	>	inputStream.getline(buffer, bufferSize);
393	>	std::string line(buffer);
394	>
395	>	if (line.find("<StuntDoubles>") == std::string::npos) {
396	>	sprintf(painCave.errMsg,
397	>	"RestReader Error: Missing <StuntDoubles>\n");
398	>	painCave.isFatal = 1;
399	>	simError();
400	>	}
401	>
402	>	while(inputStream.getline(buffer, bufferSize)) {
403	>	line = buffer;
404	>
405	>	if(line.find("</StuntDoubles>") != std::string::npos) {
406	>	break;
407		}
408	+
409	+	parseDumpLine(line);
410		}
411	<
411	>
412		}
291	–
413
414	<
414	>
415		void RestReader::readFrameProperties(std::istream& inputStream) {
416		inputStream.getline(buffer, bufferSize);
417		std::string line(buffer);
418
419		if (line.find("<FrameData>") == std::string::npos) {
420		sprintf(painCave.errMsg,
421	<	"DumpReader Error: Missing <FrameData>\n");
421	>	"RestReader Error: Missing <FrameData>\n");
422		painCave.isFatal = 1;
423		simError();
424		}
#	Line 314 | Line 435 | namespace OpenMD {
435		}
436
437		}
438	<
438	>
439		}
440
441

Comparing trunk/src/io/RestReader.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1971 by gezelter, Fri Feb 28 13:25:13 2014 UTC

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

Diff Legend

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