ViewVC Help

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

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

Comparing trunk/src/applications/staticProps/GofRAngle.cpp (file contents):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 2023 by gezelter, Thu Oct 2 14:35:14 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).
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 <algorithm>
44		#include <fstream>
45		#include "applications/staticProps/GofRAngle.hpp"
46	+	#include "primitives/Atom.hpp"
47	+	#include "types/MultipoleAdapter.hpp"
48		#include "utils/simError.h"
49
50		namespace OpenMD {
51	+
52	+	GofRAngle::GofRAngle(SimInfo* info, const std::string& filename,
53	+	const std::string& sele1,
54	+	const std::string& sele2,
55	+	RealType len, int nrbins, int nangleBins)
56	+	: RadialDistrFunc(info, filename, sele1, sele2), len_(len),
57	+	nRBins_(nrbins), nAngleBins_(nangleBins), evaluator3_(info),
58	+	seleMan3_(info), doSele3_(false) {
59	+
60	+	deltaR_ = len_ /(double) nRBins_;
61	+	deltaCosAngle_ = 2.0 / (double)nAngleBins_;
62	+	histogram_.resize(nRBins_);
63	+	avgGofr_.resize(nRBins_);
64	+	for (int i = 0 ; i < nRBins_; ++i) {
65	+	histogram_[i].resize(nAngleBins_);
66	+	avgGofr_[i].resize(nAngleBins_);
67	+	}
68	+	}
69	+
70	+	GofRAngle::GofRAngle(SimInfo* info, const std::string& filename,
71	+	const std::string& sele1,
72	+	const std::string& sele2,
73	+	const std::string& sele3,
74	+	RealType len, int nrbins, int nangleBins)
75	+	: RadialDistrFunc(info, filename, sele1, sele2), len_(len),
76	+	nRBins_(nrbins), nAngleBins_(nangleBins), selectionScript3_(sele3),
77	+	evaluator3_(info), seleMan3_(info), doSele3_(true) {
78
79	<	GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1,
80	<	const std::string& sele2, RealType len, int nrbins, int nangleBins)
81	<	: RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins), nAngleBins_(nangleBins){
79	>	deltaR_ = len_ /(double) nRBins_;
80	>	deltaCosAngle_ = 2.0 / (double)nAngleBins_;
81	>	histogram_.resize(nRBins_);
82	>	avgGofr_.resize(nRBins_);
83	>	for (int i = 0 ; i < nRBins_; ++i) {
84	>	histogram_[i].resize(nAngleBins_);
85	>	avgGofr_[i].resize(nAngleBins_);
86	>	}
87
88	<	deltaR_ = len_ /(double) nRBins_;
89	<	deltaCosAngle_ = 2.0 / (double)nAngleBins_;
90	<	histogram_.resize(nRBins_);
91	<	avgGofr_.resize(nRBins_);
92	<	for (int i = 0 ; i < nRBins_; ++i) {
93	<	histogram_[i].resize(nAngleBins_);
94	<	avgGofr_[i].resize(nAngleBins_);
88	>	evaluator3_.loadScriptString(sele3);
89	>	if (!evaluator3_.isDynamic()) {
90	>	seleMan3_.setSelectionSet(evaluator3_.evaluate());
91	>	}
92	>
93	>	}
94	>
95	>	void GofRAngle::processNonOverlapping( SelectionManager& sman1,
96	>	SelectionManager& sman2) {
97	>	StuntDouble* sd1;
98	>	StuntDouble* sd2;
99	>	StuntDouble* sd3;
100	>	int i;
101	>	int j;
102	>	int k;
103	>
104	>	// This is the same as a non-overlapping pairwise loop structure:
105	>	// for (int i = 0;  i < ni ; ++i ) {
106	>	//   for (int j = 0; j < nj; ++j) {}
107	>	// }
108	>
109	>	if (doSele3_) {
110	>	if  (evaluator3_.isDynamic()) {
111	>	seleMan3_.setSelectionSet(evaluator3_.evaluate());
112		}
113	+	if (sman1.getSelectionCount() != seleMan3_.getSelectionCount() ) {
114	+	RadialDistrFunc::processNonOverlapping( sman1, sman2 );
115	+	}
116	+
117	+	for (sd1 = sman1.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
118	+	sd1 != NULL && sd3 != NULL;
119	+	sd1 = sman1.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
120	+	for (sd2 = sman2.beginSelected(j); sd2 != NULL;
121	+	sd2 = sman2.nextSelected(j)) {
122	+	collectHistogram(sd1, sd2, sd3);
123	+	}
124	+	}
125	+	} else {
126	+	RadialDistrFunc::processNonOverlapping( sman1, sman2 );
127		}
128	+	}
129
130	+	void GofRAngle::processOverlapping( SelectionManager& sman) {
131	+	StuntDouble* sd1;
132	+	StuntDouble* sd2;
133	+	StuntDouble* sd3;
134	+	int i;
135	+	int j;
136	+	int k;
137
138	+	// This is the same as a pairwise loop structure:
139	+	// for (int i = 0;  i < n-1 ; ++i ) {
140	+	//   for (int j = i + 1; j < n; ++j) {}
141	+	// }
142	+
143	+	if (doSele3_) {
144	+	if  (evaluator3_.isDynamic()) {
145	+	seleMan3_.setSelectionSet(evaluator3_.evaluate());
146	+	}
147	+	if (sman.getSelectionCount() != seleMan3_.getSelectionCount() ) {
148	+	RadialDistrFunc::processOverlapping( sman);
149	+	}
150	+	for (sd1 = sman.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
151	+	sd1 != NULL && sd3 != NULL;
152	+	sd1 = sman.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
153	+	for (j  = i, sd2 = sman.nextSelected(j); sd2 != NULL;
154	+	sd2 = sman.nextSelected(j)) {
155	+	collectHistogram(sd1, sd2, sd3);
156	+	}
157	+	}
158	+	} else {
159	+	RadialDistrFunc::processOverlapping( sman);
160	+	}
161	+	}
162	+
163	+
164		void GofRAngle::preProcess() {
165	<	for (int i = 0; i < avgGofr_.size(); ++i) {
165	>	for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
166		std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
167		}
168		}
169
170	<	void GofRAngle::initalizeHistogram() {
170	>	void GofRAngle::initializeHistogram() {
171		npairs_ = 0;
172	<	for (int i = 0; i < histogram_.size(); ++i){
172	>	for (unsigned int i = 0; i < histogram_.size(); ++i){
173		std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
174		}
175		}
#	Line 80 | Line 180 | namespace OpenMD {
180		RealType pairDensity = nPairs /volume;
181		RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
182
183	<	for(int i = 0 ; i < histogram_.size(); ++i){
183	>	for(unsigned int i = 0 ; i < histogram_.size(); ++i){
184
185		RealType rLower = i * deltaR_;
186		RealType rUpper = rLower + deltaR_;
187	<	RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
187	>	RealType volSlice = ( rUpper * rUpper * rUpper ) -
188	>	( rLower * rLower * rLower );
189		RealType nIdeal = volSlice * pairConstant;
190
191	<	for (int j = 0; j < histogram_[i].size(); ++j){
191	>	for (unsigned int j = 0; j < histogram_[i].size(); ++j){
192		avgGofr_[i][j] += histogram_[i][j] / nIdeal;
193		}
194		}
#	Line 106 | Line 207 | namespace OpenMD {
207		currentSnapshot_->wrapVector(r12);
208
209		RealType distance = r12.length();
210	<	int whichRBin = distance / deltaR_;
210	>	int whichRBin = int(distance / deltaR_);
211
212		if (distance <= len_) {
213
214		RealType cosAngle = evaluateAngle(sd1, sd2);
215		RealType halfBin = (nAngleBins_ - 1) * 0.5;
216	<	int whichThetaBin = halfBin * (cosAngle + 1.0);
216	>	int whichThetaBin = int(halfBin * (cosAngle + 1.0));
217		++histogram_[whichRBin][whichThetaBin];
218
219		++npairs_;
220		}
221		}
222
223	+	void GofRAngle::collectHistogram(StuntDouble* sd1, StuntDouble* sd2,
224	+	StuntDouble* sd3) {
225	+
226	+	if (sd1 == sd2) {
227	+	return;
228	+	}
229	+
230	+	Vector3d p1 = sd1->getPos();
231	+	Vector3d p3 = sd3->getPos();
232	+
233	+	Vector3d c = 0.5 * (p1 + p3);
234	+	Vector3d r13 = p3 - p1;
235	+
236	+	Vector3d r12 = sd2->getPos() - c;
237	+
238	+	if (usePeriodicBoundaryConditions_) {
239	+	currentSnapshot_->wrapVector(r12);
240	+	currentSnapshot_->wrapVector(r13);
241	+	}
242	+
243	+	RealType distance = r12.length();
244	+	int whichRBin = int(distance / deltaR_);
245	+
246	+	if (distance <= len_) {
247	+
248	+	RealType cosAngle = evaluateAngle(sd1, sd2, sd3);
249	+	RealType halfBin = (nAngleBins_ - 1) * 0.5;
250	+	int whichThetaBin = int(halfBin * (cosAngle + 1.0));
251	+	++histogram_[whichRBin][whichThetaBin];
252	+
253	+	++npairs_;
254	+	}
255	+	}
256	+
257		void GofRAngle::writeRdf() {
258		std::ofstream rdfStream(outputFilename_.c_str());
259		if (rdfStream.is_open()) {
260		rdfStream << "#radial distribution function\n";
261		rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
262	<	rdfStream << "selection2: (" << selectionScript2_ << ")\n";
263	<	rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
264	<	rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = " << deltaCosAngle_ << "\n";
265	<	for (int i = 0; i < avgGofr_.size(); ++i) {
266	<	RealType r = deltaR_ * (i + 0.5);
262	>	rdfStream <<  "selection2: (" << selectionScript2_ << ")";
263	>	if (doSele3_) {
264	>	rdfStream << "\tselection3: (" << selectionScript3_ << ")\n";
265	>	} else {
266	>	rdfStream << "\n";
267	>	}
268	>	rdfStream << "#nRBins = " << nRBins_ << "\tmaxLen = "
269	>	<< len_ << "\tdeltaR = " << deltaR_ <<"\n";
270	>	rdfStream << "#nAngleBins =" << nAngleBins_ << "\tdeltaCosAngle = "
271	>	<< deltaCosAngle_ << "\n";
272	>	for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
273	>	// RealType r = deltaR_ * (i + 0.5);
274
275	<	for(int j = 0; j < avgGofr_[i].size(); ++j) {
276	<	RealType cosAngle = -1.0 + (j + 0.5)*deltaCosAngle_;
275	>	for(unsigned int j = 0; j < avgGofr_[i].size(); ++j) {
276	>	// RealType cosAngle = -1.0 + (j + 0.5)*deltaCosAngle_;
277		rdfStream << avgGofr_[i][j]/nProcessed_ << "\t";
278		}
279
#	Line 139 | Line 281 | namespace OpenMD {
281		}
282
283		} else {
284	<	sprintf(painCave.errMsg, "GofRAngle: unable to open %s\n", outputFilename_.c_str());
284	>	sprintf(painCave.errMsg, "GofRAngle: unable to open %s\n",
285	>	outputFilename_.c_str());
286		painCave.isFatal = 1;
287		simError();
288		}
#	Line 156 | Line 299 | namespace OpenMD {
299		currentSnapshot_->wrapVector(r12);
300
301		r12.normalize();
302	<	Vector3d dipole = sd1->getElectroFrame().getColumn(2);
303	<	dipole.normalize();
304	<	return dot(r12, dipole);
302	>
303	>	Vector3d vec;
304	>
305	>	if (!sd1->isDirectional()) {
306	>	sprintf(painCave.errMsg,
307	>	"GofRTheta: attempted to use a non-directional object: %s\n",
308	>	sd1->getType().c_str());
309	>	painCave.isFatal = 1;
310	>	simError();
311	>	}
312	>
313	>	if (sd1->isAtom()) {
314	>	AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
315	>	MultipoleAdapter ma1 = MultipoleAdapter(atype1);
316	>
317	>	if (ma1.isDipole() )
318	>	vec = sd1->getDipole();
319	>	else
320	>	vec = sd1->getA().transpose() * V3Z;
321	>	} else {
322	>	vec = sd1->getA().transpose() * V3Z;
323	>	}
324	>
325	>	vec.normalize();
326	>
327	>	return dot(r12, vec);
328		}
329
330	+	RealType GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2,
331	+	StuntDouble* sd3) {
332	+	Vector3d p1 = sd1->getPos();
333	+	Vector3d p3 = sd3->getPos();
334	+
335	+	Vector3d c = 0.5 * (p1 + p3);
336	+	Vector3d r13 = p3 - p1;
337	+
338	+	Vector3d r12 = sd2->getPos() - c;
339	+
340	+	if (usePeriodicBoundaryConditions_) {
341	+	currentSnapshot_->wrapVector(r12);
342	+	currentSnapshot_->wrapVector(r13);
343	+	}
344	+
345	+	r12.normalize();
346	+	r13.normalize();
347	+
348	+	return dot(r12, r13);
349	+	}
350	+
351		RealType GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
352	<	Vector3d v1 = sd1->getElectroFrame().getColumn(2);
353	<	Vector3d v2 = sd2->getElectroFrame().getColumn(2);
352	>	Vector3d v1, v2;
353	>
354	>	if (!sd1->isDirectional()) {
355	>	sprintf(painCave.errMsg,
356	>	"GofROmega: attempted to use a non-directional object: %s\n",
357	>	sd1->getType().c_str());
358	>	painCave.isFatal = 1;
359	>	simError();
360	>	}
361	>
362	>	if (sd1->isAtom()){
363	>	AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
364	>	MultipoleAdapter ma1 = MultipoleAdapter(atype1);
365	>	if (ma1.isDipole() )
366	>	v1 = sd1->getDipole();
367	>	else
368	>	v1 = sd1->getA().transpose() * V3Z;
369	>	} else {
370	>	v1 = sd1->getA().transpose() * V3Z;
371	>	}
372	>
373	>	if (!sd2->isDirectional()) {
374	>	sprintf(painCave.errMsg,
375	>	"GofROmega attempted to use a non-directional object: %s\n",
376	>	sd2->getType().c_str());
377	>	painCave.isFatal = 1;
378	>	simError();
379	>	}
380	>
381	>	if (sd2->isAtom()) {
382	>	AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
383	>	MultipoleAdapter ma2 = MultipoleAdapter(atype2);
384	>
385	>	if (ma2.isDipole() )
386	>	v2 = sd2->getDipole();
387	>	else
388	>	v2 = sd2->getA().transpose() * V3Z;
389	>	} else {
390	>	v2 = sd2->getA().transpose() * V3Z;
391	>	}
392	>
393		v1.normalize();
394		v2.normalize();
395		return dot(v1, v2);
396		}
397
398	+	RealType GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2,
399	+	StuntDouble* sd3) {
400
401	+	Vector3d v1;
402	+	Vector3d v2;
403	+
404	+	v1 = sd3->getPos() - sd1->getPos();
405	+	if (usePeriodicBoundaryConditions_)
406	+	currentSnapshot_->wrapVector(v1);
407	+
408	+	if (!sd2->isDirectional()) {
409	+	sprintf(painCave.errMsg,
410	+	"GofROmega: attempted to use a non-directional object: %s\n",
411	+	sd2->getType().c_str());
412	+	painCave.isFatal = 1;
413	+	simError();
414	+	}
415	+
416	+	if (sd2->isAtom()) {
417	+	AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
418	+	MultipoleAdapter ma2 = MultipoleAdapter(atype2);
419	+
420	+	if (ma2.isDipole() )
421	+	v2 = sd2->getDipole();
422	+	else
423	+	v2 = sd2->getA().transpose() * V3Z;
424	+	} else {
425	+	v2 = sd2->getA().transpose() * V3Z;
426	+	}
427	+
428	+	v1.normalize();
429	+	v2.normalize();
430	+	return dot(v1, v2);
431	+	}
432		}
433
434

Comparing trunk/src/applications/staticProps/GofRAngle.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 2023 by gezelter, Thu Oct 2 14:35:14 2014 UTC

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

Diff Legend

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