ViewVC Help

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

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

Comparing trunk/src/applications/staticProps/GofAngle2.cpp (file contents):
Revision 353 by tim, Wed Feb 16 19:36:30 2005 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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/GofAngle2.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	>	GofAngle2::GofAngle2(SimInfo* info, const std::string& filename,
53	>	const std::string& sele1,
54	>	const std::string& sele2, int nangleBins)
55	>	: RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins),
56	>	doSele3_(false), seleMan3_(info), evaluator3_(info) {
57	>
58	>	setOutputName(getPrefix(filename) + ".gto");
59	>
60	>	deltaCosAngle_ = 2.0 / nAngleBins_;
61	>
62	>	histogram_.resize(nAngleBins_);
63	>	avgGofr_.resize(nAngleBins_);
64	>	for (int i = 0 ; i < nAngleBins_; ++i) {
65	>	histogram_[i].resize(nAngleBins_);
66	>	avgGofr_[i].resize(nAngleBins_);
67	>	}
68	>	}
69	>
70	>	GofAngle2::GofAngle2(SimInfo* info, const std::string& filename,
71	>	const std::string& sele1,
72	>	const std::string& sele2,
73	>	const std::string& sele3, int nangleBins)
74	>	: RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins),
75	>	doSele3_(true), seleMan3_(info), evaluator3_(info),
76	>	selectionScript3_(sele3) {
77	>
78	>	setOutputName(getPrefix(filename) + ".gto");
79	>
80	>	deltaCosAngle_ = 2.0 / nAngleBins_;
81	>
82	>	histogram_.resize(nAngleBins_);
83	>	avgGofr_.resize(nAngleBins_);
84	>	for (int i = 0 ; i < nAngleBins_; ++i) {
85	>	histogram_[i].resize(nAngleBins_);
86	>	avgGofr_[i].resize(nAngleBins_);
87	>	}
88	>	evaluator3_.loadScriptString(sele3);
89	>	if (!evaluator3_.isDynamic()) {
90	>	seleMan3_.setSelectionSet(evaluator3_.evaluate());
91	>	}
92	>	}
93
94	<	GofAngle2::GofAngle2(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
95	<	: RadialDistrFunc(info, filename, sele1, sele2){
96	<	setOutputName(getPrefix(filename) + ".gfotw");
94	>	void GofAngle2::processNonOverlapping( SelectionManager& sman1,
95	>	SelectionManager& sman2) {
96	>	StuntDouble* sd1;
97	>	StuntDouble* sd2;
98	>	StuntDouble* sd3;
99	>	int i;
100	>	int j;
101	>	int k;
102	>
103	>	// This is the same as a non-overlapping pairwise loop structure:
104	>	// for (int i = 0;  i < ni ; ++i ) {
105	>	//   for (int j = 0; j < nj; ++j) {}
106	>	// }
107
108	<	}
108	>	if (doSele3_) {
109	>	if  (evaluator3_.isDynamic()) {
110	>	seleMan3_.setSelectionSet(evaluator3_.evaluate());
111	>	}
112	>	if (sman1.getSelectionCount() != seleMan3_.getSelectionCount() ) {
113	>	RadialDistrFunc::processNonOverlapping( sman1, sman2 );
114	>	}
115
116	+	for (sd1 = sman1.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
117	+	sd1 != NULL && sd3 != NULL;
118	+	sd1 = sman1.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
119	+	for (sd2 = sman2.beginSelected(j); sd2 != NULL;
120	+	sd2 = sman2.nextSelected(j)) {
121	+	collectHistogram(sd1, sd2, sd3);
122	+	}
123	+	}
124	+	} else {
125	+	RadialDistrFunc::processNonOverlapping( sman1, sman2 );
126	+	}
127	+	}
128
129	<	void GofAngle2::preProcess() {
129	>	void GofAngle2::processOverlapping( SelectionManager& sman) {
130	>	StuntDouble* sd1;
131	>	StuntDouble* sd2;
132	>	StuntDouble* sd3;
133	>	int i;
134	>	int j;
135	>	int k;
136
137	<	for (int i = 0; i < avgGofr_.size(); ++i) {
138	<	std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
137	>	// This is the same as a pairwise loop structure:
138	>	// for (int i = 0;  i < n-1 ; ++i ) {
139	>	//   for (int j = i + 1; j < n; ++j) {}
140	>	// }
141	>
142	>	if (doSele3_) {
143	>	if  (evaluator3_.isDynamic()) {
144	>	seleMan3_.setSelectionSet(evaluator3_.evaluate());
145	>	}
146	>	if (sman.getSelectionCount() != seleMan3_.getSelectionCount() ) {
147	>	RadialDistrFunc::processOverlapping( sman);
148	>	}
149	>	for (sd1 = sman.beginSelected(i), sd3 = seleMan3_.beginSelected(k);
150	>	sd1 != NULL && sd3 != NULL;
151	>	sd1 = sman.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
152	>	for (j  = i, sd2 = sman.nextSelected(j); sd2 != NULL;
153	>	sd2 = sman.nextSelected(j)) {
154	>	collectHistogram(sd1, sd2, sd3);
155	>	}
156	>	}
157	>	} else {
158	>	RadialDistrFunc::processOverlapping( sman);
159	>	}
160	>	}
161	>
162	>
163	>	void GofAngle2::preProcess() {
164	>
165	>	for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
166	>	std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
167		}
168	<	}
168	>	}
169
170	<	void GofAngle2::initalizeHistogram() {
170	>	void GofAngle2::initializeHistogram() {
171		npairs_ = 0;
172	<	for (int i = 0; i < histogram_.size(); ++i)
173	<	std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
174	<	}
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 GofAngle2::processHistogram() {
177	>	void GofAngle2::processHistogram() {
178
179		//std::for_each(avgGofr_.begin(), avgGofr_.end(), std::plus<std::vector<int>>)
180
181	<	}
181	>	}
182
183	<	void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
183	>	void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
184
185		if (sd1 == sd2) {
186	<	return;
186	>	return;
187		}
188
189		Vector3d pos1 = sd1->getPos();
190		Vector3d pos2 = sd2->getPos();
191		Vector3d r12 = pos1 - pos2;
192	<	currentSnapshot_->wrapVector(r12);
193	<	Vector3d dipole1 = sd1->getElectroFrame().getColumn(2);
194	<	Vector3d dipole2 = sd2->getElectroFrame().getColumn(2);
192	>	if (usePeriodicBoundaryConditions_)
193	>	currentSnapshot_->wrapVector(r12);
194	>
195	>	AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
196	>	AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
197	>	MultipoleAdapter ma1 = MultipoleAdapter(atype1);
198	>	MultipoleAdapter ma2 = MultipoleAdapter(atype2);
199	>
200	>	if (!sd1->isDirectional()) {
201	>	sprintf(painCave.errMsg,
202	>	"GofAngle2: attempted to use a non-directional object: %s\n",
203	>	sd1->getType().c_str());
204	>	painCave.isFatal = 1;
205	>	simError();
206	>	}
207	>
208	>	if (!sd2->isDirectional()) {
209	>	sprintf(painCave.errMsg,
210	>	"GofAngle2: attempted to use a non-directional object: %s\n",
211	>	sd2->getType().c_str());
212	>	painCave.isFatal = 1;
213	>	simError();
214	>	}
215	>
216	>	Vector3d dipole1, dipole2;
217	>	if (ma1.isDipole())
218	>	dipole1 = sd1->getDipole();
219	>	else
220	>	dipole1 = sd1->getA().transpose() * V3Z;
221	>
222	>	if (ma2.isDipole())
223	>	dipole2 = sd2->getDipole();
224	>	else
225	>	dipole2 = sd2->getA().transpose() * V3Z;
226
227		r12.normalize();
228		dipole1.normalize();
229		dipole2.normalize();
230
231
232	<	double cosAngle1 = dot(r12, dipole1);
233	<	double cosAngle2 = dot(dipole1, dipole2);
232	>	RealType cosAngle1 = dot(r12, dipole1);
233	>	RealType cosAngle2 = dot(dipole1, dipole2);
234
235	<	double halfBin = (nAngleBins_ - 1) * 0.5;
236	<	int angleBin1 = halfBin * (cosAngle1 + 1.0);
237	<	int angleBin2 = halfBin * (cosAngle1 + 1.0);
235	>	RealType halfBin = (nAngleBins_ - 1) * 0.5;
236	>	int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
237	>	int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
238
239	<	++histogram_[angleBin1][angleBin1];
239	>	++histogram_[angleBin1][angleBin2];
240		++npairs_;
241	<	}
241	>	}
242
243	<	void GofAngle2::writeRdf() {
243	>	void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2,
244	>	StuntDouble* sd3) {
245	>
246	>	if (sd1 == sd2) {
247	>	return;
248	>	}
249	>
250	>	Vector3d p1 = sd1->getPos();
251	>	Vector3d p3 = sd3->getPos();
252	>
253	>	Vector3d c = 0.5 * (p1 + p3);
254	>	Vector3d r13 = p3 - p1;
255	>
256	>	Vector3d r12 = sd2->getPos() - c;
257	>
258	>	if (usePeriodicBoundaryConditions_) {
259	>	currentSnapshot_->wrapVector(r12);
260	>	currentSnapshot_->wrapVector(r13);
261	>	}
262	>	r12.normalize();
263	>	r13.normalize();
264	>
265	>	if (!sd2->isDirectional()) {
266	>	sprintf(painCave.errMsg,
267	>	"GofAngle2: attempted to use a non-directional object: %s\n",
268	>	sd2->getType().c_str());
269	>	painCave.isFatal = 1;
270	>	simError();
271	>	}
272	>
273	>	AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
274	>	MultipoleAdapter ma2 = MultipoleAdapter(atype2);
275	>
276	>	Vector3d dipole2;
277	>	if (ma2.isDipole())
278	>	dipole2 = sd2->getDipole();
279	>	else
280	>	dipole2 = sd2->getA().transpose() * V3Z;
281	>
282	>	dipole2.normalize();
283	>
284	>	RealType cosAngle1 = dot(r12, r13);
285	>	RealType cosAngle2 = dot(r13, dipole2);
286	>
287	>	RealType halfBin = (nAngleBins_ - 1) * 0.5;
288	>	int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
289	>	int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
290	>
291	>	++histogram_[angleBin1][angleBin2];
292	>	++npairs_;
293	>
294	>	}
295	>
296	>	void GofAngle2::writeRdf() {
297		std::ofstream rdfStream(outputFilename_.c_str());
298		if (rdfStream.is_open()) {
299	<	rdfStream << "#radial distribution function\n";
300	<	rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
301	<	rdfStream << "selection2: (" << selectionScript2_ << ")\n";
302	<	rdfStream << "#r\tcorrValue\n";
303	<	for (int i = 0; i < avgGofr_.size(); ++i) {
304	<	double cosAngle1 = -1.0 + (i + 0.5)*deltaCosAngle_;
305	<
306	<	for(int j = 0; j < avgGofr_[i].size(); ++j) {
307	<	double cosAngle2 = -1.0 + (i + 0.5)*deltaCosAngle_;
308	<	rdfStream << cosAngle1 << "\t" << cosAngle2 << "\t" << avgGofr_[i][j]/nProcessed_ << "\n";
309	<	}
310	<	}
299	>	rdfStream << "#radial distribution function\n";
300	>	rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
301	>	rdfStream << "selection2: (" << selectionScript2_ << ")";
302	>	if (doSele3_) {
303	>	rdfStream << "\tselection3: (" << selectionScript3_ << ")\n";
304	>	} else {
305	>	rdfStream << "\n";
306	>	}
307	>	rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = "
308	>	<< deltaCosAngle_ << "\n";
309	>	for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
310	>	// RealType cosAngle1 = -1.0 + (i + 0.5)*deltaCosAngle_;
311
312	+	for(unsigned int j = 0; j < avgGofr_[i].size(); ++j) {
313	+	// RealType cosAngle2 = -1.0 + (j + 0.5)*deltaCosAngle_;
314	+	rdfStream <<avgGofr_[i][j]/nProcessed_ << "\t";
315	+	}
316	+	rdfStream << "\n";
317	+	}
318	+
319		} else {
320
321	<
321	>	sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n",
322	>	outputFilename_.c_str());
323	>	painCave.isFatal = 1;
324	>	simError();
325		}
326
327		rdfStream.close();
328	<	}
328	>	}
329
330		}

Comparing trunk/src/applications/staticProps/GofAngle2.cpp (property svn:keywords):
Revision 353 by tim, Wed Feb 16 19:36:30 2005 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

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

Diff Legend

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