ViewVC Help

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

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

Comparing trunk/src/applications/staticProps/Hxy.cpp (file contents):
Revision 955 by xsun, Fri May 12 21:34:43 2006 UTC vs.
Revision 1796 by gezelter, Mon Sep 10 18:38:44 2012 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 38 | Line 29
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, 24107 (2008).
39	+	* [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41		*
42	–	*  Hxy.cpp
43	–	*  OOPSE-2.0
44	–	*
42		*  Created by Xiuquan Sun on 05/09/06.
43		*  @author  Xiuquan Sun
44	<	*  @version $Id: Hxy.cpp,v 1.1 2006-05-12 21:34:43 xsun Exp $
44	>	*  @version $Id$
45		*
46		*/
47
#	Line 60 | Line 57
57		#include<string.h>
58		#include<stdlib.h>
59		#include<math.h>
63	–	#include<fftw3.h>
64	–	#include<mkl_lapack64.h>
60
61	<	namespace oopse {
61	>	namespace OpenMD {
62
63		Hxy::Hxy(SimInfo* info, const std::string& filename, const std::string& sele, int nbins_x, int nbins_y, int nrbins)
64		: StaticAnalyser(info, filename), selectionScript_(sele),  evaluator_(info), seleMan_(info), nBinsX_(nbins_x), nBinsY_(nbins_y), nbins_(nrbins){
#	Line 75 | Line 70 | namespace oopse {
70
71		gridsample_.resize(nBinsX_*nBinsY_);
72		gridZ_.resize(nBinsX_*nBinsY_);
73	+	mag.resize(nBinsX_*nBinsY_);
74	+	newmag.resize(nBinsX_*nBinsY_);
75
76	<	sum_bin.resize(nbins);
77	<	avg_bin.resize(nbins);
78	<	errbin_sum.resize(nbins);
79	<	errbin.resize(nbins);
80	<	sum_bin_sq.resize(nbins);
81	<	avg_bin_sq.resize(nbins);
82	<	errbin_sum_sq.resize(nbins);
83	<	errbin_sq.resize(nbins);
76	>	sum_bin.resize(nbins_);
77	>	avg_bin.resize(nbins_);
78	>	errbin_sum.resize(nbins_);
79	>	errbin.resize(nbins_);
80	>	sum_bin_sq.resize(nbins_);
81	>	avg_bin_sq.resize(nbins_);
82	>	errbin_sum_sq.resize(nbins_);
83	>	errbin_sq.resize(nbins_);
84
85	+	bin.resize(nbins_);
86	+	samples.resize(nbins_);
87	+
88		setOutputName(getPrefix(filename) + ".Hxy");
89		}
90
91	+	Hxy::~Hxy(){
92	+	gridsample_.clear();
93	+	gridZ_.clear();
94	+	sum_bin.clear();
95	+	avg_bin.clear();
96	+	errbin_sum.clear();
97	+	errbin.clear();
98	+	sum_bin_sq.clear();
99	+	avg_bin_sq.clear();
100	+	errbin_sum_sq.clear();
101	+	errbin_sq.clear();
102	+
103	+	for(unsigned int i=0; i < bin.size(); i++)
104	+	bin[i].clear();
105	+	for(unsigned int i=0; i < samples.size(); i++)
106	+	samples[i].clear();
107	+
108	+	mag.clear();
109	+	newmag.clear();
110	+	}
111	+
112		void Hxy::process() {
113	+	#if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
114		DumpReader reader(info_, dumpFilename_);
115		int nFrames = reader.getNFrames();
116		nProcessed_ = nFrames/step_;
117	+
118	+	for(unsigned int k=0; k < bin.size(); k++)
119	+	bin[k].resize(nFrames);
120	+	for(unsigned int k=0; k < samples.size(); k++)
121	+	samples[k].resize(nFrames);
122
123	<	std::vector<double> mag, newmag;
124	<	double lenX_, lenY_;
125	<	double gridX_, gridY_;
126	<	double halfBoxX_, halfBoxY_;
127	<	int binNoX, binNoY;
101	<	double interpsum, value;
123	>	RealType lenX_, lenY_;
124	>	RealType gridX_, gridY_;
125	>	RealType halfBoxX_, halfBoxY_;
126	>
127	>	RealType interpsum, value;
128		int ninterp, px, py, newp;
129	<	int newx, newy, newindex, index;
129	>	int newindex, index;
130		int new_i, new_j, new_index;
131	<	double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
132	<	double maxfreqx, maxfreqy, maxfreq, dfreq;
131	>
132	>	RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
133	>	RealType maxfreqx, maxfreqy, maxfreq;
134	>
135		int whichbin;
136	+
137	+	std::fill(sum_bin.begin(), sum_bin.end(), 0.0);
138	+	std::fill(avg_bin.begin(), avg_bin.end(), 0.0);
139	+	std::fill(errbin_sum.begin(), errbin_sum.end(), 0.0);
140	+	std::fill(errbin.begin(), errbin.end(), 0.0);
141	+	std::fill(sum_bin_sq.begin(), sum_bin_sq.end(), 0.0);
142	+	std::fill(avg_bin_sq.begin(), avg_bin_sq.end(), 0.0);
143	+	std::fill(errbin_sum_sq.begin(), errbin_sum_sq.end(), 0.0);
144	+	std::fill(errbin_sq.begin(), errbin_sq.end(), 0.0);
145
146	+	for(unsigned int i=0; i < bin.size(); i++)
147	+	std::fill(bin[i].begin(), bin[i].end(), 0.0);
148	+
149	+	for(unsigned int i=0; i < samples.size(); i++)
150	+	std::fill(samples[i].begin(), samples[i].end(), 0);
151	+
152		for (int istep = 0; istep < nFrames; istep += step_) {
153
154	<	int nMolecules = reader.getNMolecules();
155	<
156	<	fftw_complex *in, *out;
154	>	reader.readFrame(istep);
155	>	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
156	>
157	>	Mat3x3d hmat = currentSnapshot_->getHmat();
158	>
159	>	#ifdef HAVE_FFTW3_H
160		fftw_plan p;
161	<
162	<	in = fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
163	<	out = fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
164	<	p =  fftw_plan_dft_2d(nBinsX_,
119	<	nBinsY_,
120	<	in, out,
121	<	FFTW_FORWARD,
122	<	FFTW_ESTIMATE);
161	>	#else
162	>	fftwnd_plan p;
163	>	#endif
164	>	fftw_complex *in, *out;
165
166	<	int i, j;
166	>	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
167	>	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
168
169	<	for(i=0; i < nBinsX_*nBinsY_; i++){
170	<	gridsample_[i].clear();
171	<	gridZ_[i].clear();
172	<	}
169	>	#ifdef HAVE_FFTW3_H
170	>	p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,
171	>	FFTW_FORWARD, FFTW_ESTIMATE);
172	>	#else
173	>	p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
174	>	#endif
175
176	<	for(i=0; i < nbins; i++){
177	<	sum_bin[i].clear();
178	<	avg_bin[i].clear();
179	<	errbin_sum[i].clear();
135	<	errbin[i].clear();
136	<	sum_bin_sq[i].clear();
137	<	avg_bin_sq[i].clear();
138	<	errbin_sum_sq[i].clear();
139	<	errbin_sq[i].clear();
140	<	}
176	>	std::fill(gridsample_.begin(), gridsample_.end(), 0);
177	>	std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
178	>	std::fill(mag.begin(), mag.end(), 0.0);
179	>	std::fill(newmag.begin(), newmag.end(), 0.0);
180
181	<	mag.resize(nBinsX_*nBinsY_);
143	<	newmag.resize(nBinsX_*nBinsY_);
144	<
145	<	for(i=0; i < nBinsX_*nBinsY_; i++){
146	<	mag[i].clear();
147	<	newmag[i].clear();
148	<	}
181	>	int i, j;
182
183		StuntDouble* sd;
184	<	reader.readFrame(istep);
152	<	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
153	<
154	<	Mat3x3d hmat = currentSnapshot_->getHmat();
155	<
184	>
185		lenX_ = hmat(0,0);
186		lenY_ = hmat(1,1);
187	<
187	>
188		gridX_ = lenX_ /(nBinsX_);
189		gridY_ = lenY_ /(nBinsY_);
190	<
191	<	double halfBoxX_ = lenX_ / 2.0;
192	<	double halfBoxY_ = lenY_ / 2.0;
193	<
190	>
191	>	halfBoxX_ = lenX_ / 2.0;
192	>	halfBoxY_ = lenY_ / 2.0;
193	>
194		if (evaluator_.isDynamic()) {
195		seleMan_.setSelectionSet(evaluator_.evaluate());
196		}
#	Line 169 | Line 198 | namespace oopse {
198		//wrap the stuntdoubles into a cell
199		for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
200		Vector3d pos = sd->getPos();
201	<	currentSnapshot_->wrapVector(pos);
201	>	if (usePeriodicBoundaryConditions_)
202	>	currentSnapshot_->wrapVector(pos);
203		sd->setPos(pos);
204	<	}
204	>	}
205
206		//determine which atom belongs to which grid
207		for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
208		Vector3d pos = sd->getPos();
209		//int binNo = (pos.z() /deltaR_) - 1;
210	<	int binNoX = (pos.x() + halfBoxX_) /gridX_;
211	<	int binNoY = (pos.y() + halfBoxY_) /gridY_;
210	>	int binNoX = (int) ((pos.x() + halfBoxX_) / gridX_);
211	>	int binNoY = (int) ((pos.y() + halfBoxY_) / gridY_);
212		//std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = "  << deltaR_ << " binNo = " << binNo << "\n";
213		gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
214		gridsample_[binNoX*nBinsY_+binNoY]++;
#	Line 187 | Line 217 | namespace oopse {
217		// FFT stuff depends on nx and ny, so delay allocation until we have
218		// that information
219
190	–	for (i=0; i< nBinsX_; i++) {
191	–	for(j=0; j< nBinsY_; j++) {
192	–	newindex = i*nBinsY_ + j;
193	–	mag[newindex] = 0.0;
194	–	newmag[newindex] = 0.0;
195	–	}
196	–	}
197	–
220		for(i = 0; i < nBinsX_; i++){
221		for(j = 0; j < nBinsY_; j++){
222		newindex = i * nBinsY_ + j;
223		if(gridsample_[newindex] > 0){
224	<	gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
224	>	gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
225		}
226		}
227		}
228	<
228	>
229		for (i=0; i< nBinsX_; i++) {
230		for(j=0; j< nBinsY_; j++) {
231		newindex = i*nBinsY_ + j;
#	Line 248 | Line 270 | namespace oopse {
270		px = i + 1;
271		py = j;
272		newp = px*nBinsY_ + py;
273	<	if ((px < nBinsX_) && (gridsample_[newp] > 0)) {
274	<	interpsum += gridZ_[newp];
275	<	ninterp++;
276	<	}
277	<
278	<	value = interpsum / (double)ninterp;
273	>	if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
274	>	interpsum += gridZ_[newp];
275	>	ninterp++;
276	>	}
277	>
278	>	value = interpsum / (RealType)ninterp;
279
280		gridZ_[newindex] = value;
281		}
#	Line 268 | Line 290 | namespace oopse {
290		c_im(in[newindex]) = 0.0;
291		}
292		}
293	+
294	+	#ifdef HAVE_FFTW3_H
295	+	fftw_execute(p);
296	+	#else
297	+	fftwnd_one(p, in, out);
298	+	#endif
299
272	–	fftw_execute(p);
273	–
300		for (i=0; i< nBinsX_; i++) {
301		for(j=0; j< nBinsY_; j++) {
302		newindex = i*nBinsY_ + j;
303		mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
304		}
305		}
306	<
306	>
307	>	#ifdef HAVE_FFTW3_H
308		fftw_destroy_plan(p);
309	+	#else
310	+	fftwnd_destroy_plan(p);
311	+	#endif
312		fftw_free(out);
313		fftw_free(in);
314
#	Line 328 | Line 358 | namespace oopse {
358		//  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
359
360		maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
361	<	dfreq = maxfreq/(double)(nbins-1);
361	>	dfreq = maxfreq/(RealType)(nbins_-1);
362
363		//printf("%lf\n", dfreq);
364
#	Line 338 | Line 368 | namespace oopse {
368		for (i=0; i< nBinsX_; i++) {
369		for(j=0; j< nBinsY_; j++) {
370
371	<	freq_x = (double)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
372	<	freq_y = (double)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
371	>	freq_x = (RealType)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
372	>	freq_y = (RealType)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
373
374		freq = sqrt(freq_x*freq_x + freq_y*freq_y);
375
#	Line 351 | Line 381 | namespace oopse {
381		}
382		}
383
384	<	for ( i = 0; i < nbins; i++) {
384	>	for ( i = 0; i < nbins_; i++) {
385		if ( samples[i][istep] > 0) {
386	<	bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
386	>	bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nBinsX_ / (RealType)nBinsY_;
387		}
388		}
359	–
389		}
390	<
391	<	for (i = 0; i < nbins; i++) {
392	<	for (j = 0; j < nFrames; j++) {
390	>
391	>	for (int i = 0; i < nbins_; i++) {
392	>	for (int j = 0; j < nFrames; j++) {
393		sum_bin[i] += bin[i][j];
394		sum_bin_sq[i] += bin[i][j] * bin[i][j];
395		}
396	<	avg_bin[i] = sum_bin[i] / (double)nFrames;
397	<	avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
398	<	for (j = 0; j < nFrames; j++) {
396	>	avg_bin[i] = sum_bin[i] / (RealType)nFrames;
397	>	avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
398	>	for (int j = 0; j < nFrames; j++) {
399		errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
400		errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
401		}
402	<	errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
403	<	errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
402	>	errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
403	>	errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
404		}
405
406		printSpectrum();
407	+
408	+	#else
409	+	sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
410	+	painCave.isFatal = 1;
411	+	simError();
412	+
413	+	#endif
414		}
415
416		void Hxy::printSpectrum() {
417		std::ofstream rdfStream(outputFilename_.c_str());
418		if (rdfStream.is_open()) {
419	<
420	<	for (int i = 0; i < nbins; i++) {
419	>
420	>	for (int i = 0; i < nbins_; ++i) {
421		if ( avg_bin[i] > 0 ){
422	<	rdfStream << i*dfreq << "\t"
422	>	rdfStream << (RealType)i * dfreq << "\t"
423		<<pow(avg_bin[i], 2)<<"\t"
424		<<errbin_sq[i]<<"\t"
425		<<avg_bin[i]<<"\t"
#	Line 396 | Line 432 | namespace oopse {
432		painCave.isFatal = 1;
433		simError();
434		}
435	<
435	>
436		rdfStream.close();
437	+
438		}
439
440		}

Comparing trunk/src/applications/staticProps/Hxy.cpp (property svn:keywords):
Revision 955 by xsun, Fri May 12 21:34:43 2006 UTC vs.
Revision 1796 by gezelter, Mon Sep 10 18:38:44 2012 UTC

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

Diff Legend

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