ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/Hxy.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/Hxy.cpp (file contents):
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 969 by xsun, Tue May 23 21:12:45 2006 UTC

# Line 44 | Line 44
44   *
45   *  Created by Xiuquan Sun on 05/09/06.
46   *  @author  Xiuquan Sun
47 < *  @version $Id: Hxy.cpp,v 1.4 2006-05-17 21:51:42 tim Exp $
47 > *  @version $Id: Hxy.cpp,v 1.6 2006-05-23 21:12:45 xsun Exp $
48   *
49   */
50  
# Line 73 | Line 73 | namespace oopse {
73  
74      gridsample_.resize(nBinsX_*nBinsY_);
75      gridZ_.resize(nBinsX_*nBinsY_);
76 +    mag.resize(nBinsX_*nBinsY_);
77 +    newmag.resize(nBinsX_*nBinsY_);    
78  
79      sum_bin.resize(nbins_);
80      avg_bin.resize(nbins_);
# Line 83 | Line 85 | namespace oopse {
85      errbin_sum_sq.resize(nbins_);
86      errbin_sq.resize(nbins_);
87  
88 +    bin.resize(nbins_);
89 +    samples.resize(nbins_);
90 +
91      setOutputName(getPrefix(filename) + ".Hxy");
92    }
93  
94 +  Hxy::~Hxy(){
95 +      gridsample_.clear();
96 +      gridZ_.clear();
97 +      sum_bin.clear();
98 +      avg_bin.clear();
99 +      errbin_sum.clear();
100 +      errbin.clear();
101 +      sum_bin_sq.clear();
102 +      avg_bin_sq.clear();
103 +      errbin_sum_sq.clear();
104 +      errbin_sq.clear();
105 +      
106 +      for(int i=0; i < bin.size(); i++)
107 +        bin[i].clear();
108 +      for(int i=0; i < samples.size(); i++)
109 +        samples[i].clear();
110 +
111 +      mag.clear();
112 +      newmag.clear();
113 +  }
114 +
115    void Hxy::process() {
116   #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
117      DumpReader reader(info_, dumpFilename_);    
118      int nFrames = reader.getNFrames();
119      nProcessed_ = nFrames/step_;
120      
121 <    std::vector<RealType> mag, newmag;
121 >    for(int k=0; k < bin.size(); k++)
122 >      bin[k].resize(nFrames);
123 >    for(int k=0; k < samples.size(); k++)
124 >      samples[k].resize(nFrames);
125 >
126      RealType lenX_, lenY_;
127      RealType gridX_, gridY_;
128      RealType halfBoxX_, halfBoxY_;
129 +
130      int binNoX, binNoY;
131      RealType interpsum, value;
132      int ninterp, px, py, newp;
133      int newx, newy, newindex, index;
134      int new_i, new_j, new_index;
135 +
136      RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
137 <    RealType maxfreqx, maxfreqy, maxfreq, dfreq;
137 >    RealType maxfreqx, maxfreqy, maxfreq;
138 >
139      int whichbin;
140      int nMolecules;
141 +
142 +    std::fill(sum_bin.begin(), sum_bin.end(), 0.0);
143 +    std::fill(avg_bin.begin(), avg_bin.end(), 0.0);
144 +    std::fill(errbin_sum.begin(), errbin_sum.end(), 0.0);
145 +    std::fill(errbin.begin(), errbin.end(), 0.0);
146 +    std::fill(sum_bin_sq.begin(), sum_bin_sq.end(), 0.0);
147 +    std::fill(avg_bin_sq.begin(), avg_bin_sq.end(), 0.0);
148 +    std::fill(errbin_sum_sq.begin(), errbin_sum_sq.end(), 0.0);
149 +    std::fill(errbin_sq.begin(), errbin_sq.end(), 0.0);
150      
151 +    for(int i=0; i < bin.size(); i++)
152 +      std::fill(bin[i].begin(), bin[i].end(), 0.0);
153 +    
154 +    for(int i=0; i < samples.size(); i++)
155 +      std::fill(samples[i].begin(), samples[i].end(), 0);
156 +    
157      for (int istep = 0; istep < nFrames; istep += step_) {
158        
159        reader.readFrame(istep);
# Line 120 | Line 168 | namespace oopse {
168        fftwnd_plan p;
169   #endif
170        fftw_complex *in, *out;
123 –
171        
172        in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
173        out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
# Line 131 | Line 178 | namespace oopse {
178   #else
179        p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
180   #endif
181 <      
181 >
182 >      std::fill(gridsample_.begin(), gridsample_.end(), 0);
183 >      std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
184 >      std::fill(mag.begin(), mag.end(), 0.0);
185 >      std::fill(newmag.begin(), newmag.end(), 0.0);
186 >
187        int i, j;  
136 –      
137 –      gridsample_.clear();
138 –      gridZ_.clear();
139 –      sum_bin.clear();
140 –      avg_bin.clear();
141 –      errbin_sum.clear();
142 –      errbin.clear();
143 –      sum_bin_sq.clear();
144 –      avg_bin_sq.clear();
145 –      errbin_sum_sq.clear();
146 –      errbin_sq.clear();
188        
148 –      mag.resize(nBinsX_*nBinsY_);
149 –      newmag.resize(nBinsX_*nBinsY_);    
150 –      mag.clear();
151 –      newmag.clear();
152 –      
189        StuntDouble* sd;
190        
191        lenX_ = hmat(0,0);
# Line 158 | Line 194 | namespace oopse {
194        gridX_ = lenX_ /(nBinsX_);
195        gridY_ = lenY_ /(nBinsY_);
196        
197 <      RealType halfBoxX_ = lenX_ / 2.0;      
198 <      RealType halfBoxY_ = lenY_ / 2.0;      
197 >      halfBoxX_ = lenX_ / 2.0;      
198 >      halfBoxY_ = lenY_ / 2.0;      
199        
200        if (evaluator_.isDynamic()) {
201          seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 170 | Line 206 | namespace oopse {
206          Vector3d pos = sd->getPos();
207          currentSnapshot_->wrapVector(pos);
208          sd->setPos(pos);
209 <      }
209 >      }
210        
211        //determine which atom belongs to which grid
212        for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
213          Vector3d pos = sd->getPos();
214          //int binNo = (pos.z() /deltaR_) - 1;
215 <        int binNoX = (pos.x() + halfBoxX_) /gridX_;
216 <        int binNoY = (pos.y() + halfBoxY_) /gridY_;
215 >        int binNoX = (int) ((pos.x() + halfBoxX_) / gridX_);
216 >        int binNoY = (int) ((pos.y() + halfBoxY_) / gridY_);
217          //std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = "  << deltaR_ << " binNo = " << binNo << "\n";
218          gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
219          gridsample_[binNoX*nBinsY_+binNoY]++;
# Line 185 | Line 221 | namespace oopse {
221        
222        // FFT stuff depends on nx and ny, so delay allocation until we have
223        // that information
188 –      
189 –      for (i=0; i< nBinsX_; i++) {
190 –        for(j=0; j< nBinsY_; j++) {
191 –          newindex = i*nBinsY_ + j;
192 –          mag[newindex] = 0.0;
193 –          newmag[newindex] = 0.0;
194 –        }
195 –      }
224        
225        for(i = 0; i < nBinsX_; i++){
226          for(j = 0; j < nBinsY_; j++){
# Line 267 | Line 295 | namespace oopse {
295            c_im(in[newindex]) = 0.0;
296          }
297        }
298 +
299   #ifdef HAVE_FFTW3_H
300        fftw_execute(p);
301   #else
# Line 279 | Line 308 | namespace oopse {
308            mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
309          }
310        }
311 +
312   #ifdef HAVE_FFTW3_H
313        fftw_destroy_plan(p);
314   #else
# Line 358 | Line 388 | namespace oopse {
388        
389        for ( i = 0; i < nbins_; i++) {
390          if ( samples[i][istep] > 0) {
391 <          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nMolecules;
391 >          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nBinsX_ / (RealType)nBinsY_;
392          }
393        }    
364 –      
394      }
395 <  
395 >
396      for (int i = 0; i < nbins_; i++) {
397        for (int j = 0; j < nFrames; j++) {
398          sum_bin[i] += bin[i][j];
# Line 378 | Line 407 | namespace oopse {
407        errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
408        errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
409      }
410 <    
410 >
411      printSpectrum();
412 +
413   #else
414      sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
415      painCave.isFatal = 1;
416      simError();  
417  
418   #endif
419 <
390 < }
419 >  }
420      
421    void Hxy::printSpectrum() {
422      std::ofstream rdfStream(outputFilename_.c_str());
423      if (rdfStream.is_open()) {
424 <      
425 <      for (int i = 0; i < nbins_; i++) {
424 >
425 >      for (int i = 0; i < nbins_; ++i) {
426          if ( avg_bin[i] > 0 ){
427 <          rdfStream << i*dfreq << "\t"
427 >          rdfStream << (RealType)i * dfreq << "\t"
428                      <<pow(avg_bin[i], 2)<<"\t"
429                      <<errbin_sq[i]<<"\t"
430                      <<avg_bin[i]<<"\t"
# Line 408 | Line 437 | namespace oopse {
437        painCave.isFatal = 1;
438        simError();  
439      }
440 <    
440 >
441      rdfStream.close();
442 +
443    }
444    
445   }

Diff Legend

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