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

Comparing trunk/src/applications/staticProps/BondOrderParameter.cpp (file contents):
Revision 1051 by gezelter, Mon Sep 25 22:08:33 2006 UTC vs.
Revision 1052 by gezelter, Tue Sep 26 01:30:32 2006 UTC

# Line 57 | Line 57 | namespace oopse {
57    BondOrderParameter::BondOrderParameter(SimInfo* info,
58                                           const std::string& filename,
59                                           const std::string& sele,
60 <                                         double rCut, int lMax, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
60 >                                         double rCut, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
61      
62      setOutputName(getPrefix(filename) + ".bo");
63  
# Line 68 | Line 68 | namespace oopse {
68  
69      // Set up cutoff radius and order of the Legendre Polynomial:
70  
71 –    lMax_ = lMax;
71      rCut_ = rCut;
72      nBins_ = nbins;
73      Qcount_.resize(lMax_+1);
# Line 83 | Line 82 | namespace oopse {
82      // W_6 for icosahedral clusters is 11 / sqrt(4199) = 0.169754, so we'll
83      // use values for MinW_ and MaxW_ that are slightly larger than this:
84  
85 <    MinW_ = -0.18;
86 <    MaxW_ = 0.18;
85 >    MinW_ = -0.25;
86 >    MaxW_ = 0.25;
87      deltaW_ = (MaxW_ - MinW_) / nbins;
88    }
89  
# Line 118 | Line 117 | namespace oopse {
117      RealType dist;
118      std::map<std::pair<int,int>,ComplexType> q;
119      std::vector<RealType> q_l;
120 +    std::vector<RealType> q2;
121 +    std::vector<ComplexType> w;
122 +    std::vector<ComplexType> w_hat;
123      std::map<std::pair<int,int>,ComplexType> QBar;
124      std::vector<RealType> Q2;
125      std::vector<RealType> Q;
# Line 138 | Line 140 | namespace oopse {
140      frameCounter_ = 0;
141  
142      q_l.resize(lMax_+1);
143 +    q2.resize(lMax_+1);
144 +    w.resize(lMax_+1);
145 +    w_hat.resize(lMax_+1);
146 +
147      Q2.resize(lMax_+1);
148      Q.resize(lMax_+1);
149      W.resize(lMax_+1);
# Line 223 | Line 229 | namespace oopse {
229            q_l[l] *= 4.0*NumericConstant::PI/(RealType)(2*l + 1);
230            q_l[l] = sqrt(q_l[l])/(RealType)nBonds;
231          }
232 <        collectHistogram(q_l);
232 >
233 >        // Find second order invariant Q_l
234          
235 +        for (int l = 0; l <= lMax_; l++) {
236 +          q2[l] = 0.0;
237 +          for (int m = -l; m <= l; m++){
238 +            q2[l] += norm(q[std::make_pair(l,m)]);
239 +          }
240 +          q_l[l] = sqrt(q2[l] * 4.0 * NumericConstant::PI /
241 +                        (RealType)(2*l + 1))/(RealType)nBonds;
242 +        }
243 +
244 +        // Find Third Order Invariant W_l
245 +    
246 +        for (int l = 0; l <= lMax_; l++) {
247 +          w[l] = 0.0;
248 +          lPass = (double)l;
249 +          for (int m1 = -l; m1 <= l; m1++) {
250 +            // Zero work array
251 +            for (int ii = 0; ii < 2*l + 1; ii++){
252 +              THRCOF[ii] = 0.0;
253 +            }
254 +            // Get Wigner coefficients
255 +            m1Pass = (double)m1;
256 +            
257 +            Wigner3jm(&lPass, &lPass, &lPass,
258 +                      &m1Pass, &m2Min, &m2Max,
259 +                      THRCOF, &mSize, &error);
260 +            
261 +            for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
262 +              m2 = (int)floor(m2Min) + mmm;
263 +              m3 = -m1-m2;
264 +              w[l] += THRCOF[mmm] *
265 +                q[std::make_pair(l,m1)] *
266 +                q[std::make_pair(l,m2)] *
267 +                q[std::make_pair(l,m3)];
268 +            }
269 +          }
270 +          
271 +          w_hat[l] = w[l] / pow(q2[l], 1.5);
272 +        }
273 +
274 +        collectHistogram(q_l, w_hat);
275 +        
276          Nbonds += nBonds;
277          for (int l = 0; l <= lMax_;  l++) {
278            for (int m = -l; m <= l; m++) {
# Line 251 | Line 299 | namespace oopse {
299        Q[l] = sqrt(Q2[l] * 4.0 * NumericConstant::PI / (RealType)(2*l + 1));
300      }
301      
254 –
255 –    
302      // Find Third Order Invariant W_l
303      
304      for (int l = 0; l <= lMax_; l++) {
# Line 286 | Line 332 | namespace oopse {
332      writeOrderParameter(Q, W_hat);    
333    }
334  
335 <  void BondOrderParameter::collectHistogram(std::vector<RealType> q) {
335 >  void BondOrderParameter::collectHistogram(std::vector<RealType> q,
336 >                                            std::vector<ComplexType> what) {
337  
338      for (int l = 0; l <= lMax_; l++) {
339        if (q[l] >= MinQ_ && q[l] < MaxQ_) {
# Line 302 | Line 349 | namespace oopse {
349        }
350      }
351  
352 +    for (int l = 0; l <= lMax_; l++) {
353 +      if (real(what[l]) >= MinW_ && real(what[l]) < MaxW_) {
354 +        int wbin = (real(what[l]) - MinW_) / deltaW_;
355 +        W_histogram_[std::make_pair(wbin,l)] += 1;
356 +        Wcount_[l]++;      
357 +      } else {
358 +        sprintf( painCave.errMsg,
359 +                 "Re[w_hat] value outside reasonable range\n");
360 +        painCave.severity = OOPSE_ERROR;
361 +        painCave.isFatal = 1;
362 +        simError();  
363 +      }
364 +    }
365 +
366    }  
367  
368  
369 <  void BondOrderParameter::writeOrderParameter(std::vector<RealType> Q, std::vector<ComplexType> What) {
369 >  void BondOrderParameter::writeOrderParameter(std::vector<RealType> Q,
370 >                                               std::vector<ComplexType> What) {
371      
372 <    std::ofstream os(getOutputFileName().c_str());
373 <    
374 <    if (os.is_open()) {
372 >    std::ofstream osq((getOutputFileName() + "q").c_str());
373 >
374 >    if (osq.is_open()) {
375        
376 <      os << "# Bond Order Parameters\n";
377 <      os << "# selection: (" << selectionScript_ << ")\n";
376 >      osq << "# Bond Order Parameters\n";
377 >      osq << "# selection: (" << selectionScript_ << ")\n";
378 >      osq << "# \n";
379        for (int l = 0; l <= lMax_; l++) {
380 <        os << "# \n";
318 <        os << "# <Q_" << l << ">: " << Q[l] << "\n";
319 <        os << "# <W_" << l << ">: " << real(What[l]) << "\n";
380 >        osq << "# <Q_" << l << ">: " << Q[l] << "\n";
381        }
382        // Normalize by number of frames and write it out:
383        for (int i = 0; i < nBins_; ++i) {
384          RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;              
385 <        os << Qval;
385 >        osq << Qval;
386          for (int l = 0; l <= lMax_; l++) {
387 <          os << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)] / (RealType)Qcount_[l];
387 >          osq << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)] / (RealType)Qcount_[l];
388          }
389 <        os << "\n";
389 >        osq << "\n";
390        }
391  
392 <      os.close();
392 >      osq.close();
393  
394      } else {
395        sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
396 <              getOutputFileName().c_str());
396 >              (getOutputFileName() + "q").c_str());
397        painCave.isFatal = 1;
398        simError();  
399      }
400 +
401 +    std::ofstream osw((getOutputFileName() + "w").c_str());
402 +
403 +    if (osw.is_open()) {
404 +      osw << "# Bond Order Parameters\n";
405 +      osw << "# selection: (" << selectionScript_ << ")\n";
406 +      osw << "# \n";
407 +      for (int l = 0; l <= lMax_; l++) {
408 +        osw << "# <W_" << l << ">: " << real(What[l]) << "\n";
409 +      }
410 +      // Normalize by number of frames and write it out:
411 +      for (int i = 0; i < nBins_; ++i) {
412 +        RealType Wval = MinW_ + (i + 0.5) * deltaW_;              
413 +        osw << Wval;
414 +        for (int l = 0; l <= lMax_; l++) {
415 +          osw << "\t" << (RealType)W_histogram_[std::make_pair(i,l)] / (RealType)Wcount_[l];
416 +        }
417 +        osw << "\n";
418 +      }
419 +
420 +      osw.close();
421 +    } else {
422 +      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
423 +              (getOutputFileName() + "w").c_str());
424 +      painCave.isFatal = 1;
425 +      simError();  
426 +    }
427 +      
428    }
429   }

Diff Legend

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