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 957 by gezelter, Tue May 16 20:38:23 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.3 2006-05-16 20:38:23 gezelter 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 <
94 >    
95      std::vector<double> mag, newmag;
96      double lenX_, lenY_;
97      double gridX_, gridY_;
# Line 105 | Line 104 | namespace oopse {
104      double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
105      double 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();
134 <        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 <      }
141 <
142 <      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 <      }
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 <
160 >      
161        double halfBoxX_ = lenX_ / 2.0;      
162        double halfBoxY_ = lenY_ / 2.0;      
163 <        
163 >      
164        if (evaluator_.isDynamic()) {
165          seleMan_.setSelectionSet(evaluator_.evaluate());
166        }
# Line 203 | Line 202 | namespace oopse {
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 <            }
254 <            
250 >            if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
251 >              interpsum += gridZ_[newp];
252 >              ninterp++;
253 >            }
254 >        
255              value = interpsum / (double)ninterp;
256              
257              gridZ_[newindex] = value;
# Line 268 | Line 267 | namespace oopse {
267            c_im(in[newindex]) = 0.0;
268          }
269        }
270 + #ifdef HAVE_FFTW3_H
271 +      fftw_execute(p);
272 + #else
273 +      fftwnd_one(p, in, out);
274 + #endif
275        
272 –      fftw_execute(p);
273 –      
276        for (i=0; i< nBinsX_; i++) {
277          for(j=0; j< nBinsY_; j++) {
278            newindex = i*nBinsY_ + j;
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/(double)(nbins_-1);
337      
338        //printf("%lf\n", dfreq);
339        
# 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;
362          }
# Line 359 | Line 364 | namespace oopse {
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++) {
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 );
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