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 956 by gezelter, Tue May 16 02:06:37 2006 UTC vs.
Revision 967 by xsun, Mon May 22 15:30:42 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.2 2006-05-16 02:06:37 gezelter Exp $
47 > *  @version $Id: Hxy.cpp,v 1.5 2006-05-22 15:30:42 xsun Exp $
48   *
49   */
50  
# Line 60 | Line 60
60   #include<string.h>
61   #include<stdlib.h>
62   #include<math.h>
63 – #ifndef WITHOUT_FFTW
64 – #include<fftw3.h>
65 – #endif
63  
64   namespace oopse {
65    
# Line 76 | 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 86 | 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 <  void Hxy::process() {
95 < #ifndef WITHOUT_FFTW  
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<double> mag, newmag;
122 <    double lenX_, lenY_;
123 <    double gridX_, gridY_;
124 <    double halfBoxX_, halfBoxY_;
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 <    double interpsum, value;
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 <    double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
136 <    double maxfreqx, maxfreqy, maxfreq, dfreq;
135 >
136 >    RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
137 >    RealType maxfreqx, maxfreqy, maxfreq;
138 >
139      int whichbin;
140      int nMolecules;
141 <    
141 >
142      for (int istep = 0; istep < nFrames; istep += step_) {
143        
144        reader.readFrame(istep);
145        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
146        nMolecules = info_->getNGlobalMolecules();
147 <
147 >      
148        Mat3x3d hmat = currentSnapshot_->getHmat();
120 –
149        
150 <      fftw_complex *in, *out;
150 > #ifdef HAVE_FFTW3_H
151        fftw_plan p;
152 + #else
153 +      fftwnd_plan p;
154 + #endif
155 +      fftw_complex *in, *out;
156        
157        in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
158        out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
159 <      p =  fftw_plan_dft_2d(nBinsX_,
160 <                            nBinsY_,
161 <                            in, out,
162 <                            FFTW_FORWARD,
163 <                            FFTW_ESTIMATE);
159 >
160 > #ifdef HAVE_FFTW3_H
161 >      p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,
162 >                           FFTW_FORWARD, FFTW_ESTIMATE);
163 > #else
164 >      p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
165 > #endif
166        
167        int i, j;  
168        
169 <      gridsample_.clear();
170 <      gridZ_.clear();
171 <      sum_bin.clear();
172 <      avg_bin.clear();
173 <      errbin_sum.clear();
174 <      errbin.clear();
175 <      sum_bin_sq.clear();
176 <      avg_bin_sq.clear();
177 <      errbin_sum_sq.clear();
178 <      errbin_sq.clear();
179 <      
180 <      mag.resize(nBinsX_*nBinsY_);
181 <      newmag.resize(nBinsX_*nBinsY_);    
182 <      mag.clear();
183 <      newmag.clear();
169 >      std::fill(gridsample_.begin(), gridsample_.end(), 0);
170 >      std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
171 >      std::fill(sum_bin.begin(), sum_bin.end(), 0.0);
172 >      std::fill(avg_bin.begin(), avg_bin.end(), 0.0);
173 >      std::fill(errbin_sum.begin(), errbin_sum.end(), 0.0);
174 >      std::fill(errbin.begin(), errbin.end(), 0.0);
175 >      std::fill(sum_bin_sq.begin(), sum_bin_sq.end(), 0.0);
176 >      std::fill(avg_bin_sq.begin(), avg_bin_sq.end(), 0.0);
177 >      std::fill(errbin_sum_sq.begin(), errbin_sum_sq.end(), 0.0);
178 >      std::fill(errbin_sq.begin(), errbin_sq.end(), 0.0);
179 >      std::fill(mag.begin(), mag.end(), 0.0);
180 >      std::fill(newmag.begin(), newmag.end(), 0.0);
181 >
182 >      for(i=0; i < bin.size(); i++)
183 >        std::fill(bin[i].begin(), bin[i].end(), 0.0);
184 >
185 >      for(i=0; i < samples.size(); i++)
186 >        std::fill(samples[i].begin(), samples[i].end(), 0);
187        
188        StuntDouble* sd;
189        
# Line 156 | Line 193 | namespace oopse {
193        gridX_ = lenX_ /(nBinsX_);
194        gridY_ = lenY_ /(nBinsY_);
195        
196 <      double halfBoxX_ = lenX_ / 2.0;      
197 <      double halfBoxY_ = lenY_ / 2.0;      
196 >      halfBoxX_ = lenX_ / 2.0;      
197 >      halfBoxY_ = lenY_ / 2.0;      
198        
199        if (evaluator_.isDynamic()) {
200          seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 168 | Line 205 | namespace oopse {
205          Vector3d pos = sd->getPos();
206          currentSnapshot_->wrapVector(pos);
207          sd->setPos(pos);
208 <      }
208 >      }
209        
210        //determine which atom belongs to which grid
211        for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
212          Vector3d pos = sd->getPos();
213          //int binNo = (pos.z() /deltaR_) - 1;
214 <        int binNoX = (pos.x() + halfBoxX_) /gridX_;
215 <        int binNoY = (pos.y() + halfBoxY_) /gridY_;
214 >        int binNoX = (int) ((pos.x() + halfBoxX_) / gridX_);
215 >        int binNoY = (int) ((pos.y() + halfBoxY_) / gridY_);
216          //std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = "  << deltaR_ << " binNo = " << binNo << "\n";
217          gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
218          gridsample_[binNoX*nBinsY_+binNoY]++;
# Line 184 | Line 221 | namespace oopse {
221        // FFT stuff depends on nx and ny, so delay allocation until we have
222        // that information
223        
187 –      for (i=0; i< nBinsX_; i++) {
188 –        for(j=0; j< nBinsY_; j++) {
189 –          newindex = i*nBinsY_ + j;
190 –          mag[newindex] = 0.0;
191 –          newmag[newindex] = 0.0;
192 –        }
193 –      }
194 –      
224        for(i = 0; i < nBinsX_; i++){
225          for(j = 0; j < nBinsY_; j++){
226            newindex = i * nBinsY_ + j;
227            if(gridsample_[newindex] > 0){
228 <            gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
228 >            gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
229            }
230          }
231        }
# Line 250 | Line 279 | namespace oopse {
279                ninterp++;
280              }
281          
282 <            value = interpsum / (double)ninterp;
282 >            value = interpsum / (RealType)ninterp;
283              
284              gridZ_[newindex] = value;
285            }
# Line 265 | Line 294 | namespace oopse {
294            c_im(in[newindex]) = 0.0;
295          }
296        }
297 +
298 + #ifdef HAVE_FFTW3_H
299 +      fftw_execute(p);
300 + #else
301 +      fftwnd_one(p, in, out);
302 + #endif
303        
269 –      fftw_execute(p);
270 –      
304        for (i=0; i< nBinsX_; i++) {
305          for(j=0; j< nBinsY_; j++) {
306            newindex = i*nBinsY_ + j;
307            mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
308          }
309        }
310 <      
310 >
311 > #ifdef HAVE_FFTW3_H
312        fftw_destroy_plan(p);
313 + #else
314 +      fftwnd_destroy_plan(p);
315 + #endif      
316        fftw_free(out);
317        fftw_free(in);
318  
# Line 325 | Line 362 | namespace oopse {
362        //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
363        
364        maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
365 <      dfreq = maxfreq/(double)(nbins_-1);
365 >      dfreq = maxfreq/(RealType)(nbins_-1);
366      
367        //printf("%lf\n", dfreq);
368        
# Line 335 | Line 372 | namespace oopse {
372        for (i=0; i< nBinsX_; i++) {
373          for(j=0; j< nBinsY_; j++) {
374            
375 <          freq_x = (double)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
376 <          freq_y = (double)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
375 >          freq_x = (RealType)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
376 >          freq_y = (RealType)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
377            
378            freq = sqrt(freq_x*freq_x + freq_y*freq_y);
379            
# Line 350 | Line 387 | namespace oopse {
387        
388        for ( i = 0; i < nbins_; i++) {
389          if ( samples[i][istep] > 0) {
390 <          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
390 >          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nMolecules;
391          }
392        }    
356 –      
393      }
394 <  
394 >
395      for (int i = 0; i < nbins_; i++) {
396        for (int j = 0; j < nFrames; j++) {
397          sum_bin[i] += bin[i][j];
398          sum_bin_sq[i] += bin[i][j] * bin[i][j];
399        }
400 <      avg_bin[i] = sum_bin[i] / (double)nFrames;
401 <      avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
400 >      avg_bin[i] = sum_bin[i] / (RealType)nFrames;
401 >      avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
402        for (int j = 0; j < nFrames; j++) {
403          errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
404          errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
405        }
406 <      errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
407 <      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
406 >      errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
407 >      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
408      }
409 <    
409 >
410      printSpectrum();
411 +
412   #else
413      sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
414      painCave.isFatal = 1;
415      simError();  
416  
417   #endif
418 <
382 < }
418 >  }
419      
420    void Hxy::printSpectrum() {
421      std::ofstream rdfStream(outputFilename_.c_str());
422      if (rdfStream.is_open()) {
423 <      
424 <      for (int i = 0; i < nbins_; i++) {
423 >
424 >      for (int i = 0; i < nbins_; ++i) {
425          if ( avg_bin[i] > 0 ){
426 <          rdfStream << i*dfreq << "\t"
426 >          rdfStream << (RealType)i * dfreq << "\t"
427                      <<pow(avg_bin[i], 2)<<"\t"
428                      <<errbin_sq[i]<<"\t"
429                      <<avg_bin[i]<<"\t"
# Line 400 | Line 436 | namespace oopse {
436        painCave.isFatal = 1;
437        simError();  
438      }
439 <    
439 >
440      rdfStream.close();
441 +
442    }
443    
444   }

Diff Legend

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