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 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 38 | Line 29
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31   *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   *
42 *  Hxy.cpp
43 *  OOPSE-2.0
41   *
42   *  Created by Xiuquan Sun on 05/09/06.
43   *  @author  Xiuquan Sun
44 < *  @version $Id: Hxy.cpp,v 1.1 2006-05-12 21:34:43 xsun Exp $
44 > *  @version $Id: Hxy.cpp,v 1.8 2009-11-25 20:01:59 gezelter Exp $
45   *
46   */
47  
# Line 60 | Line 57
57   #include<string.h>
58   #include<stdlib.h>
59   #include<math.h>
63 #include<fftw3.h>
64 #include<mkl_lapack64.h>
60  
61 < namespace oopse {
61 > namespace OpenMD {
62    
63    Hxy::Hxy(SimInfo* info, const std::string& filename, const std::string& sele, int nbins_x, int nbins_y, int nrbins)
64      : StaticAnalyser(info, filename), selectionScript_(sele),  evaluator_(info), seleMan_(info), nBinsX_(nbins_x), nBinsY_(nbins_y), nbins_(nrbins){
# Line 75 | Line 70 | namespace oopse {
70  
71      gridsample_.resize(nBinsX_*nBinsY_);
72      gridZ_.resize(nBinsX_*nBinsY_);
73 +    mag.resize(nBinsX_*nBinsY_);
74 +    newmag.resize(nBinsX_*nBinsY_);    
75  
76 <    sum_bin.resize(nbins);
77 <    avg_bin.resize(nbins);
78 <    errbin_sum.resize(nbins);
79 <    errbin.resize(nbins);
80 <    sum_bin_sq.resize(nbins);
81 <    avg_bin_sq.resize(nbins);
82 <    errbin_sum_sq.resize(nbins);
83 <    errbin_sq.resize(nbins);
76 >    sum_bin.resize(nbins_);
77 >    avg_bin.resize(nbins_);
78 >    errbin_sum.resize(nbins_);
79 >    errbin.resize(nbins_);
80 >    sum_bin_sq.resize(nbins_);
81 >    avg_bin_sq.resize(nbins_);
82 >    errbin_sum_sq.resize(nbins_);
83 >    errbin_sq.resize(nbins_);
84  
85 +    bin.resize(nbins_);
86 +    samples.resize(nbins_);
87 +
88      setOutputName(getPrefix(filename) + ".Hxy");
89    }
90  
91 +  Hxy::~Hxy(){
92 +      gridsample_.clear();
93 +      gridZ_.clear();
94 +      sum_bin.clear();
95 +      avg_bin.clear();
96 +      errbin_sum.clear();
97 +      errbin.clear();
98 +      sum_bin_sq.clear();
99 +      avg_bin_sq.clear();
100 +      errbin_sum_sq.clear();
101 +      errbin_sq.clear();
102 +      
103 +      for(int i=0; i < bin.size(); i++)
104 +        bin[i].clear();
105 +      for(int i=0; i < samples.size(); i++)
106 +        samples[i].clear();
107 +
108 +      mag.clear();
109 +      newmag.clear();
110 +  }
111 +
112    void Hxy::process() {
113 + #if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
114      DumpReader reader(info_, dumpFilename_);    
115      int nFrames = reader.getNFrames();
116      nProcessed_ = nFrames/step_;
117 +    
118 +    for(int k=0; k < bin.size(); k++)
119 +      bin[k].resize(nFrames);
120 +    for(int k=0; k < samples.size(); k++)
121 +      samples[k].resize(nFrames);
122  
123 <    std::vector<double> mag, newmag;
124 <    double lenX_, lenY_;
125 <    double gridX_, gridY_;
126 <    double halfBoxX_, halfBoxY_;
123 >    RealType lenX_, lenY_;
124 >    RealType gridX_, gridY_;
125 >    RealType halfBoxX_, halfBoxY_;
126 >
127      int binNoX, binNoY;
128 <    double interpsum, value;
128 >    RealType interpsum, value;
129      int ninterp, px, py, newp;
130      int newx, newy, newindex, index;
131      int new_i, new_j, new_index;
132 <    double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
133 <    double maxfreqx, maxfreqy, maxfreq, dfreq;
132 >
133 >    RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
134 >    RealType maxfreqx, maxfreqy, maxfreq;
135 >
136      int whichbin;
137 +    int nMolecules;
138 +
139 +    std::fill(sum_bin.begin(), sum_bin.end(), 0.0);
140 +    std::fill(avg_bin.begin(), avg_bin.end(), 0.0);
141 +    std::fill(errbin_sum.begin(), errbin_sum.end(), 0.0);
142 +    std::fill(errbin.begin(), errbin.end(), 0.0);
143 +    std::fill(sum_bin_sq.begin(), sum_bin_sq.end(), 0.0);
144 +    std::fill(avg_bin_sq.begin(), avg_bin_sq.end(), 0.0);
145 +    std::fill(errbin_sum_sq.begin(), errbin_sum_sq.end(), 0.0);
146 +    std::fill(errbin_sq.begin(), errbin_sq.end(), 0.0);
147      
148 +    for(int i=0; i < bin.size(); i++)
149 +      std::fill(bin[i].begin(), bin[i].end(), 0.0);
150 +    
151 +    for(int i=0; i < samples.size(); i++)
152 +      std::fill(samples[i].begin(), samples[i].end(), 0);
153 +    
154      for (int istep = 0; istep < nFrames; istep += step_) {
155        
156 <      int nMolecules = reader.getNMolecules();
157 <
158 <      fftw_complex *in, *out;
156 >      reader.readFrame(istep);
157 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
158 >      nMolecules = info_->getNGlobalMolecules();
159 >      
160 >      Mat3x3d hmat = currentSnapshot_->getHmat();
161 >      
162 > #ifdef HAVE_FFTW3_H
163        fftw_plan p;
164 <
165 <      in = fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
166 <      out = fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
167 <      p =  fftw_plan_dft_2d(nBinsX_,
119 <                            nBinsY_,
120 <                            in, out,
121 <                            FFTW_FORWARD,
122 <                            FFTW_ESTIMATE);
164 > #else
165 >      fftwnd_plan p;
166 > #endif
167 >      fftw_complex *in, *out;
168        
169 <      int i, j;  
169 >      in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
170 >      out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
171  
172 <      for(i=0; i < nBinsX_*nBinsY_; i++){
173 <        gridsample_[i].clear();
174 <        gridZ_[i].clear();
175 <      }
172 > #ifdef HAVE_FFTW3_H
173 >      p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,
174 >                           FFTW_FORWARD, FFTW_ESTIMATE);
175 > #else
176 >      p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
177 > #endif
178 >
179 >      std::fill(gridsample_.begin(), gridsample_.end(), 0);
180 >      std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
181 >      std::fill(mag.begin(), mag.end(), 0.0);
182 >      std::fill(newmag.begin(), newmag.end(), 0.0);
183  
184 <      for(i=0; i < nbins; i++){
132 <        sum_bin[i].clear();
133 <        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 <      }
184 >      int i, j;  
185        
186        StuntDouble* sd;
187 <      reader.readFrame(istep);
152 <      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
153 <        
154 <      Mat3x3d hmat = currentSnapshot_->getHmat();
155 <
187 >      
188        lenX_ = hmat(0,0);
189        lenY_ = hmat(1,1);
190 <
190 >      
191        gridX_ = lenX_ /(nBinsX_);
192        gridY_ = lenY_ /(nBinsY_);
193 <
194 <      double halfBoxX_ = lenX_ / 2.0;      
195 <      double halfBoxY_ = lenY_ / 2.0;      
196 <        
193 >      
194 >      halfBoxX_ = lenX_ / 2.0;      
195 >      halfBoxY_ = lenY_ / 2.0;      
196 >      
197        if (evaluator_.isDynamic()) {
198          seleMan_.setSelectionSet(evaluator_.evaluate());
199        }
# Line 169 | Line 201 | namespace oopse {
201        //wrap the stuntdoubles into a cell    
202        for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
203          Vector3d pos = sd->getPos();
204 <        currentSnapshot_->wrapVector(pos);
204 >        if (usePeriodicBoundaryConditions_)
205 >          currentSnapshot_->wrapVector(pos);
206          sd->setPos(pos);
207 <      }
207 >      }
208        
209        //determine which atom belongs to which grid
210        for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
211          Vector3d pos = sd->getPos();
212          //int binNo = (pos.z() /deltaR_) - 1;
213 <        int binNoX = (pos.x() + halfBoxX_) /gridX_;
214 <        int binNoY = (pos.y() + halfBoxY_) /gridY_;
213 >        int binNoX = (int) ((pos.x() + halfBoxX_) / gridX_);
214 >        int binNoY = (int) ((pos.y() + halfBoxY_) / gridY_);
215          //std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = "  << deltaR_ << " binNo = " << binNo << "\n";
216          gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
217          gridsample_[binNoX*nBinsY_+binNoY]++;
# Line 187 | Line 220 | namespace oopse {
220        // FFT stuff depends on nx and ny, so delay allocation until we have
221        // that information
222        
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      
223        for(i = 0; i < nBinsX_; i++){
224          for(j = 0; j < nBinsY_; j++){
225            newindex = i * nBinsY_ + j;
226            if(gridsample_[newindex] > 0){
227 <            gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
227 >            gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
228            }
229          }
230        }
231 <        
231 >      
232        for (i=0; i< nBinsX_; i++) {
233          for(j=0; j< nBinsY_; j++) {
234            newindex = i*nBinsY_ + j;
# Line 248 | Line 273 | namespace oopse {
273              px = i + 1;
274              py = j;
275              newp = px*nBinsY_ + py;
276 <            if ((px < nBinsX_) && (gridsample_[newp] > 0)) {
277 <              interpsum += gridZ_[newp];
278 <              ninterp++;
279 <            }
280 <            
281 <            value = interpsum / (double)ninterp;
276 >            if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
277 >              interpsum += gridZ_[newp];
278 >              ninterp++;
279 >            }
280 >        
281 >            value = interpsum / (RealType)ninterp;
282              
283              gridZ_[newindex] = value;
284            }
# Line 268 | Line 293 | namespace oopse {
293            c_im(in[newindex]) = 0.0;
294          }
295        }
296 +
297 + #ifdef HAVE_FFTW3_H
298 +      fftw_execute(p);
299 + #else
300 +      fftwnd_one(p, in, out);
301 + #endif
302        
272      fftw_execute(p);
273      
303        for (i=0; i< nBinsX_; i++) {
304          for(j=0; j< nBinsY_; j++) {
305            newindex = i*nBinsY_ + j;
306            mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
307          }
308        }
309 <      
309 >
310 > #ifdef HAVE_FFTW3_H
311        fftw_destroy_plan(p);
312 + #else
313 +      fftwnd_destroy_plan(p);
314 + #endif      
315        fftw_free(out);
316        fftw_free(in);
317  
# Line 328 | Line 361 | namespace oopse {
361        //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
362        
363        maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
364 <      dfreq = maxfreq/(double)(nbins-1);
364 >      dfreq = maxfreq/(RealType)(nbins_-1);
365      
366        //printf("%lf\n", dfreq);
367        
# Line 338 | Line 371 | namespace oopse {
371        for (i=0; i< nBinsX_; i++) {
372          for(j=0; j< nBinsY_; j++) {
373            
374 <          freq_x = (double)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
375 <          freq_y = (double)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
374 >          freq_x = (RealType)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
375 >          freq_y = (RealType)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
376            
377            freq = sqrt(freq_x*freq_x + freq_y*freq_y);
378            
# Line 351 | Line 384 | namespace oopse {
384          }
385        }
386        
387 <      for ( i = 0; i < nbins; i++) {
387 >      for ( i = 0; i < nbins_; i++) {
388          if ( samples[i][istep] > 0) {
389 <          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
389 >          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nBinsX_ / (RealType)nBinsY_;
390          }
391        }    
359      
392      }
393 <  
394 <    for (i = 0; i < nbins; i++) {
395 <      for (j = 0; j < nFrames; j++) {
393 >
394 >    for (int i = 0; i < nbins_; i++) {
395 >      for (int j = 0; j < nFrames; j++) {
396          sum_bin[i] += bin[i][j];
397          sum_bin_sq[i] += bin[i][j] * bin[i][j];
398        }
399 <      avg_bin[i] = sum_bin[i] / (double)nFrames;
400 <      avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
401 <      for (j = 0; j < nFrames; j++) {
399 >      avg_bin[i] = sum_bin[i] / (RealType)nFrames;
400 >      avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
401 >      for (int j = 0; j < nFrames; j++) {
402          errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
403          errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
404        }
405 <      errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
406 <      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
405 >      errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
406 >      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
407      }
408  
409      printSpectrum();
410 +
411 + #else
412 +    sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
413 +    painCave.isFatal = 1;
414 +    simError();  
415 +
416 + #endif
417    }
418      
419    void Hxy::printSpectrum() {
420      std::ofstream rdfStream(outputFilename_.c_str());
421      if (rdfStream.is_open()) {
422 <      
423 <      for (int i = 0; i < nbins; i++) {
422 >
423 >      for (int i = 0; i < nbins_; ++i) {
424          if ( avg_bin[i] > 0 ){
425 <          rdfStream << i*dfreq << "\t"
425 >          rdfStream << (RealType)i * dfreq << "\t"
426                      <<pow(avg_bin[i], 2)<<"\t"
427                      <<errbin_sq[i]<<"\t"
428                      <<avg_bin[i]<<"\t"
# Line 396 | Line 435 | namespace oopse {
435        painCave.isFatal = 1;
436        simError();  
437      }
438 <    
438 >
439      rdfStream.close();
440 +
441    }
442    
443   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines