ViewVC Help

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

root/OpenMD/trunk/src/applications/staticProps/P2OrderParameter.cpp

Comparing trunk/src/applications/staticProps/P2OrderParameter.cpp (file contents):
Revision 1078 by gezelter, Wed Oct 18 21:58:48 2006 UTC vs.
Revision 1542 by gezelter, Thu Mar 3 20:32:49 2011 UTC

#	Line 6 | Line 6
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
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
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.
#	Line 37 | Line 28
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 "applications/staticProps/P2OrderParameter.hpp"
#	Line 44 | Line 44
44		#include "io/DumpReader.hpp"
45		#include "primitives/Molecule.hpp"
46		#include "utils/NumericConstant.hpp"
47	–	namespace oopse {
47
48	+	using namespace std;
49	+	namespace OpenMD {
50
51	<	P2OrderParameter::P2OrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
52	<	: StaticAnalyser(info, filename),
53	<	selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
54	<	seleMan1_(info), seleMan2_(info){
51	>	P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename,
52	>	const string& sele1)
53	>	: StaticAnalyser(info, filename), doVect_(true),
54	>	selectionScript1_(sele1), evaluator1_(info),
55	>	evaluator2_(info), seleMan1_(info), seleMan2_(info) {
56	>
57	>	setOutputName(getPrefix(filename) + ".p2");
58	>
59	>	evaluator1_.loadScriptString(sele1);
60	>
61	>	if (!evaluator1_.isDynamic()) {
62	>	seleMan1_.setSelectionSet(evaluator1_.evaluate());
63	>	}
64	>
65	>	}
66
67	+	P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename,
68	+	const string& sele1, const string& sele2)
69	+	: StaticAnalyser(info, filename), doVect_(false),
70	+	selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info),
71	+	evaluator2_(info), seleMan1_(info), seleMan2_(info) {
72	+
73		setOutputName(getPrefix(filename) + ".p2");
74	<
74	>
75		evaluator1_.loadScriptString(sele1);
76		evaluator2_.loadScriptString(sele2);
77	<
77	>
78		if (!evaluator1_.isDynamic()) {
79		seleMan1_.setSelectionSet(evaluator1_.evaluate());
80		}else {
81		sprintf( painCave.errMsg,
82		"--sele1 must be static selection\n");
83	<	painCave.severity = OOPSE_ERROR;
83	>	painCave.severity = OPENMD_ERROR;
84		painCave.isFatal = 1;
85		simError();
86		}
87	<
87	>
88		if (!evaluator2_.isDynamic()) {
89		seleMan2_.setSelectionSet(evaluator2_.evaluate());
90		}else {
91		sprintf( painCave.errMsg,
92		"--sele2 must be static selection\n");
93	<	painCave.severity = OOPSE_ERROR;
93	>	painCave.severity = OPENMD_ERROR;
94		painCave.isFatal = 1;
95		simError();
96		}
97	<
97	>
98		if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
99		sprintf( painCave.errMsg,
100		"The number of selected Stuntdoubles are not the same in --sele1 and sele2\n");
101	<	painCave.severity = OOPSE_ERROR;
101	>	painCave.severity = OPENMD_ERROR;
102		painCave.isFatal = 1;
103		simError();
104	<
104	>
105		}
106	<
106	>
107		int i;
108		int j;
109		StuntDouble* sd1;
#	Line 94 | Line 112 | namespace oopse {
112		sd1 != NULL && sd2 != NULL;
113		sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
114
115	<	sdPairs_.push_back(std::make_pair(sd1, sd2));
116	<	}
99	<
100	<
115	>	sdPairs_.push_back(make_pair(sd1, sd2));
116	>	}
117		}
118
119		void P2OrderParameter::process() {
#	Line 105 | Line 121 | namespace oopse {
121		RigidBody* rb;
122		SimInfo::MoleculeIterator mi;
123		Molecule::RigidBodyIterator rbIter;
124	+	StuntDouble* sd;
125	+	int i;
126	+
127
128		DumpReader reader(info_, dumpFilename_);
129		int nFrames = reader.getNFrames();
#	Line 113 | Line 132 | namespace oopse {
132		reader.readFrame(i);
133		currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
134
135	<
136	<	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
135	>	for (mol = info_->beginMolecule(mi); mol != NULL;
136	>	mol = info_->nextMolecule(mi)) {
137		//change the positions of atoms which belong to the rigidbodies
138	<	for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
138	>	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
139	>	rb = mol->nextRigidBody(rbIter)) {
140		rb->updateAtoms();
141	<	}
122	<
141	>	}
142		}
143
144		Mat3x3d orderTensor(0.0);
145	<	for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
146	<	Vector3d vec = j->first->getPos() - j->second->getPos();
147	<	if (usePeriodicBoundaryConditions_)
148	<	currentSnapshot_->wrapVector(vec);
149	<	vec.normalize();
150	<	orderTensor +=outProduct(vec, vec);
145	>
146	>	if (doVect_) {
147	>
148	>	if  (evaluator1_.isDynamic())
149	>	seleMan1_.setSelectionSet(evaluator1_.evaluate());
150	>
151	>	for (sd = seleMan1_.beginSelected(i); sd != NULL;
152	>	sd = seleMan1_.nextSelected(i)) {
153	>	if (sd->isDirectional()) {
154	>	Vector3d vec = sd->getA().getColumn(2);
155	>	vec.normalize();
156	>	orderTensor += outProduct(vec, vec);
157	>	}
158	>	}
159	>
160	>	orderTensor /= seleMan1_.getSelectionCount();
161	>
162	>	} else {
163	>
164	>	for (vector<pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin();
165	>	j != sdPairs_.end(); ++j) {
166	>	Vector3d vec = j->first->getPos() - j->second->getPos();
167	>	if (usePeriodicBoundaryConditions_)
168	>	currentSnapshot_->wrapVector(vec);
169	>	vec.normalize();
170	>	orderTensor +=outProduct(vec, vec);
171	>	}
172	>
173	>	orderTensor /= sdPairs_.size();
174		}
175
134	–	orderTensor /= sdPairs_.size();
176		orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();
177
178		Vector3d eigenvalues;
#	Line 155 | Line 196 | namespace oopse {
196		}
197
198		RealType angle = 0.0;
158	–	for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
159	–	Vector3d vec = j->first->getPos() - j->second->getPos();
160	–	if (usePeriodicBoundaryConditions_)
161	–	currentSnapshot_->wrapVector(vec);
162	–	vec.normalize();
199
200	<	angle += acos(dot(vec, director)) ;
200	>
201	>	if (doVect_) {
202	>	for (sd = seleMan1_.beginSelected(i); sd != NULL;
203	>	sd = seleMan1_.nextSelected(i)) {
204	>	if (sd->isDirectional()) {
205	>	Vector3d vec = sd->getA().getColumn(2);
206	>	vec.normalize();
207	>	angle += acos(dot(vec, director));
208	>	}
209	>	}
210	>	angle = angle/(seleMan1_.getSelectionCount()*NumericConstant::PI)*180.0;
211	>
212	>	} else {
213	>	for (vector<pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
214	>	Vector3d vec = j->first->getPos() - j->second->getPos();
215	>	if (usePeriodicBoundaryConditions_)
216	>	currentSnapshot_->wrapVector(vec);
217	>	vec.normalize();
218	>
219	>	angle += acos(dot(vec, director)) ;
220	>	}
221	>	angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
222		}
166	–	angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
223
224		OrderParam param;
225		param.p2 = p2;
#	Line 173 | Line 229 | namespace oopse {
229		orderParams_.push_back(param);
230
231		}
232	<
232	>
233		writeP2();
234	<
234	>
235		}
236
237		void P2OrderParameter::writeP2() {
238
239	<	std::ofstream os(getOutputFileName().c_str());
239	>	ofstream os(getOutputFileName().c_str());
240		os << "#radial distribution function\n";
241		os<< "#selection1: (" << selectionScript1_ << ")\t";
242	<	os << "selection2: (" << selectionScript2_ << ")\n";
242	>	if (!doVect_) {
243	>	os << "selection2: (" << selectionScript2_ << ")\n";
244	>	}
245		os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";
246
247	<	for (std::size_t i = 0; i < orderParams_.size(); ++i) {
247	>	for (size_t i = 0; i < orderParams_.size(); ++i) {
248		os <<  orderParams_[i].p2 << "\t"
249		<<  orderParams_[i].director[0] << "\t"
250		<<  orderParams_[i].director[1] << "\t"

Comparing trunk/src/applications/staticProps/P2OrderParameter.cpp (property svn:keywords):
Revision 1078 by gezelter, Wed Oct 18 21:58:48 2006 UTC vs.
Revision 1542 by gezelter, Thu Mar 3 20:32:49 2011 UTC

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

Diff Legend

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