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 1988 by gezelter, Fri Apr 18 19:49:54 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) {
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		void pAngle::process() {
116		Molecule* mol;
117		RigidBody* rb;
72	–	StuntDouble* sd;
118		SimInfo::MoleculeIterator mi;
119		Molecule::RigidBodyIterator rbIter;
120	<	int i;
120	>	StuntDouble* sd1;
121	>	StuntDouble* sd2;
122	>	int ii;
123	>	int jj;
124
125	+	Thermo thermo(info_);
126		DumpReader reader(info_, dumpFilename_);
127		int nFrames = reader.getNFrames();
128	+
129		nProcessed_ = nFrames/step_;
130
131		std::fill(histogram_.begin(), histogram_.end(), 0.0);
#	Line 94 | Line 144 | namespace OpenMD {
144		}
145		}
146
147	<	Vector3d CenterOfMass = info_->getCom();
147	>	Vector3d CenterOfMass = thermo.getCom();
148
149	<	if  (evaluator_.isDynamic()) {
150	<	seleMan_.setSelectionSet(evaluator_.evaluate());
149	>	if  (evaluator1_.isDynamic()) {
150	>	seleMan1_.setSelectionSet(evaluator1_.evaluate());
151		}
152
153	<	for (sd = seleMan_.beginSelected(i); sd != NULL;
154	<	sd = seleMan_.nextSelected(i)) {
153	>	if (doVect_) {
154	>
155
156	<	Vector3d pos = sd->getPos();
157	<
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();
156	>	for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
157	>	sd1 = seleMan1_.nextSelected(ii)) {
158
159	<	dipole.normalize();
160	<	r1.normalize();
161	<	RealType cosangle = dot(r1, dipole);
162	<
163	<	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
164	<	count_[binNo]++;
159	>	Vector3d pos = sd1->getPos();
160	>
161	>	Vector3d r1 = CenterOfMass - pos;
162	>	// only do this if the stunt double actually has a vector associated
163	>	// with it
164	>	if (sd1->isDirectional()) {
165	>	Vector3d vec = sd1->getA().getColumn(2);
166	>	RealType distance = r1.length();
167	>
168	>	vec.normalize();
169	>	r1.normalize();
170	>	RealType cosangle = dot(r1, vec);
171	>
172	>	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
173	>	count_[binNo]++;
174	>	}
175		}
176	<
176	>	} else {
177	>	if (doOffset_) {
178	>
179	>	for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
180	>	sd1 = seleMan1_.nextSelected(ii)) {
181	>
182	>	// This will require careful rewriting if StaticProps is
183	>	// ever parallelized.  For an example, see
184	>	// Thermo::getTaggedAtomPairDistance
185	>
186	>	int sd2Index = sd1->getGlobalIndex() + seleOffset_;
187	>	sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
188	>
189	>	Vector3d r1 = CenterOfMass - sd1->getPos();
190	>	if (usePeriodicBoundaryConditions_)
191	>	currentSnapshot_->wrapVector(r1);
192	>
193	>	Vector3d r2 = CenterOfMass - sd2->getPos();
194	>	if (usePeriodicBoundaryConditions_)
195	>	currentSnapshot_->wrapVector(r1);
196	>
197	>	Vector3d rc = 0.5*(r1 + r2);
198	>	if (usePeriodicBoundaryConditions_)
199	>	currentSnapshot_->wrapVector(rc);
200	>
201	>	Vector3d vec = r1-r2;
202	>	if (usePeriodicBoundaryConditions_)
203	>	currentSnapshot_->wrapVector(vec);
204	>
205	>	rc.normalize();
206	>	vec.normalize();
207	>	RealType cosangle = dot(rc, vec);
208	>	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
209	>	count_[binNo]++;
210	>	}
211	>	} else {
212	>
213	>	if  (evaluator2_.isDynamic()) {
214	>	seleMan2_.setSelectionSet(evaluator2_.evaluate());
215	>	}
216	>
217	>	if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
218	>	sprintf( painCave.errMsg,
219	>	"In frame %d, the number of selected StuntDoubles are\n"
220	>	"\tnot the same in --sele1 and sele2\n", istep);
221	>	painCave.severity = OPENMD_INFO;
222	>	painCave.isFatal = 0;
223	>	simError();
224	>	}
225	>
226	>	for (sd1 = seleMan1_.beginSelected(ii),
227	>	sd2 = seleMan2_.beginSelected(jj);
228	>	sd1 != NULL && sd2 != NULL;
229	>	sd1 = seleMan1_.nextSelected(ii),
230	>	sd2 = seleMan2_.nextSelected(jj)) {
231	>
232	>	Vector3d r1 = CenterOfMass - sd1->getPos();
233	>	if (usePeriodicBoundaryConditions_)
234	>	currentSnapshot_->wrapVector(r1);
235	>
236	>	Vector3d r2 = CenterOfMass - sd2->getPos();
237	>	if (usePeriodicBoundaryConditions_)
238	>	currentSnapshot_->wrapVector(r1);
239	>
240	>	Vector3d rc = 0.5*(r1 + r2);
241	>	if (usePeriodicBoundaryConditions_)
242	>	currentSnapshot_->wrapVector(rc);
243	>
244	>	Vector3d vec = r1-r2;
245	>	if (usePeriodicBoundaryConditions_)
246	>	currentSnapshot_->wrapVector(vec);
247	>
248	>	rc.normalize();
249	>	vec.normalize();
250	>	RealType cosangle = dot(rc, vec);
251	>	int binNo = int(nThetaBins_ * (1.0 + cosangle) / 2.0);
252	>	count_[binNo]++;
253	>
254	>	}
255	>	}
256		}
257		}
258	+
259		processHistogram();
260		writeProbs();
261
#	Line 130 | Line 264 | namespace OpenMD {
264		void pAngle::processHistogram() {
265
266		int atot = 0;
267	<	for(int i = 0; i < count_.size(); ++i)
267	>	for(unsigned int i = 0; i < count_.size(); ++i)
268		atot += count_[i];
269
270	<	for(int i = 0; i < count_.size(); ++i) {
271	<	histogram_[i] = double(count_[i] / atot);
270	>	for(unsigned int i = 0; i < count_.size(); ++i) {
271	>	histogram_[i] = double(count_[i] / double(atot));
272		}
273		}
274
#	Line 145 | Line 279 | namespace OpenMD {
279		if (rdfStream.is_open()) {
280		rdfStream << "#pAngle\n";
281		rdfStream << "#nFrames:\t" << nProcessed_ << "\n";
282	<	rdfStream << "#selection: (" << selectionScript_ << ")\n";
282	>	rdfStream << "#selection1: (" << selectionScript1_ << ")";
283	>	if (!doVect_) {
284	>	rdfStream << "\tselection2: (" << selectionScript2_ << ")";
285	>	}
286	>	rdfStream << "\n";
287		rdfStream << "#cos(theta)\tp(cos(theta))\n";
288	<	for (int i = 0; i < histogram_.size(); ++i) {
289	<	RealType ct = -1.0 + i / histogram_.size();
290	<	rdfStream << ct << "\t" << histogram_[i] << "\n";
288	>	RealType dct = 2.0 / histogram_.size();
289	>	for (unsigned int i = 0; i < histogram_.size(); ++i) {
290	>	RealType ct = -1.0 + (2.0 * i + 1) / (histogram_.size());
291	>	rdfStream << ct << "\t" << histogram_[i]/dct << "\n";
292		}
293
294		} else {
295
296	<	sprintf(painCave.errMsg, "pAngle: unable to open %s\n", outputFilename_.c_str());
296	>	sprintf(painCave.errMsg, "pAngle: unable to open %s\n",
297	>	outputFilename_.c_str());
298		painCave.isFatal = 1;
299		simError();
300		}

Comparing trunk/src/applications/staticProps/pAngle.cpp (property svn:keywords):
Revision 1413 by gezelter, Mon Mar 22 19:21:22 2010 UTC vs.
Revision 1988 by gezelter, Fri Apr 18 19:49:54 2014 UTC

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

Diff Legend

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