ViewVC Help

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

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

Comparing trunk/src/applications/staticProps/pAngle.cpp (file contents):
Revision 1413 by gezelter, Mon Mar 22 19:21:22 2010 UTC vs.
Revision 1991 by gezelter, Wed Apr 23 20:34:17 2014 UTC

#	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).
40	<	*
38	>	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	>	* [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	>	* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41		*/
42
43		/* Calculates Rho(theta) */
44
45		#include <algorithm>
46	<	#include <fstream>
46	>	#include <fstream>
47		#include "applications/staticProps/pAngle.hpp"
48		#include "utils/simError.h"
49		#include "io/DumpReader.hpp"
50		#include "primitives/Molecule.hpp"
51	+	#include "brains/Thermo.hpp"
52	+
53		namespace OpenMD {
54
55		pAngle::pAngle(SimInfo* info, const std::string& filename,
56	<	const std::string& sele, int nthetabins)
57	<	: StaticAnalyser(info, filename), selectionScript_(sele),
58	<	evaluator_(info), seleMan_(info), nThetaBins_(nthetabins){
56	>	const std::string& sele1, int nthetabins)
57	>	: StaticAnalyser(info, filename), selectionScript1_(sele1),
58	>	evaluator1_(info),  evaluator2_(info), seleMan1_(info), seleMan2_(info),
59	>	nThetaBins_(nthetabins),
60	>	doVect_(true), doOffset_(false) {
61
62	<	evaluator_.loadScriptString(sele);
63	<	if (!evaluator_.isDynamic()) {
64	<	seleMan_.setSelectionSet(evaluator_.evaluate());
62	>	setOutputName(getPrefix(filename) + ".pAngle");
63	>
64	>	evaluator1_.loadScriptString(sele1);
65	>	if (!evaluator1_.isDynamic()) {
66	>	seleMan1_.setSelectionSet(evaluator1_.evaluate());
67		}
68
69		count_.resize(nThetaBins_);
70	<	histogram_.resize(nThetaBins_);
70	>	histogram_.resize(nThetaBins_);
71	>	}
72	>
73	>	pAngle::pAngle(SimInfo* info, const std::string& filename,
74	>	const std::string& sele1, const std::string& sele2,
75	>	int nthetabins)
76	>	: StaticAnalyser(info, filename), selectionScript1_(sele1),
77	>	selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
78	>	seleMan1_(info), seleMan2_(info), nThetaBins_(nthetabins),
79	>	doVect_(false), doOffset_(false) {
80
81		setOutputName(getPrefix(filename) + ".pAngle");
82	+
83	+	evaluator1_.loadScriptString(sele1);
84	+	if (!evaluator1_.isDynamic()) {
85	+	seleMan1_.setSelectionSet(evaluator1_.evaluate());
86	+	}
87	+
88	+	evaluator2_.loadScriptString(sele2);
89	+	if (!evaluator2_.isDynamic()) {
90	+	seleMan2_.setSelectionSet(evaluator2_.evaluate());
91	+	}
92	+
93	+	count_.resize(nThetaBins_);
94	+	histogram_.resize(nThetaBins_);
95		}
96	+
97	+	pAngle::pAngle(SimInfo* info, const std::string& filename,
98	+	const std::string& sele1, int seleOffset, int nthetabins)
99	+	: StaticAnalyser(info, filename), selectionScript1_(sele1),
100	+	evaluator1_(info), evaluator2_(info), seleMan1_(info), seleMan2_(info),
101	+	nThetaBins_(nthetabins), seleOffset_(seleOffset),
102	+	doVect_(false), doOffset_(true), doOffset2_(false) {
103	+
104	+	setOutputName(getPrefix(filename) + ".pAngle");
105	+
106	+	evaluator1_.loadScriptString(sele1);
107	+	if (!evaluator1_.isDynamic()) {
108	+	seleMan1_.setSelectionSet(evaluator1_.evaluate());
109	+	}
110	+
111	+	count_.resize(nThetaBins_);
112	+	histogram_.resize(nThetaBins_);
113	+	}
114	+
115	+	pAngle::pAngle(SimInfo* info, const std::string& filename,
116	+	const std::string& sele1, int seleOffset, int seleOffset2,
117	+	int nthetabins)
118	+	: StaticAnalyser(info, filename), selectionScript1_(sele1),
119	+	evaluator1_(info), evaluator2_(info), seleMan1_(info), seleMan2_(info),
120	+	nThetaBins_(nthetabins), seleOffset_(seleOffset),
121	+	seleOffset2_(seleOffset2),
122	+	doVect_(false), doOffset_(true), doOffset2_(true) {
123	+
124	+	setOutputName(getPrefix(filename) + ".pAngle");
125	+
126	+	evaluator1_.loadScriptString(sele1);
127	+	if (!evaluator1_.isDynamic()) {
128	+	seleMan1_.setSelectionSet(evaluator1_.evaluate());
129	+	}
130	+
131	+	count_.resize(nThetaBins_);
132	+	histogram_.resize(nThetaBins_);
133	+	}
134
135		void pAngle::process() {
136		Molecule* mol;
137		RigidBody* rb;
72	–	StuntDouble* sd;
138		SimInfo::MoleculeIterator mi;
139		Molecule::RigidBodyIterator rbIter;
140	<	int i;
140	>	StuntDouble* sd1;
141	>	StuntDouble* sd2;
142	>	int ii;
143	>	int jj;
144
145	+	Thermo thermo(info_);
146		DumpReader reader(info_, dumpFilename_);
147		int nFrames = reader.getNFrames();
148	+
149		nProcessed_ = nFrames/step_;
150
151		std::fill(histogram_.begin(), histogram_.end(), 0.0);
#	Line 94 | Line 164 | namespace OpenMD {
164		}
165		}
166
167	<	Vector3d CenterOfMass = info_->getCom();
167	>	Vector3d CenterOfMass = thermo.getCom();
168
169	<	if  (evaluator_.isDynamic()) {
170	<	seleMan_.setSelectionSet(evaluator_.evaluate());
169	>	if  (evaluator1_.isDynamic()) {
170	>	seleMan1_.setSelectionSet(evaluator1_.evaluate());
171		}
172
173	<	for (sd = seleMan_.beginSelected(i); sd != NULL;
174	<	sd = seleMan_.nextSelected(i)) {
173	>	if (doVect_) {
174	>
175
176	<	Vector3d pos = sd->getPos();
177	<
108	<	Vector3d r1 = CenterOfMass - pos;
109	<	// only do this if the stunt double actually has a vector associated
110	<	// with it
111	<	if (sd->isDirectional()) {
112	<	Vector3d dipole = sd->getA().getColumn(2);
113	<	RealType distance = r1.length();
176	>	for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
177	>	sd1 = seleMan1_.nextSelected(ii)) {
178
179	<	dipole.normalize();
180	<	r1.normalize();
181	<	RealType cosangle = dot(r1, dipole);
182	<
183	<	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
184	<	count_[binNo]++;
179	>	Vector3d pos = sd1->getPos();
180	>
181	>	Vector3d r1 = CenterOfMass - pos;
182	>	// only do this if the stunt double actually has a vector associated
183	>	// with it
184	>	if (sd1->isDirectional()) {
185	>	Vector3d vec = sd1->getA().getColumn(2);
186	>	RealType distance = r1.length();
187	>
188	>	vec.normalize();
189	>	r1.normalize();
190	>	RealType cosangle = dot(r1, vec);
191	>
192	>	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
193	>	count_[binNo]++;
194	>	}
195		}
196	<
196	>	} else {
197	>	if (doOffset_) {
198	>
199	>	for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
200	>	sd1 = seleMan1_.nextSelected(ii)) {
201	>
202	>	// This will require careful rewriting if StaticProps is
203	>	// ever parallelized.  For an example, see
204	>	// Thermo::getTaggedAtomPairDistance
205	>	Vector3d r1;
206	>
207	>	if (doOffset2_) {
208	>	int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
209	>	StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
210	>	r1 = CenterOfMass - sd1A->getPos();
211	>	} else {
212	>	r1 = CenterOfMass - sd1->getPos();
213	>	}
214	>
215	>	if (usePeriodicBoundaryConditions_)
216	>	currentSnapshot_->wrapVector(r1);
217	>
218	>
219	>	int sd2Index = sd1->getGlobalIndex() + seleOffset_;
220	>	sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
221	>
222	>	Vector3d r2 = CenterOfMass - sd2->getPos();
223	>	if (usePeriodicBoundaryConditions_)
224	>	currentSnapshot_->wrapVector(r1);
225	>
226	>	Vector3d rc = 0.5*(r1 + r2);
227	>	if (usePeriodicBoundaryConditions_)
228	>	currentSnapshot_->wrapVector(rc);
229	>
230	>	Vector3d vec = r1-r2;
231	>	if (usePeriodicBoundaryConditions_)
232	>	currentSnapshot_->wrapVector(vec);
233	>
234	>	rc.normalize();
235	>	vec.normalize();
236	>	RealType cosangle = dot(rc, vec);
237	>	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
238	>	count_[binNo]++;
239	>	}
240	>	} else {
241	>
242	>	if  (evaluator2_.isDynamic()) {
243	>	seleMan2_.setSelectionSet(evaluator2_.evaluate());
244	>	}
245	>
246	>	if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
247	>	sprintf( painCave.errMsg,
248	>	"In frame %d, the number of selected StuntDoubles are\n"
249	>	"\tnot the same in --sele1 and sele2\n", istep);
250	>	painCave.severity = OPENMD_INFO;
251	>	painCave.isFatal = 0;
252	>	simError();
253	>	}
254	>
255	>	for (sd1 = seleMan1_.beginSelected(ii),
256	>	sd2 = seleMan2_.beginSelected(jj);
257	>	sd1 != NULL && sd2 != NULL;
258	>	sd1 = seleMan1_.nextSelected(ii),
259	>	sd2 = seleMan2_.nextSelected(jj)) {
260	>
261	>	Vector3d r1 = CenterOfMass - sd1->getPos();
262	>	if (usePeriodicBoundaryConditions_)
263	>	currentSnapshot_->wrapVector(r1);
264	>
265	>	Vector3d r2 = CenterOfMass - sd2->getPos();
266	>	if (usePeriodicBoundaryConditions_)
267	>	currentSnapshot_->wrapVector(r1);
268	>
269	>	Vector3d rc = 0.5*(r1 + r2);
270	>	if (usePeriodicBoundaryConditions_)
271	>	currentSnapshot_->wrapVector(rc);
272	>
273	>	Vector3d vec = r1-r2;
274	>	if (usePeriodicBoundaryConditions_)
275	>	currentSnapshot_->wrapVector(vec);
276	>
277	>	rc.normalize();
278	>	vec.normalize();
279	>	RealType cosangle = dot(rc, vec);
280	>	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
281	>	count_[binNo]++;
282	>
283	>	}
284	>	}
285		}
286		}
287	+
288		processHistogram();
289		writeProbs();
290
#	Line 130 | Line 293 | namespace OpenMD {
293		void pAngle::processHistogram() {
294
295		int atot = 0;
296	<	for(int i = 0; i < count_.size(); ++i)
296	>	for(unsigned int i = 0; i < count_.size(); ++i)
297		atot += count_[i];
298
299	<	for(int i = 0; i < count_.size(); ++i) {
300	<	histogram_[i] = double(count_[i] / atot);
299	>	for(unsigned int i = 0; i < count_.size(); ++i) {
300	>	histogram_[i] = double(count_[i] / double(atot));
301		}
302		}
303
#	Line 145 | Line 308 | namespace OpenMD {
308		if (rdfStream.is_open()) {
309		rdfStream << "#pAngle\n";
310		rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
311	<	rdfStream << "#selection: (" << selectionScript_ << ")\n";
311	>	rdfStream << "#selection1: (" << selectionScript1_ << ")";
312	>	if (!doVect_) {
313	>	rdfStream << "\tselection2: (" << selectionScript2_ << ")";
314	>	}
315	>	rdfStream << "\n";
316		rdfStream << "#cos(theta)\tp(cos(theta))\n";
317	<	for (int i = 0; i < histogram_.size(); ++i) {
318	<	RealType ct = -1.0 + i / histogram_.size();
319	<	rdfStream << ct << "\t" << histogram_[i] << "\n";
317	>	RealType dct = 2.0 / histogram_.size();
318	>	for (unsigned int i = 0; i < histogram_.size(); ++i) {
319	>	RealType ct = -1.0 + (2.0 * i + 1) / (histogram_.size());
320	>	rdfStream << ct << "\t" << histogram_[i]/dct << "\n";
321		}
322
323		} else {
324
325	<	sprintf(painCave.errMsg, "pAngle: unable to open %s\n", outputFilename_.c_str());
325	>	sprintf(painCave.errMsg, "pAngle: unable to open %s\n",
326	>	outputFilename_.c_str());
327		painCave.isFatal = 1;
328		simError();
329		}

Comparing trunk/src/applications/staticProps/pAngle.cpp (property svn:keywords):
Revision 1413 by gezelter, Mon Mar 22 19:21:22 2010 UTC vs.
Revision 1991 by gezelter, Wed Apr 23 20:34:17 2014 UTC

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

Diff Legend

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