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

Comparing trunk/OOPSE-4/src/applications/staticProps/Hxy.cpp (file contents):
Revision 2751 by xsun, Fri May 12 21:34:43 2006 UTC vs.
Revision 2752 by gezelter, Tue May 16 02:06:37 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.2 2006-05-16 02:06:37 gezelter 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 < #include<mkl_lapack64.h>
65 > #endif
66  
67   namespace oopse {
68    
# Line 76 | Line 77 | namespace oopse {
77      gridsample_.resize(nBinsX_*nBinsY_);
78      gridZ_.resize(nBinsX_*nBinsY_);
79  
80 <    sum_bin.resize(nbins);
81 <    avg_bin.resize(nbins);
82 <    errbin_sum.resize(nbins);
83 <    errbin.resize(nbins);
84 <    sum_bin_sq.resize(nbins);
85 <    avg_bin_sq.resize(nbins);
86 <    errbin_sum_sq.resize(nbins);
87 <    errbin_sq.resize(nbins);
80 >    sum_bin.resize(nbins_);
81 >    avg_bin.resize(nbins_);
82 >    errbin_sum.resize(nbins_);
83 >    errbin.resize(nbins_);
84 >    sum_bin_sq.resize(nbins_);
85 >    avg_bin_sq.resize(nbins_);
86 >    errbin_sum_sq.resize(nbins_);
87 >    errbin_sq.resize(nbins_);
88  
89      setOutputName(getPrefix(filename) + ".Hxy");
90    }
91  
92    void Hxy::process() {
93 + #ifndef WITHOUT_FFTW  
94 +
95      DumpReader reader(info_, dumpFilename_);    
96      int nFrames = reader.getNFrames();
97      nProcessed_ = nFrames/step_;
98 <
98 >    
99      std::vector<double> mag, newmag;
100      double lenX_, lenY_;
101      double gridX_, gridY_;
# Line 105 | Line 108 | namespace oopse {
108      double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
109      double maxfreqx, maxfreqy, maxfreq, dfreq;
110      int whichbin;
111 +    int nMolecules;
112      
113      for (int istep = 0; istep < nFrames; istep += step_) {
114        
115 <      int nMolecules = reader.getNMolecules();
115 >      reader.readFrame(istep);
116 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
117 >      nMolecules = info_->getNGlobalMolecules();
118  
119 +      Mat3x3d hmat = currentSnapshot_->getHmat();
120 +
121 +      
122        fftw_complex *in, *out;
123        fftw_plan p;
124 <
125 <      in = fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
126 <      out = fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
124 >      
125 >      in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
126 >      out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
127        p =  fftw_plan_dft_2d(nBinsX_,
128                              nBinsY_,
129                              in, out,
# Line 122 | Line 131 | namespace oopse {
131                              FFTW_ESTIMATE);
132        
133        int i, j;  
134 <
135 <      for(i=0; i < nBinsX_*nBinsY_; i++){
136 <        gridsample_[i].clear();
137 <        gridZ_[i].clear();
138 <      }
139 <
140 <      for(i=0; i < nbins; i++){
141 <        sum_bin[i].clear();
142 <        avg_bin[i].clear();
143 <        errbin_sum[i].clear();
144 <        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 <      }
134 >      
135 >      gridsample_.clear();
136 >      gridZ_.clear();
137 >      sum_bin.clear();
138 >      avg_bin.clear();
139 >      errbin_sum.clear();
140 >      errbin.clear();
141 >      sum_bin_sq.clear();
142 >      avg_bin_sq.clear();
143 >      errbin_sum_sq.clear();
144 >      errbin_sq.clear();
145        
146 +      mag.resize(nBinsX_*nBinsY_);
147 +      newmag.resize(nBinsX_*nBinsY_);    
148 +      mag.clear();
149 +      newmag.clear();
150 +      
151        StuntDouble* sd;
152 <      reader.readFrame(istep);
152 <      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
153 <        
154 <      Mat3x3d hmat = currentSnapshot_->getHmat();
155 <
152 >      
153        lenX_ = hmat(0,0);
154        lenY_ = hmat(1,1);
155 <
155 >      
156        gridX_ = lenX_ /(nBinsX_);
157        gridY_ = lenY_ /(nBinsY_);
158 <
158 >      
159        double halfBoxX_ = lenX_ / 2.0;      
160        double halfBoxY_ = lenY_ / 2.0;      
161 <        
161 >      
162        if (evaluator_.isDynamic()) {
163          seleMan_.setSelectionSet(evaluator_.evaluate());
164        }
# Line 203 | Line 200 | namespace oopse {
200            }
201          }
202        }
203 <        
203 >      
204        for (i=0; i< nBinsX_; i++) {
205          for(j=0; j< nBinsY_; j++) {
206            newindex = i*nBinsY_ + j;
# Line 248 | Line 245 | namespace oopse {
245              px = i + 1;
246              py = j;
247              newp = px*nBinsY_ + py;
248 <            if ((px < nBinsX_) && (gridsample_[newp] > 0)) {
249 <              interpsum += gridZ_[newp];
250 <              ninterp++;
251 <            }
252 <            
248 >            if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
249 >              interpsum += gridZ_[newp];
250 >              ninterp++;
251 >            }
252 >        
253              value = interpsum / (double)ninterp;
254              
255              gridZ_[newindex] = value;
# Line 328 | Line 325 | namespace oopse {
325        //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
326        
327        maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
328 <      dfreq = maxfreq/(double)(nbins-1);
328 >      dfreq = maxfreq/(double)(nbins_-1);
329      
330        //printf("%lf\n", dfreq);
331        
# Line 351 | Line 348 | namespace oopse {
348          }
349        }
350        
351 <      for ( i = 0; i < nbins; i++) {
351 >      for ( i = 0; i < nbins_; i++) {
352          if ( samples[i][istep] > 0) {
353            bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
354          }
# Line 359 | Line 356 | namespace oopse {
356        
357      }
358    
359 <    for (i = 0; i < nbins; i++) {
360 <      for (j = 0; j < nFrames; j++) {
359 >    for (int i = 0; i < nbins_; i++) {
360 >      for (int j = 0; j < nFrames; j++) {
361          sum_bin[i] += bin[i][j];
362          sum_bin_sq[i] += bin[i][j] * bin[i][j];
363        }
364        avg_bin[i] = sum_bin[i] / (double)nFrames;
365        avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
366 <      for (j = 0; j < nFrames; j++) {
366 >      for (int j = 0; j < nFrames; j++) {
367          errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
368          errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
369        }
370        errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
371        errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
372      }
373 <
373 >    
374      printSpectrum();
375 <  }
375 > #else
376 >    sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
377 >    painCave.isFatal = 1;
378 >    simError();  
379 >
380 > #endif
381 >
382 > }
383      
384    void Hxy::printSpectrum() {
385      std::ofstream rdfStream(outputFilename_.c_str());
386      if (rdfStream.is_open()) {
387        
388 <      for (int i = 0; i < nbins; i++) {
388 >      for (int i = 0; i < nbins_; i++) {
389          if ( avg_bin[i] > 0 ){
390            rdfStream << i*dfreq << "\t"
391                      <<pow(avg_bin[i], 2)<<"\t"

Diff Legend

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