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 1053 by gezelter, Tue Sep 26 14:33:56 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;
89 –  }
88  
89 +    // Make arrays for Wigner3jm
90 +    double* THRCOF = new double[2*lMax_+1];
91 +    // Variables for Wigner routine
92 +    double lPass, m1Pass, m2m, m2M;
93 +    int error, mSize;
94 +    mSize = 2*lMax_+1;
95 +
96 +    for (int l = 0; l <= lMax_; l++) {
97 +      lPass = (double)l;
98 +      for (int m1 = -l; m1 <= l; m1++) {
99 +        m1Pass = (double)m1;
100 +
101 +        std::pair<int,int> lm = std::make_pair(l, m1);
102 +        
103 +        // Zero work array
104 +        for (int ii = 0; ii < 2*l + 1; ii++){
105 +          THRCOF[ii] = 0.0;
106 +        }
107 +            
108 +        // Get Wigner coefficients
109 +        Wigner3jm(&lPass, &lPass, &lPass,
110 +                  &m1Pass, &m2m, &m2M,
111 +                  THRCOF, &mSize, &error);
112 +        
113 +        m2Min[lm] = (int)floor(m2m);
114 +        m2Max[lm] = (int)floor(m2M);
115 +        
116 +        for (int mmm = 0; mmm < (int)(m2M - m2m); mmm++) {
117 +          w3j[lm].push_back(THRCOF[mmm]);
118 +        }
119 +      }
120 +    }
121 +  }
122 +  
123    BondOrderParameter::~BondOrderParameter() {
124      Q_histogram_.clear();
125      W_histogram_.clear();
# Line 118 | Line 150 | namespace oopse {
150      RealType dist;
151      std::map<std::pair<int,int>,ComplexType> q;
152      std::vector<RealType> q_l;
153 +    std::vector<RealType> q2;
154 +    std::vector<ComplexType> w;
155 +    std::vector<ComplexType> w_hat;
156      std::map<std::pair<int,int>,ComplexType> QBar;
157      std::vector<RealType> Q2;
158      std::vector<RealType> Q;
# Line 126 | Line 161 | namespace oopse {
161      int nBonds, Nbonds;
162      SphericalHarmonic sphericalHarmonic;
163      int i, j;
129 –    // Make arrays for Wigner3jm
130 –    double* THRCOF = new double[2*lMax_+1];
131 –    // Variables for Wigner routine
132 –    double lPass, m1Pass, m2Min, m2Max;
133 –    int error, m1, m2, m3, mSize;
134 –    mSize = 2*lMax_+1;
164  
165      DumpReader reader(info_, dumpFilename_);    
166      int nFrames = reader.getNFrames();
167      frameCounter_ = 0;
168  
169      q_l.resize(lMax_+1);
170 +    q2.resize(lMax_+1);
171 +    w.resize(lMax_+1);
172 +    w_hat.resize(lMax_+1);
173 +
174      Q2.resize(lMax_+1);
175      Q.resize(lMax_+1);
176      W.resize(lMax_+1);
# Line 223 | Line 256 | namespace oopse {
256            q_l[l] *= 4.0*NumericConstant::PI/(RealType)(2*l + 1);
257            q_l[l] = sqrt(q_l[l])/(RealType)nBonds;
258          }
259 <        collectHistogram(q_l);
259 >
260 >        // Find second order invariant Q_l
261          
262 +        for (int l = 0; l <= lMax_; l++) {
263 +          q2[l] = 0.0;
264 +          for (int m = -l; m <= l; m++){
265 +            q2[l] += norm(q[std::make_pair(l,m)]);
266 +          }
267 +          q_l[l] = sqrt(q2[l] * 4.0 * NumericConstant::PI /
268 +                        (RealType)(2*l + 1))/(RealType)nBonds;
269 +        }
270 +
271 +        // Find Third Order Invariant W_l
272 +    
273 +        for (int l = 0; l <= lMax_; l++) {
274 +          w[l] = 0.0;
275 +          for (int m1 = -l; m1 <= l; m1++) {
276 +            std::pair<int,int> lm = std::make_pair(l, m1);
277 +            for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) {
278 +              int m2 = m2Min[lm] + mmm;
279 +              int m3 = -m1-m2;
280 +              w[l] += w3j[lm][mmm] * q[lm] *
281 +                q[std::make_pair(l,m2)] *  q[std::make_pair(l,m3)];
282 +            }
283 +          }
284 +          
285 +          w_hat[l] = w[l] / pow(q2[l], 1.5);
286 +        }
287 +
288 +        collectHistogram(q_l, w_hat);
289 +        
290          Nbonds += nBonds;
291          for (int l = 0; l <= lMax_;  l++) {
292            for (int m = -l; m <= l; m++) {
# Line 251 | Line 313 | namespace oopse {
313        Q[l] = sqrt(Q2[l] * 4.0 * NumericConstant::PI / (RealType)(2*l + 1));
314      }
315      
254 –
255 –    
316      // Find Third Order Invariant W_l
317      
318      for (int l = 0; l <= lMax_; l++) {
319        W[l] = 0.0;
260 –      lPass = (double)l;
320        for (int m1 = -l; m1 <= l; m1++) {
321 <        // Zero work array
322 <        for (int ii = 0; ii < 2*l + 1; ii++){
323 <          THRCOF[ii] = 0.0;
324 <        }
325 <        // Get Wigner coefficients
326 <        m1Pass = (double)m1;
268 <        
269 <        Wigner3jm(&lPass, &lPass, &lPass,
270 <                  &m1Pass, &m2Min, &m2Max,
271 <                  THRCOF, &mSize, &error);
272 <        
273 <        for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
274 <          m2 = (int)floor(m2Min) + mmm;
275 <          m3 = -m1-m2;
276 <          W[l] += THRCOF[mmm] *
277 <            QBar[std::make_pair(l,m1)] *
278 <            QBar[std::make_pair(l,m2)] *
279 <            QBar[std::make_pair(l,m3)];
321 >        std::pair<int,int> lm = std::make_pair(l, m1);
322 >        for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) {
323 >          int m2 = m2Min[lm] + mmm;
324 >          int m3 = -m1-m2;
325 >          W[l] += w3j[lm][mmm] * QBar[lm] *
326 >            QBar[std::make_pair(l,m2)] * QBar[std::make_pair(l,m3)];
327          }
328        }
329        
# Line 286 | Line 333 | namespace oopse {
333      writeOrderParameter(Q, W_hat);    
334    }
335  
336 <  void BondOrderParameter::collectHistogram(std::vector<RealType> q) {
336 >  void BondOrderParameter::collectHistogram(std::vector<RealType> q,
337 >                                            std::vector<ComplexType> what) {
338  
339      for (int l = 0; l <= lMax_; l++) {
340        if (q[l] >= MinQ_ && q[l] < MaxQ_) {
# Line 302 | Line 350 | namespace oopse {
350        }
351      }
352  
353 +    for (int l = 0; l <= lMax_; l++) {
354 +      if (real(what[l]) >= MinW_ && real(what[l]) < MaxW_) {
355 +        int wbin = (real(what[l]) - MinW_) / deltaW_;
356 +        W_histogram_[std::make_pair(wbin,l)] += 1;
357 +        Wcount_[l]++;      
358 +      } else {
359 +        sprintf( painCave.errMsg,
360 +                 "Re[w_hat] value outside reasonable range\n");
361 +        painCave.severity = OOPSE_ERROR;
362 +        painCave.isFatal = 1;
363 +        simError();  
364 +      }
365 +    }
366 +
367    }  
368  
369  
370 <  void BondOrderParameter::writeOrderParameter(std::vector<RealType> Q, std::vector<ComplexType> What) {
370 >  void BondOrderParameter::writeOrderParameter(std::vector<RealType> Q,
371 >                                               std::vector<ComplexType> What) {
372      
373 <    std::ofstream os(getOutputFileName().c_str());
374 <    
375 <    if (os.is_open()) {
373 >    std::ofstream osq((getOutputFileName() + "q").c_str());
374 >
375 >    if (osq.is_open()) {
376        
377 <      os << "# Bond Order Parameters\n";
378 <      os << "# selection: (" << selectionScript_ << ")\n";
377 >      osq << "# Bond Order Parameters\n";
378 >      osq << "# selection: (" << selectionScript_ << ")\n";
379 >      osq << "# \n";
380        for (int l = 0; l <= lMax_; l++) {
381 <        os << "# \n";
318 <        os << "# <Q_" << l << ">: " << Q[l] << "\n";
319 <        os << "# <W_" << l << ">: " << real(What[l]) << "\n";
381 >        osq << "# <Q_" << l << ">: " << Q[l] << "\n";
382        }
383        // Normalize by number of frames and write it out:
384        for (int i = 0; i < nBins_; ++i) {
385          RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;              
386 <        os << Qval;
386 >        osq << Qval;
387          for (int l = 0; l <= lMax_; l++) {
388 <          os << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)] / (RealType)Qcount_[l];
388 >          osq << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)] / (RealType)Qcount_[l];
389          }
390 <        os << "\n";
390 >        osq << "\n";
391        }
392  
393 <      os.close();
393 >      osq.close();
394  
395      } else {
396        sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
397 <              getOutputFileName().c_str());
397 >              (getOutputFileName() + "q").c_str());
398        painCave.isFatal = 1;
399        simError();  
400      }
401 +
402 +    std::ofstream osw((getOutputFileName() + "w").c_str());
403 +
404 +    if (osw.is_open()) {
405 +      osw << "# Bond Order Parameters\n";
406 +      osw << "# selection: (" << selectionScript_ << ")\n";
407 +      osw << "# \n";
408 +      for (int l = 0; l <= lMax_; l++) {
409 +        osw << "# <W_" << l << ">: " << real(What[l]) << "\n";
410 +      }
411 +      // Normalize by number of frames and write it out:
412 +      for (int i = 0; i < nBins_; ++i) {
413 +        RealType Wval = MinW_ + (i + 0.5) * deltaW_;              
414 +        osw << Wval;
415 +        for (int l = 0; l <= lMax_; l++) {
416 +          osw << "\t" << (RealType)W_histogram_[std::make_pair(i,l)] / (RealType)Wcount_[l];
417 +        }
418 +        osw << "\n";
419 +      }
420 +
421 +      osw.close();
422 +    } else {
423 +      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
424 +              (getOutputFileName() + "w").c_str());
425 +      painCave.isFatal = 1;
426 +      simError();  
427 +    }
428 +      
429    }
430   }

Diff Legend

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