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 353 by tim, Wed Feb 16 19:36:30 2005 UTC vs.
Revision 2023 by gezelter, Thu Oct 2 14:35:14 2014 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, 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 oopse {
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, const std::string& sele2)
80	<	: RadialDistrFunc(info, filename, sele1, sele2){
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	<	}
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	<	void GofRAngle::preProcess() {
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 (int i = 0; i < avgGofr_.size(); ++i) {
118	<	std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
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	<	}
128	>	}
129
130	<	void GofRAngle::initalizeHistogram() {
131	<	npairs_ = 0;
132	<	for (int i = 0; i < histogram_.size(); ++i)
133	<	std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
134	<	}
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::processHistogram() {
164	>	void GofRAngle::preProcess() {
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::initializeHistogram() {
171	+	npairs_ = 0;
172	+	for (unsigned int i = 0; i < histogram_.size(); ++i){
173	+	std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
174	+	}
175	+	}
176	+
177	+	void GofRAngle::processHistogram() {
178		int nPairs = getNPairs();
179	<	double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
180	<	double pairDensity = nPairs /volume;
181	<	double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
179	>	RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
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	<	double rLower = i * deltaR_;
186	<	double rUpper = rLower + deltaR_;
187	<	double volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
188	<	double nIdeal = volSlice * pairConstant;
185	>	RealType rLower = i * deltaR_;
186	>	RealType rUpper = rLower + deltaR_;
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){
192	<	avgGofr_[i][j] += histogram_[i][j] / nIdeal;
193	<	}
191	>	for (unsigned int j = 0; j < histogram_[i].size(); ++j){
192	>	avgGofr_[i][j] += histogram_[i][j] / nIdeal;
193	>	}
194		}
195
196	<	}
196	>	}
197
198	<	void GofRAngle::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
198	>	void GofRAngle::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
199
200		if (sd1 == sd2) {
201	<	return;
201	>	return;
202		}
95	–
203		Vector3d pos1 = sd1->getPos();
204		Vector3d pos2 = sd2->getPos();
205	<	Vector3d r12 = pos1 - pos2;
206	<	currentSnapshot_->wrapVector(r12);
205	>	Vector3d r12 = pos2 - pos1;
206	>	if (usePeriodicBoundaryConditions_)
207	>	currentSnapshot_->wrapVector(r12);
208
209	<	double distance = r12.length();
210	<	int whichRBin = distance / deltaR_;
209	>	RealType distance = r12.length();
210	>	int whichRBin = int(distance / deltaR_);
211
212		if (distance <= len_) {
213	<	double cosAngle = evaluateAngle(sd1, sd2);
214	<	double halfBin = (nAngleBins_ - 1) * 0.5;
215	<	int whichThetaBin = halfBin * (cosAngle + 1.0);
216	<	++histogram_[whichRBin][whichThetaBin];
213	>
214	>	RealType cosAngle = evaluateAngle(sd1, sd2);
215	>	RealType halfBin = (nAngleBins_ - 1) * 0.5;
216	>	int whichThetaBin = int(halfBin * (cosAngle + 1.0));
217	>	++histogram_[whichRBin][whichThetaBin];
218
219	<	++npairs_;
219	>	++npairs_;
220		}
221	<	}
221	>	}
222
223	<	void GofRAngle::writeRdf() {
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 << "#r\tcorrValue\n";
264	<	for (int i = 0; i < avgGofr_.size(); ++i) {
265	<	double r = deltaR_ * (i + 0.5);
260	>	rdfStream << "#radial distribution function\n";
261	>	rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
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	<	double cosAngle = -1.0 + (i + 0.5)*deltaCosAngle_;
277	<	rdfStream << r << "\t" << cosAngle << "\t" << avgGofr_[i][j]/nProcessed_ << "\n";
278	<	}
279	<	}
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	>
280	>	rdfStream << "\n";
281	>	}
282
283		} else {
284	<
285	<
284	>	sprintf(painCave.errMsg, "GofRAngle: unable to open %s\n",
285	>	outputFilename_.c_str());
286	>	painCave.isFatal = 1;
287	>	simError();
288		}
289
290		rdfStream.close();
291	<	}
291	>	}
292
293	<	double GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
293	>	RealType GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
294		Vector3d pos1 = sd1->getPos();
295		Vector3d pos2 = sd2->getPos();
296	<	Vector3d r12 = pos1 - pos2;
297	<	currentSnapshot_->wrapVector(r12);
296	>	Vector3d r12 = pos2 - pos1;
297	>
298	>	if (usePeriodicBoundaryConditions_)
299	>	currentSnapshot_->wrapVector(r12);
300	>
301		r12.normalize();
144	–	Vector3d dipole = sd1->getElectroFrame().getColumn(2);
145	–	dipole.normalize();
146	–	return dot(r12, dipole);
147	–	}
302
303	<	double GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
304	<	Vector3d v1 = sd1->getElectroFrame().getColumn(2);
305	<	Vector3d v2 = sd1->getElectroFrame().getColumn(2);
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, 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	<	}
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 353 by tim, Wed Feb 16 19:36:30 2005 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