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 955 by xsun, Fri May 12 21:34:43 2006 UTC vs.
Revision 963 by tim, Wed May 17 21:51: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.1 2006-05-12 21:34:43 xsun Exp $
47 > *  @version $Id: Hxy.cpp,v 1.4 2006-05-17 21:51:42 tim Exp $
48   *
49   */
50  
# Line 60 | Line 60
60   #include<string.h>
61   #include<stdlib.h>
62   #include<math.h>
63 – #include<fftw3.h>
64 – #include<mkl_lapack64.h>
63  
64   namespace oopse {
65    
# Line 76 | Line 74 | namespace oopse {
74      gridsample_.resize(nBinsX_*nBinsY_);
75      gridZ_.resize(nBinsX_*nBinsY_);
76  
77 <    sum_bin.resize(nbins);
78 <    avg_bin.resize(nbins);
79 <    errbin_sum.resize(nbins);
80 <    errbin.resize(nbins);
81 <    sum_bin_sq.resize(nbins);
82 <    avg_bin_sq.resize(nbins);
83 <    errbin_sum_sq.resize(nbins);
84 <    errbin_sq.resize(nbins);
77 >    sum_bin.resize(nbins_);
78 >    avg_bin.resize(nbins_);
79 >    errbin_sum.resize(nbins_);
80 >    errbin.resize(nbins_);
81 >    sum_bin_sq.resize(nbins_);
82 >    avg_bin_sq.resize(nbins_);
83 >    errbin_sum_sq.resize(nbins_);
84 >    errbin_sq.resize(nbins_);
85  
86      setOutputName(getPrefix(filename) + ".Hxy");
87    }
88  
89    void Hxy::process() {
90 + #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
91      DumpReader reader(info_, dumpFilename_);    
92      int nFrames = reader.getNFrames();
93      nProcessed_ = nFrames/step_;
94 <
95 <    std::vector<double> mag, newmag;
96 <    double lenX_, lenY_;
97 <    double gridX_, gridY_;
98 <    double halfBoxX_, halfBoxY_;
94 >    
95 >    std::vector<RealType> mag, newmag;
96 >    RealType lenX_, lenY_;
97 >    RealType gridX_, gridY_;
98 >    RealType halfBoxX_, halfBoxY_;
99      int binNoX, binNoY;
100 <    double interpsum, value;
100 >    RealType interpsum, value;
101      int ninterp, px, py, newp;
102      int newx, newy, newindex, index;
103      int new_i, new_j, new_index;
104 <    double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
105 <    double maxfreqx, maxfreqy, maxfreq, dfreq;
104 >    RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
105 >    RealType maxfreqx, maxfreqy, maxfreq, dfreq;
106      int whichbin;
107 +    int nMolecules;
108      
109      for (int istep = 0; istep < nFrames; istep += step_) {
110        
111 <      int nMolecules = reader.getNMolecules();
112 <
113 <      fftw_complex *in, *out;
111 >      reader.readFrame(istep);
112 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
113 >      nMolecules = info_->getNGlobalMolecules();
114 >      
115 >      Mat3x3d hmat = currentSnapshot_->getHmat();
116 >      
117 > #ifdef HAVE_FFTW3_H
118        fftw_plan p;
119 + #else
120 +      fftwnd_plan p;
121 + #endif
122 +      fftw_complex *in, *out;
123  
116 –      in = fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
117 –      out = fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
118 –      p =  fftw_plan_dft_2d(nBinsX_,
119 –                            nBinsY_,
120 –                            in, out,
121 –                            FFTW_FORWARD,
122 –                            FFTW_ESTIMATE);
124        
125 <      int i, j;  
125 >      in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
126 >      out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
127  
128 <      for(i=0; i < nBinsX_*nBinsY_; i++){
129 <        gridsample_[i].clear();
130 <        gridZ_[i].clear();
131 <      }
132 <
133 <      for(i=0; i < nbins; i++){
134 <        sum_bin[i].clear();
135 <        avg_bin[i].clear();
136 <        errbin_sum[i].clear();
137 <        errbin[i].clear();
138 <        sum_bin_sq[i].clear();
139 <        avg_bin_sq[i].clear();
140 <        errbin_sum_sq[i].clear();
141 <        errbin_sq[i].clear();
142 <      }
143 <
144 <      mag.resize(nBinsX_*nBinsY_);
145 <      newmag.resize(nBinsX_*nBinsY_);
146 <
145 <      for(i=0; i < nBinsX_*nBinsY_; i++){
146 <        mag[i].clear();
147 <        newmag[i].clear();
148 <      }
128 > #ifdef HAVE_FFTW3_H
129 >      p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,
130 >                           FFTW_FORWARD, FFTW_ESTIMATE);
131 > #else
132 >      p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
133 > #endif
134 >      
135 >      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();
147        
148 +      mag.resize(nBinsX_*nBinsY_);
149 +      newmag.resize(nBinsX_*nBinsY_);    
150 +      mag.clear();
151 +      newmag.clear();
152 +      
153        StuntDouble* sd;
154 <      reader.readFrame(istep);
152 <      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
153 <        
154 <      Mat3x3d hmat = currentSnapshot_->getHmat();
155 <
154 >      
155        lenX_ = hmat(0,0);
156        lenY_ = hmat(1,1);
157 <
157 >      
158        gridX_ = lenX_ /(nBinsX_);
159        gridY_ = lenY_ /(nBinsY_);
160 <
161 <      double halfBoxX_ = lenX_ / 2.0;      
162 <      double halfBoxY_ = lenY_ / 2.0;      
163 <        
160 >      
161 >      RealType halfBoxX_ = lenX_ / 2.0;      
162 >      RealType halfBoxY_ = lenY_ / 2.0;      
163 >      
164        if (evaluator_.isDynamic()) {
165          seleMan_.setSelectionSet(evaluator_.evaluate());
166        }
# Line 199 | Line 198 | namespace oopse {
198          for(j = 0; j < nBinsY_; j++){
199            newindex = i * nBinsY_ + j;
200            if(gridsample_[newindex] > 0){
201 <            gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
201 >            gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
202            }
203          }
204        }
205 <        
205 >      
206        for (i=0; i< nBinsX_; i++) {
207          for(j=0; j< nBinsY_; j++) {
208            newindex = i*nBinsY_ + j;
# Line 248 | Line 247 | namespace oopse {
247              px = i + 1;
248              py = j;
249              newp = px*nBinsY_ + py;
250 <            if ((px < nBinsX_) && (gridsample_[newp] > 0)) {
251 <              interpsum += gridZ_[newp];
252 <              ninterp++;
253 <            }
250 >            if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
251 >              interpsum += gridZ_[newp];
252 >              ninterp++;
253 >            }
254 >        
255 >            value = interpsum / (RealType)ninterp;
256              
256 –            value = interpsum / (double)ninterp;
257 –            
257              gridZ_[newindex] = value;
258            }
259          }
# Line 268 | Line 267 | namespace oopse {
267            c_im(in[newindex]) = 0.0;
268          }
269        }
270 <      
271 <      fftw_execute(p);
270 > #ifdef HAVE_FFTW3_H
271 >      fftw_execute(p);
272 > #else
273 >      fftwnd_one(p, in, out);
274 > #endif
275        
276        for (i=0; i< nBinsX_; i++) {
277          for(j=0; j< nBinsY_; j++) {
# Line 277 | Line 279 | namespace oopse {
279            mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
280          }
281        }
282 <      
282 > #ifdef HAVE_FFTW3_H
283        fftw_destroy_plan(p);
284 + #else
285 +      fftwnd_destroy_plan(p);
286 + #endif      
287        fftw_free(out);
288        fftw_free(in);
289  
# Line 328 | Line 333 | namespace oopse {
333        //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
334        
335        maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
336 <      dfreq = maxfreq/(double)(nbins-1);
336 >      dfreq = maxfreq/(RealType)(nbins_-1);
337      
338        //printf("%lf\n", dfreq);
339        
# Line 338 | Line 343 | namespace oopse {
343        for (i=0; i< nBinsX_; i++) {
344          for(j=0; j< nBinsY_; j++) {
345            
346 <          freq_x = (double)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
347 <          freq_y = (double)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
346 >          freq_x = (RealType)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
347 >          freq_y = (RealType)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
348            
349            freq = sqrt(freq_x*freq_x + freq_y*freq_y);
350            
# Line 351 | Line 356 | namespace oopse {
356          }
357        }
358        
359 <      for ( i = 0; i < nbins; i++) {
359 >      for ( i = 0; i < nbins_; i++) {
360          if ( samples[i][istep] > 0) {
361 <          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
361 >          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nMolecules;
362          }
363        }    
364        
365      }
366    
367 <    for (i = 0; i < nbins; i++) {
368 <      for (j = 0; j < nFrames; j++) {
367 >    for (int i = 0; i < nbins_; i++) {
368 >      for (int j = 0; j < nFrames; j++) {
369          sum_bin[i] += bin[i][j];
370          sum_bin_sq[i] += bin[i][j] * bin[i][j];
371        }
372 <      avg_bin[i] = sum_bin[i] / (double)nFrames;
373 <      avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
374 <      for (j = 0; j < nFrames; j++) {
372 >      avg_bin[i] = sum_bin[i] / (RealType)nFrames;
373 >      avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
374 >      for (int j = 0; j < nFrames; j++) {
375          errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
376          errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
377        }
378 <      errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
379 <      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
378 >      errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
379 >      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
380      }
381 <
381 >    
382      printSpectrum();
383 <  }
383 > #else
384 >    sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
385 >    painCave.isFatal = 1;
386 >    simError();  
387 >
388 > #endif
389 >
390 > }
391      
392    void Hxy::printSpectrum() {
393      std::ofstream rdfStream(outputFilename_.c_str());
394      if (rdfStream.is_open()) {
395        
396 <      for (int i = 0; i < nbins; i++) {
396 >      for (int i = 0; i < nbins_; i++) {
397          if ( avg_bin[i] > 0 ){
398            rdfStream << i*dfreq << "\t"
399                      <<pow(avg_bin[i], 2)<<"\t"

Diff Legend

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