ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/TetrahedralityParamXYZ.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/TetrahedralityParamXYZ.cpp (file contents):
Revision 2016 by plouden, Thu Aug 14 19:04:30 2014 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 56 | Line 56 | namespace OpenMD {
56                                                   const std::string& filename,
57                                                   const std::string& sele1,
58                                                   const std::string& sele2,
59 <                                                 RealType rCut, RealType voxelSize,
59 >                                                 RealType rCut,
60 >                                                 RealType voxelSize,
61                                                   RealType gaussWidth)
62 <    : StaticAnalyser(info, filename), selectionScript1_(sele1),
63 <      evaluator1_(info), seleMan1_(info), selectionScript2_(sele2),
64 <      evaluator2_(info), seleMan2_(info), rCut_(rCut), voxelSize_(voxelSize),
65 <      gaussWidth_(gaussWidth) {
62 >    : StaticAnalyser(info, filename),
63 >      selectionScript1_(sele1), selectionScript2_(sele2),
64 >      seleMan1_(info),  seleMan2_(info), evaluator1_(info), evaluator2_(info),
65 >      rCut_(rCut), voxelSize_(voxelSize), gaussWidth_(gaussWidth) {
66      
67      evaluator1_.loadScriptString(sele1);
68      if (!evaluator1_.isDynamic()) {
# Line 114 | Line 115 | namespace OpenMD {
115      RealType cospsi;
116      RealType Qk;
117      std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
118 <    std::vector<std::pair<Vector3d, RealType> > qvals;
119 <    std::vector<std::pair<Vector3d, RealType> >::iterator qiter;
118 >    //std::vector<std::pair<Vector3d, RealType> > qvals;
119 >    //std::vector<std::pair<Vector3d, RealType> >::iterator qiter;
120      int isd1;
121      int isd2;
122  
123 +
124 +    int kMax = int(5.0 * gaussWidth_ / voxelSize_);
125 +    int kSqLim = kMax*kMax;
126 +    cerr << "gw = " << gaussWidth_ << " vS = " << voxelSize_ << " kMax = "
127 +         << kMax << " kSqLim = " << kSqLim << "\n";
128 +    
129      DumpReader reader(info_, dumpFilename_);    
130      int nFrames = reader.getNFrames();
131  
# Line 146 | Line 153 | namespace OpenMD {
153          }
154        }
155        
156 <      qvals.clear();
156 >      //qvals.clear();
157  
158        // outer loop is over the selected StuntDoubles:
159        for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
# Line 217 | Line 224 | namespace OpenMD {
224          if (nang > 0) {
225            if (usePeriodicBoundaryConditions_)
226              currentSnapshot_->wrapVector(rk);
227 <          qvals.push_back(std::make_pair(rk, Qk));
228 <        }  
229 <      }
227 >          //qvals.push_back(std::make_pair(rk, Qk));
228 >          
229 >          Vector3d pos = rk + halfBox;
230  
231 <      for (int i = 0; i < nBins_(0); ++i) {
232 <        for(int j = 0; j < nBins_(1); ++j) {
233 <          for(int k = 0; k < nBins_(2); ++k) {
234 <            Vector3d pos = Vector3d(i, j, k) * voxelSize_ - halfBox;
235 <            for(qiter = qvals.begin(); qiter != qvals.end(); ++qiter) {
236 <              Vector3d d = pos - (*qiter).first;
237 <              currentSnapshot_->wrapVector(d);
238 <              RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3);
239 <              RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2);
240 <              RealType weight = exp(exponent) / denom;
241 <              count_[i][j][k] += weight;
242 <              hist_[i][j][k] += weight * (*qiter).second;
231 >
232 >          Vector3i whichVoxel(int(pos[0] / voxelSize_),
233 >                              int(pos[1] / voxelSize_),
234 >                              int(pos[2] / voxelSize_));
235 >          
236 >          for (int l = -kMax; l <= kMax; l++) {
237 >            for (int m = -kMax; m <= kMax; m++) {
238 >              for (int n = -kMax; n <= kMax; n++) {
239 >                int kk = l*l + m*m + n*n;
240 >                if(kk <= kSqLim) {
241 >
242 >                  int ll = (whichVoxel[0] + l) % nBins_(0);
243 >                  ll = ll < 0 ? nBins_(0) + ll : ll;
244 >                  int mm = (whichVoxel[1] + m) % nBins_(1);
245 >                  mm = mm < 0 ? nBins_(1) + mm : mm;
246 >                  int nn = (whichVoxel[2] + n) % nBins_(2);
247 >                  nn = nn < 0 ? nBins_(2) + nn : nn;
248 >
249 >                  Vector3d bPos = Vector3d(ll,mm,nn) * voxelSize_ - halfBox;
250 >                  Vector3d d = bPos - rk;
251 >                  currentSnapshot_->wrapVector(d);
252 >                  RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3);
253 >                  RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2);
254 >                  RealType weight = exp(exponent) / denom;
255 >                  count_[ll][mm][nn] += weight;
256 >                  hist_[ll][mm][nn] += weight * Qk;
257 >                }
258 >              }
259              }
260            }
261          }
262        }
263 +
264 +      // for (int i = 0; i < nBins_(0); ++i) {
265 +      //   for(int j = 0; j < nBins_(1); ++j) {
266 +      //     for(int k = 0; k < nBins_(2); ++k) {
267 +      //       Vector3d pos = Vector3d(i, j, k) * voxelSize_ - halfBox;
268 +      //       for(qiter = qvals.begin(); qiter != qvals.end(); ++qiter) {
269 +      //         Vector3d d = pos - (*qiter).first;
270 +      //         currentSnapshot_->wrapVector(d);
271 +      //         RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3);
272 +      //         RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2);
273 +      //         RealType weight = exp(exponent) / denom;
274 +      //         count_[i][j][k] += weight;
275 +      //         hist_[i][j][k] += weight * (*qiter).second;
276 +      //       }
277 +      //     }
278 +      //   }
279 +      // }
280      }
281      writeQxyz();
282    }
# Line 258 | Line 298 | namespace OpenMD {
298      if (qXYZstream.is_open()) {
299        qXYZstream <<  "# AmiraMesh ASCII 1.0\n\n";
300        qXYZstream <<  "# Dimensions in x-, y-, and z-direction\n";
301 <      qXYZstream <<  "  define Lattice " << hist_.size() << " " << hist_[0].size() << " " << hist_[0][0].size() << "\n";
301 >      qXYZstream <<  "  define Lattice " << hist_.size() << " "
302 >                 << hist_[0].size() << " " << hist_[0][0].size() << "\n";
303        
304        qXYZstream <<  "Parameters {\n";
305        qXYZstream <<  "    CoordType \"uniform\",\n";
# Line 272 | Line 313 | namespace OpenMD {
313        
314        qXYZstream <<  "@1\n";
315  
316 <      int xsize = hist_.size();
317 <      int ysize = hist_[0].size();
318 <      int zsize = hist_[0][0].size();
278 <      
279 <      for (unsigned int k = 0; k < zsize; ++k) {
280 <        for(unsigned int j = 0; j < ysize; ++j) {
281 <          for(unsigned int i = 0; i < xsize; ++i) {
316 >      for (std::size_t k = 0; k < hist_[0][0].size(); ++k) {
317 >        for(std::size_t j = 0; j < hist_[0].size(); ++j) {
318 >          for(std::size_t i = 0; i < hist_.size(); ++i) {
319              qXYZstream << hist_[i][j][k] << " ";
320  
321              //qXYZstream.write(reinterpret_cast<char *>( &hist_[i][j][k] ),

Diff Legend

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