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 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/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, const std::string& sele1, const std::string& sele2)
53	<	: RadialDistrFunc(info, filename, sele1, sele2){
54	<	setOutputName(getPrefix(filename) + ".gfotw");
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	>	evaluator3_(info),
57	>	seleMan3_(info), doSele3_(false) {
58	>
59	>	setOutputName(getPrefix(filename) + ".gto");
60	>
61	>	deltaCosAngle_ = 2.0 / nAngleBins_;
62	>
63	>	histogram_.resize(nAngleBins_);
64	>	avgGofr_.resize(nAngleBins_);
65	>	for (int i = 0 ; i < nAngleBins_; ++i) {
66	>	histogram_[i].resize(nAngleBins_);
67	>	avgGofr_[i].resize(nAngleBins_);
68	>	}
69	>	}
70
71	<	}
71	>	GofAngle2::GofAngle2(SimInfo* info, const std::string& filename,
72	>	const std::string& sele1,
73	>	const std::string& sele2,
74	>	const std::string& sele3, int nangleBins)
75	>	: RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins),
76	>	evaluator3_(info), selectionScript3_(sele3),
77	>	seleMan3_(info), doSele3_(true) {
78	>
79	>	setOutputName(getPrefix(filename) + ".gto");
80
81	+	deltaCosAngle_ = 2.0 / nAngleBins_;
82	+
83	+	histogram_.resize(nAngleBins_);
84	+	avgGofr_.resize(nAngleBins_);
85	+	for (int i = 0 ; i < nAngleBins_; ++i) {
86	+	histogram_[i].resize(nAngleBins_);
87	+	avgGofr_[i].resize(nAngleBins_);
88	+	}
89	+	evaluator3_.loadScriptString(sele3);
90	+	if (!evaluator3_.isDynamic()) {
91	+	seleMan3_.setSelectionSet(evaluator3_.evaluate());
92	+	}
93	+	}
94
95	<	void GofAngle2::preProcess() {
95	>	void GofAngle2::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	<	for (int i = 0; i < avgGofr_.size(); ++i) {
110	<	std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
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	<	}
128	>	}
129
130	<	void GofAngle2::initalizeHistogram() {
130	>	void GofAngle2::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 GofAngle2::preProcess() {
165	>
166	>	for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
167	>	std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
168	>	}
169	>	}
170	>
171	>	void GofAngle2::initializeHistogram() {
172		npairs_ = 0;
173	<	for (int i = 0; i < histogram_.size(); ++i)
174	<	std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
175	<	}
173	>	for (unsigned int i = 0; i < histogram_.size(); ++i)
174	>	std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
175	>	}
176
177
178	<	void GofAngle2::processHistogram() {
178	>	void GofAngle2::processHistogram() {
179
180		//std::for_each(avgGofr_.begin(), avgGofr_.end(), std::plus<std::vector<int>>)
181
182	<	}
182	>	}
183
184	<	void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
184	>	void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
185
186		if (sd1 == sd2) {
187	<	return;
187	>	return;
188		}
189
190		Vector3d pos1 = sd1->getPos();
191		Vector3d pos2 = sd2->getPos();
192		Vector3d r12 = pos1 - pos2;
193	<	currentSnapshot_->wrapVector(r12);
194	<	Vector3d dipole1 = sd1->getElectroFrame().getColumn(2);
195	<	Vector3d dipole2 = sd2->getElectroFrame().getColumn(2);
193	>	if (usePeriodicBoundaryConditions_)
194	>	currentSnapshot_->wrapVector(r12);
195	>
196	>	AtomType* atype1 = static_cast<Atom*>(sd1)->getAtomType();
197	>	AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
198	>	MultipoleAdapter ma1 = MultipoleAdapter(atype1);
199	>	MultipoleAdapter ma2 = MultipoleAdapter(atype2);
200	>
201	>	if (!sd1->isDirectional()) {
202	>	sprintf(painCave.errMsg,
203	>	"GofAngle2: attempted to use a non-directional object: %s\n",
204	>	sd1->getType().c_str());
205	>	painCave.isFatal = 1;
206	>	simError();
207	>	}
208	>
209	>	if (!sd2->isDirectional()) {
210	>	sprintf(painCave.errMsg,
211	>	"GofAngle2: attempted to use a non-directional object: %s\n",
212	>	sd2->getType().c_str());
213	>	painCave.isFatal = 1;
214	>	simError();
215	>	}
216	>
217	>	Vector3d dipole1, dipole2;
218	>	if (ma1.isDipole())
219	>	dipole1 = sd1->getDipole();
220	>	else
221	>	dipole1 = sd1->getA().transpose() * V3Z;
222	>
223	>	if (ma2.isDipole())
224	>	dipole2 = sd2->getDipole();
225	>	else
226	>	dipole2 = sd2->getA().transpose() * V3Z;
227
228		r12.normalize();
229		dipole1.normalize();
230		dipole2.normalize();
231
232
233	<	double cosAngle1 = dot(r12, dipole1);
234	<	double cosAngle2 = dot(dipole1, dipole2);
233	>	RealType cosAngle1 = dot(r12, dipole1);
234	>	RealType cosAngle2 = dot(dipole1, dipole2);
235
236	<	double halfBin = (nAngleBins_ - 1) * 0.5;
237	<	int angleBin1 = halfBin * (cosAngle1 + 1.0);
238	<	int angleBin2 = halfBin * (cosAngle1 + 1.0);
236	>	RealType halfBin = (nAngleBins_ - 1) * 0.5;
237	>	int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
238	>	int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
239
240	<	++histogram_[angleBin1][angleBin1];
240	>	++histogram_[angleBin1][angleBin2];
241		++npairs_;
242	<	}
242	>	}
243
244	<	void GofAngle2::writeRdf() {
244	>	void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2,
245	>	StuntDouble* sd3) {
246	>
247	>	if (sd1 == sd2) {
248	>	return;
249	>	}
250	>
251	>	Vector3d p1 = sd1->getPos();
252	>	Vector3d p3 = sd3->getPos();
253	>
254	>	Vector3d c = 0.5 * (p1 + p3);
255	>	Vector3d r13 = p3 - p1;
256	>
257	>	Vector3d r12 = sd2->getPos() - c;
258	>
259	>	if (usePeriodicBoundaryConditions_) {
260	>	currentSnapshot_->wrapVector(r12);
261	>	currentSnapshot_->wrapVector(r13);
262	>	}
263	>	r12.normalize();
264	>	r13.normalize();
265	>
266	>	if (!sd2->isDirectional()) {
267	>	sprintf(painCave.errMsg,
268	>	"GofAngle2: attempted to use a non-directional object: %s\n",
269	>	sd2->getType().c_str());
270	>	painCave.isFatal = 1;
271	>	simError();
272	>	}
273	>
274	>	AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
275	>	MultipoleAdapter ma2 = MultipoleAdapter(atype2);
276	>
277	>	Vector3d dipole2;
278	>	if (ma2.isDipole())
279	>	dipole2 = sd2->getDipole();
280	>	else
281	>	dipole2 = sd2->getA().transpose() * V3Z;
282	>
283	>	dipole2.normalize();
284	>
285	>	RealType cosAngle1 = dot(r12, r13);
286	>	RealType cosAngle2 = dot(r13, dipole2);
287	>
288	>	RealType halfBin = (nAngleBins_ - 1) * 0.5;
289	>	int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
290	>	int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
291	>
292	>	++histogram_[angleBin1][angleBin2];
293	>	++npairs_;
294	>
295	>	}
296	>
297	>	void GofAngle2::writeRdf() {
298		std::ofstream rdfStream(outputFilename_.c_str());
299		if (rdfStream.is_open()) {
300	<	rdfStream << "#radial distribution function\n";
301	<	rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
302	<	rdfStream << "selection2: (" << selectionScript2_ << ")\n";
303	<	rdfStream << "#r\tcorrValue\n";
304	<	for (int i = 0; i < avgGofr_.size(); ++i) {
305	<	double cosAngle1 = -1.0 + (i + 0.5)*deltaCosAngle_;
306	<
307	<	for(int j = 0; j < avgGofr_[i].size(); ++j) {
308	<	double cosAngle2 = -1.0 + (i + 0.5)*deltaCosAngle_;
309	<	rdfStream << cosAngle1 << "\t" << cosAngle2 << "\t" << avgGofr_[i][j]/nProcessed_ << "\n";
310	<	}
311	<	}
300	>	rdfStream << "#radial distribution function\n";
301	>	rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
302	>	rdfStream << "selection2: (" << selectionScript2_ << ")";
303	>	if (doSele3_) {
304	>	rdfStream << "\tselection3: (" << selectionScript3_ << ")\n";
305	>	} else {
306	>	rdfStream << "\n";
307	>	}
308	>	rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = "
309	>	<< deltaCosAngle_ << "\n";
310	>	for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
311	>	// RealType cosAngle1 = -1.0 + (i + 0.5)*deltaCosAngle_;
312
313	+	for(unsigned int j = 0; j < avgGofr_[i].size(); ++j) {
314	+	// RealType cosAngle2 = -1.0 + (j + 0.5)*deltaCosAngle_;
315	+	rdfStream <<avgGofr_[i][j]/nProcessed_ << "\t";
316	+	}
317	+	rdfStream << "\n";
318	+	}
319	+
320		} else {
321
322	<
322	>	sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n",
323	>	outputFilename_.c_str());
324	>	painCave.isFatal = 1;
325	>	simError();
326		}
327
328		rdfStream.close();
329	<	}
329	>	}
330
331		}

Comparing trunk/src/applications/staticProps/GofAngle2.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