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 957 by gezelter, Tue May 16 20:38:23 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.3 2006-05-16 20:38:23 gezelter Exp $
44 > *  @version $Id: Hxy.cpp,v 1.8 2009-11-25 20:01:59 gezelter Exp $
45   *
46   */
47  
# Line 61 | Line 58
58   #include<stdlib.h>
59   #include<math.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 73 | 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_);
# Line 83 | Line 82 | namespace oopse {
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 <    std::vector<double> mag, newmag;
119 <    double lenX_, lenY_;
120 <    double gridX_, gridY_;
121 <    double halfBoxX_, halfBoxY_;
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 >    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        reader.readFrame(istep);
# Line 120 | Line 165 | namespace oopse {
165        fftwnd_plan p;
166   #endif
167        fftw_complex *in, *out;
123 –
168        
169        in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
170        out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
# Line 131 | Line 175 | namespace oopse {
175   #else
176        p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
177   #endif
178 <      
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        int i, j;  
185        
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 –      
186        StuntDouble* sd;
187        
188        lenX_ = hmat(0,0);
# Line 158 | Line 191 | namespace oopse {
191        gridX_ = lenX_ /(nBinsX_);
192        gridY_ = lenY_ /(nBinsY_);
193        
194 <      double halfBoxX_ = lenX_ / 2.0;      
195 <      double halfBoxY_ = lenY_ / 2.0;      
194 >      halfBoxX_ = lenX_ / 2.0;      
195 >      halfBoxY_ = lenY_ / 2.0;      
196        
197        if (evaluator_.isDynamic()) {
198          seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 168 | 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 186 | Line 220 | namespace oopse {
220        // FFT stuff depends on nx and ny, so delay allocation until we have
221        // that information
222        
189 –      for (i=0; i< nBinsX_; i++) {
190 –        for(j=0; j< nBinsY_; j++) {
191 –          newindex = i*nBinsY_ + j;
192 –          mag[newindex] = 0.0;
193 –          newmag[newindex] = 0.0;
194 –        }
195 –      }
196 –      
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        }
# Line 252 | Line 278 | namespace oopse {
278                ninterp++;
279              }
280          
281 <            value = interpsum / (double)ninterp;
281 >            value = interpsum / (RealType)ninterp;
282              
283              gridZ_[newindex] = value;
284            }
# Line 267 | 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
# Line 279 | Line 306 | namespace oopse {
306            mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
307          }
308        }
309 +
310   #ifdef HAVE_FFTW3_H
311        fftw_destroy_plan(p);
312   #else
# Line 333 | 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 343 | 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 358 | Line 386 | namespace oopse {
386        
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        }    
364 –      
392      }
393 <  
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;
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 <    
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 <
390 < }
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 408 | 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