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 1360 by cli2, Mon Sep 7 16:31:51 2009 UTC vs.
branches/development/src/io/RestReader.cpp (file contents), Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 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
6	<	* redistribute this software in source and binary code form, provided
7	<	* that the following conditions are met:
8	<	*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
19	<	*    notice, this list of conditions and the following disclaimer.
20	<	*
21	<	* 3. Redistributions in binary form must reproduce the above copyright
22	<	*    notice, this list of conditions and the following disclaimer in the
23	<	*    documentation and/or other materials provided with the
24	<	*    distribution.
25	<	*
26	<	* This software is provided "AS IS," without a warranty of any
27	<	* kind. All express or implied conditions, representations and
28	<	* warranties, including any implied warranty of merchantability,
29	<	* fitness for a particular purpose or non-infringement, are hereby
30	<	* excluded.  The University of Notre Dame and its licensors shall not
31	<	* be liable for any damages suffered by licensee as a result of
32	<	* using, modifying or distributing the software or its
33	<	* derivatives. In no event will the University of Notre Dame or its
34	<	* licensors be liable for any lost revenue, profit or data, or for
35	<	* direct, indirect, special, consequential, incidental or punitive
36	<	* damages, however caused and regardless of the theory of liability,
37	<	* arising out of the use of or inability to use software, even if the
38	<	* University of Notre Dame has been advised of the possibility of
39	<	* such damages.
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	<	#include "io/DumpReader.hpp"
43	<	#include "io/RestReader.hpp"
44	<	#include "primitives/Molecule.hpp"
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"
60	>	#include "restraints/MolecularRestraint.hpp"
61	>
62	>	#ifdef IS_MPI
63	>
64	>	#include <mpi.h>
65	>	#endif
66
67	<	namespace oopse {
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	<	void RestReader::readReferenceStructure() {
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
52	–	// some of this comes directly from DumpReader.
103
104	<	if (!isScanned_)
105	<	scanFile();
106	<
107	<	int storageLayout = info_->getSnapshotManager()->getStorageLayout();
108	<
109	<	if (storageLayout & DataStorage::dslPosition) {
110	<	needPos_ = true;
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	<	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, 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	<	needVel_ = false;
163	<
164	<	if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
165	<	needQuaternion_ = true;
166	<	} else {
167	<	needQuaternion_ = false;
168	<	}
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
73	–	needAngMom_ = false;
74	–
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	<
82	<
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	<	readSet(0);
193	<
192	>	scanFile();
193	>
194	>	readSet();
195	>
196	>
197		// all ObjectRestraints have been handled, now we have to worry about
198		// molecular restraints:
199
#	Line 94 | Line 201 | namespace oopse {
201		Molecule::IntegrableObjectIterator j;
202		Molecule* mol;
203		StuntDouble* sd;
204	<
204	>
205		// no need to worry about parallel molecules, as molecules are not
206		// split across processor boundaries.  Just loop over all molecules
207		// we know about:
208	<
208	>
209		for (mol = info_->beginMolecule(i); mol != NULL;
210		mol = info_->nextMolecule(i)) {
211	<
211	>
212		// is this molecule restrained?
213		GenericData* data = mol->getPropertyByName("Restraint");
214
215		if (data != NULL) {
216	<
216	>
217		// make sure we can reinterpret the generic data as restraint data:
218	<
219	<	RestraintData* restData= dynamic_cast<RestraintData*>(data);
220	<
218	>
219	>	RestraintData* restData= dynamic_cast<RestraintData*>(data);
220	>
221		if (restData != NULL) {
222	<
222	>
223		// make sure we can reinterpet the restraint data as a
224		// pointer to a MolecularRestraint:
225	<
225	>
226		MolecularRestraint* mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
227	<
227	>
228		if (mRest != NULL) {
229	<
229	>
230		// now we need to pack the stunt doubles for the reference
231		// structure:
232	<
232	>
233		std::vector<Vector3d> ref;
234		int count = 0;
235		RealType mass, mTot;
#	Line 140 | Line 247 | namespace oopse {
247		// doubles from the *globally* sorted array of
248		// positions:
249
250	<	ref.push_back( all_pos_[sd->getGlobalIndex()] );
251	<
145	<	mass = sd->getMass();
146	<
250	>	ref.push_back( all_pos_[sd->getGlobalIntegrableObjectIndex()] );
251	>	mass = sd->getMass();
252		COM = COM + mass * ref[count];
253		mTot = mTot + mass;
254		count = count + 1;
255		}
151	–
256		COM /= mTot;
153	–
257		mRest->setReferenceStructure(ref, COM);
155	–
258		}
259		}
260		}
#	Line 162 | Line 264 | namespace oopse {
264
265
266		void RestReader::parseDumpLine(const std::string& line) {
267	+
268		StringTokenizer tokenizer(line);
269		int nTokens;
270
#	Line 169 | Line 272 | namespace oopse {
272
273		if (nTokens < 2) {
274		sprintf(painCave.errMsg,
275	<	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
275	>	"RestReader Error: Not enough Tokens.\n%s\n", line.c_str());
276		painCave.isFatal = 1;
277		simError();
278		}
279
280		int index = tokenizer.nextTokenAsInt();
281	<
281	>
282		StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
283
284		if (integrableObject == NULL) {
285		return;
286	<	}
287	<
286	>	}
287	>
288		std::string type = tokenizer.nextToken();
289		int size = type.size();
290	+
291		Vector3d pos;
188	–	Vector3d vel;
292		Quat4d q;
293	<	Vector3d ji;
191	<	Vector3d force;
192	<	Vector3d torque;
193	<
293	>
294		for(int i = 0; i < size; ++i) {
295		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	–	}
296
297	<	case 'q' : {
298	<	if (integrableObject->isDirectional()) {
299	<
300	<	q[0] = tokenizer.nextTokenAsDouble();
301	<	q[1] = tokenizer.nextTokenAsDouble();
302	<	q[2] = tokenizer.nextTokenAsDouble();
303	<	q[3] = tokenizer.nextTokenAsDouble();
304	<
305	<	RealType qlen = q.length();
306	<	if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
307	<
308	<	sprintf(painCave.errMsg,
309	<	"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();
297	>	case 'p': {
298	>	pos[0] = tokenizer.nextTokenAsDouble();
299	>	pos[1] = tokenizer.nextTokenAsDouble();
300	>	pos[2] = tokenizer.nextTokenAsDouble();
301	>	break;
302	>	}
303	>	case 'v' : {
304	>	Vector3d vel;
305	>	vel[0] = tokenizer.nextTokenAsDouble();
306	>	vel[1] = tokenizer.nextTokenAsDouble();
307	>	vel[2] = tokenizer.nextTokenAsDouble();
308	>	break;
309	>	}
310
311	<	break;
311	>	case 'q' : {
312	>	if (integrableObject->isDirectional()) {
313	>
314	>	q[0] = tokenizer.nextTokenAsDouble();
315	>	q[1] = tokenizer.nextTokenAsDouble();
316	>	q[2] = tokenizer.nextTokenAsDouble();
317	>	q[3] = tokenizer.nextTokenAsDouble();
318	>
319	>	RealType qlen = q.length();
320	>	if (qlen < OpenMD::epsilon) { //check quaternion is not equal to 0
321	>
322	>	sprintf(painCave.errMsg,
323	>	"RestReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
324	>	painCave.isFatal = 1;
325	>	simError();
326	>	}
327	>
328	>	q.normalize();
329	>	}
330	>	break;
331	>	}
332	>	case 'j' : {
333	>	Vector3d ji;
334	>	if (integrableObject->isDirectional()) {
335	>	ji[0] = tokenizer.nextTokenAsDouble();
336	>	ji[1] = tokenizer.nextTokenAsDouble();
337	>	ji[2] = tokenizer.nextTokenAsDouble();
338		}
339	<	case 't' : {
340	<	torque[0] = tokenizer.nextTokenAsDouble();
341	<	torque[1] = tokenizer.nextTokenAsDouble();
342	<	torque[2] = tokenizer.nextTokenAsDouble();
339	>	break;
340	>	}
341	>	case 'f': {
342	>	Vector3d force;
343	>	force[0] = tokenizer.nextTokenAsDouble();
344	>	force[1] = tokenizer.nextTokenAsDouble();
345	>	force[2] = tokenizer.nextTokenAsDouble();
346	>	break;
347	>	}
348	>	case 't' : {
349	>	Vector3d torque;
350	>	torque[0] = tokenizer.nextTokenAsDouble();
351	>	torque[1] = tokenizer.nextTokenAsDouble();
352	>	torque[2] = tokenizer.nextTokenAsDouble();
353	>	break;
354	>	}
355	>	default: {
356	>	sprintf(painCave.errMsg,
357	>	"RestReader Error: %s is an unrecognized type\n", type.c_str());
358	>	painCave.isFatal = 1;
359	>	simError();
360	>	break;
361	>	}
362	>	}
363	>	// keep the position in case we need it for a molecular restraint:
364
365	<	break;
365	>	all_pos_[index] = pos;
366	>
367	>	// is this io restrained?
368	>	GenericData* data = integrableObject->getPropertyByName("Restraint");
369	>	ObjectRestraint* oRest;
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	>	oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
378	>	if (oRest != NULL) {
379	>	if (integrableObject->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 oopse {
435		}
436
437		}
438	<
438	>
439		}
440
441
442	<	}//end namespace oopse
442	>	}//end namespace OpenMD

Comparing:
trunk/src/io/RestReader.cpp (property svn:keywords), Revision 1360 by cli2, Mon Sep 7 16:31:51 2009 UTC vs.
branches/development/src/io/RestReader.cpp (property svn:keywords), Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC

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

Diff Legend

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