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 1078 by gezelter, Wed Oct 18 21:58:48 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.7 2006-10-18 21:58:47 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 75 | 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);
81 <    errbin_sum.resize(nbins);
82 <    errbin.resize(nbins);
83 <    sum_bin_sq.resize(nbins);
84 <    avg_bin_sq.resize(nbins);
85 <    errbin_sum_sq.resize(nbins);
86 <    errbin_sq.resize(nbins);
79 >    sum_bin.resize(nbins_);
80 >    avg_bin.resize(nbins_);
81 >    errbin_sum.resize(nbins_);
82 >    errbin.resize(nbins_);
83 >    sum_bin_sq.resize(nbins_);
84 >    avg_bin_sq.resize(nbins_);
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 +    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 <    std::vector<double> mag, newmag;
127 <    double lenX_, lenY_;
128 <    double gridX_, gridY_;
129 <    double halfBoxX_, halfBoxY_;
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 +
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 <      int nMolecules = reader.getNMolecules();
160 <
161 <      fftw_complex *in, *out;
159 >      reader.readFrame(istep);
160 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
161 >      nMolecules = info_->getNGlobalMolecules();
162 >      
163 >      Mat3x3d hmat = currentSnapshot_->getHmat();
164 >      
165 > #ifdef HAVE_FFTW3_H
166        fftw_plan p;
167 <
168 <      in = fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
169 <      out = fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
170 <      p =  fftw_plan_dft_2d(nBinsX_,
119 <                            nBinsY_,
120 <                            in, out,
121 <                            FFTW_FORWARD,
122 <                            FFTW_ESTIMATE);
167 > #else
168 >      fftwnd_plan p;
169 > #endif
170 >      fftw_complex *in, *out;
171        
172 <      int i, j;  
172 >      in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
173 >      out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
174  
175 <      for(i=0; i < nBinsX_*nBinsY_; i++){
176 <        gridsample_[i].clear();
177 <        gridZ_[i].clear();
178 <      }
175 > #ifdef HAVE_FFTW3_H
176 >      p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,
177 >                           FFTW_FORWARD, FFTW_ESTIMATE);
178 > #else
179 >      p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
180 > #endif
181  
182 <      for(i=0; i < nbins; i++){
183 <        sum_bin[i].clear();
184 <        avg_bin[i].clear();
185 <        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 <      }
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 <      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 <      }
187 >      int i, j;  
188        
189        StuntDouble* sd;
190 <      reader.readFrame(istep);
152 <      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
153 <        
154 <      Mat3x3d hmat = currentSnapshot_->getHmat();
155 <
190 >      
191        lenX_ = hmat(0,0);
192        lenY_ = hmat(1,1);
193 <
193 >      
194        gridX_ = lenX_ /(nBinsX_);
195        gridY_ = lenY_ /(nBinsY_);
196 <
197 <      double halfBoxX_ = lenX_ / 2.0;      
198 <      double halfBoxY_ = lenY_ / 2.0;      
199 <        
196 >      
197 >      halfBoxX_ = lenX_ / 2.0;      
198 >      halfBoxY_ = lenY_ / 2.0;      
199 >      
200        if (evaluator_.isDynamic()) {
201          seleMan_.setSelectionSet(evaluator_.evaluate());
202        }
# Line 169 | Line 204 | namespace oopse {
204        //wrap the stuntdoubles into a cell    
205        for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
206          Vector3d pos = sd->getPos();
207 <        currentSnapshot_->wrapVector(pos);
207 >        if (usePeriodicBoundaryConditions_)
208 >          currentSnapshot_->wrapVector(pos);
209          sd->setPos(pos);
210 <      }
210 >      }
211        
212        //determine which atom belongs to which grid
213        for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
214          Vector3d pos = sd->getPos();
215          //int binNo = (pos.z() /deltaR_) - 1;
216 <        int binNoX = (pos.x() + halfBoxX_) /gridX_;
217 <        int binNoY = (pos.y() + halfBoxY_) /gridY_;
216 >        int binNoX = (int) ((pos.x() + halfBoxX_) / gridX_);
217 >        int binNoY = (int) ((pos.y() + halfBoxY_) / gridY_);
218          //std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = "  << deltaR_ << " binNo = " << binNo << "\n";
219          gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
220          gridsample_[binNoX*nBinsY_+binNoY]++;
# Line 187 | Line 223 | namespace oopse {
223        // FFT stuff depends on nx and ny, so delay allocation until we have
224        // that information
225        
190 –      for (i=0; i< nBinsX_; i++) {
191 –        for(j=0; j< nBinsY_; j++) {
192 –          newindex = i*nBinsY_ + j;
193 –          mag[newindex] = 0.0;
194 –          newmag[newindex] = 0.0;
195 –        }
196 –      }
197 –      
226        for(i = 0; i < nBinsX_; i++){
227          for(j = 0; j < nBinsY_; j++){
228            newindex = i * nBinsY_ + j;
229            if(gridsample_[newindex] > 0){
230 <            gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
230 >            gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
231            }
232          }
233        }
234 <        
234 >      
235        for (i=0; i< nBinsX_; i++) {
236          for(j=0; j< nBinsY_; j++) {
237            newindex = i*nBinsY_ + j;
# Line 248 | Line 276 | namespace oopse {
276              px = i + 1;
277              py = j;
278              newp = px*nBinsY_ + py;
279 <            if ((px < nBinsX_) && (gridsample_[newp] > 0)) {
280 <              interpsum += gridZ_[newp];
281 <              ninterp++;
282 <            }
279 >            if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
280 >              interpsum += gridZ_[newp];
281 >              ninterp++;
282 >            }
283 >        
284 >            value = interpsum / (RealType)ninterp;
285              
256 –            value = interpsum / (double)ninterp;
257 –            
286              gridZ_[newindex] = value;
287            }
288          }
# Line 268 | Line 296 | namespace oopse {
296            c_im(in[newindex]) = 0.0;
297          }
298        }
299 <      
300 <      fftw_execute(p);
299 >
300 > #ifdef HAVE_FFTW3_H
301 >      fftw_execute(p);
302 > #else
303 >      fftwnd_one(p, in, out);
304 > #endif
305        
306        for (i=0; i< nBinsX_; i++) {
307          for(j=0; j< nBinsY_; j++) {
# Line 277 | Line 309 | namespace oopse {
309            mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
310          }
311        }
312 <      
312 >
313 > #ifdef HAVE_FFTW3_H
314        fftw_destroy_plan(p);
315 + #else
316 +      fftwnd_destroy_plan(p);
317 + #endif      
318        fftw_free(out);
319        fftw_free(in);
320  
# Line 328 | Line 364 | namespace oopse {
364        //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
365        
366        maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
367 <      dfreq = maxfreq/(double)(nbins-1);
367 >      dfreq = maxfreq/(RealType)(nbins_-1);
368      
369        //printf("%lf\n", dfreq);
370        
# Line 338 | Line 374 | namespace oopse {
374        for (i=0; i< nBinsX_; i++) {
375          for(j=0; j< nBinsY_; j++) {
376            
377 <          freq_x = (double)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
378 <          freq_y = (double)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
377 >          freq_x = (RealType)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
378 >          freq_y = (RealType)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
379            
380            freq = sqrt(freq_x*freq_x + freq_y*freq_y);
381            
# Line 351 | Line 387 | namespace oopse {
387          }
388        }
389        
390 <      for ( i = 0; i < nbins; i++) {
390 >      for ( i = 0; i < nbins_; i++) {
391          if ( samples[i][istep] > 0) {
392 <          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
392 >          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nBinsX_ / (RealType)nBinsY_;
393          }
394        }    
359 –      
395      }
396 <  
397 <    for (i = 0; i < nbins; i++) {
398 <      for (j = 0; j < nFrames; j++) {
396 >
397 >    for (int i = 0; i < nbins_; i++) {
398 >      for (int j = 0; j < nFrames; j++) {
399          sum_bin[i] += bin[i][j];
400          sum_bin_sq[i] += bin[i][j] * bin[i][j];
401        }
402 <      avg_bin[i] = sum_bin[i] / (double)nFrames;
403 <      avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
404 <      for (j = 0; j < nFrames; j++) {
402 >      avg_bin[i] = sum_bin[i] / (RealType)nFrames;
403 >      avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
404 >      for (int j = 0; j < nFrames; j++) {
405          errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
406          errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
407        }
408 <      errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
409 <      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
408 >      errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
409 >      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
410      }
411  
412      printSpectrum();
413 +
414 + #else
415 +    sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
416 +    painCave.isFatal = 1;
417 +    simError();  
418 +
419 + #endif
420    }
421      
422    void Hxy::printSpectrum() {
423      std::ofstream rdfStream(outputFilename_.c_str());
424      if (rdfStream.is_open()) {
425 <      
426 <      for (int i = 0; i < nbins; i++) {
425 >
426 >      for (int i = 0; i < nbins_; ++i) {
427          if ( avg_bin[i] > 0 ){
428 <          rdfStream << i*dfreq << "\t"
428 >          rdfStream << (RealType)i * dfreq << "\t"
429                      <<pow(avg_bin[i], 2)<<"\t"
430                      <<errbin_sq[i]<<"\t"
431                      <<avg_bin[i]<<"\t"
# Line 396 | Line 438 | namespace oopse {
438        painCave.isFatal = 1;
439        simError();  
440      }
441 <    
441 >
442      rdfStream.close();
443 +
444    }
445    
446   }

Diff Legend

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