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 1044 by gezelter, Thu Sep 21 18:04:52 2006 UTC vs.
Revision 1049 by gezelter, Fri Sep 22 01:41:11 2006 UTC

# Line 71 | Line 71 | namespace oopse {
71  
72      lNumber_ = lNumber;
73      rCut_ = rCut;
74 <    mSize_ = 2*lNumber_+1;    
74 >    mSize_ = 2*lNumber_+1;  
75  
76      // Q can take values from 0 to 1
77  
78      MinQ_ = 0.0;
79 <    MaxQ_ = 1.0;
79 >    MaxQ_ = 3.0;
80      deltaQ_ = (MaxQ_ - MinQ_) / nbins;
81      Q_histogram_.resize(nbins);
82  
# Line 114 | Line 114 | namespace oopse {
114      RealType phi;
115      RealType r;
116      RealType dist;
117 <    std::map<int, ComplexType> QBar_lm;
117 >    std::map<int,ComplexType> q_lm;
118 >    std::map<int,ComplexType> QBar_lm;
119      RealType QSq_l;
120      RealType Q_l;
121      ComplexType W_l;
122      ComplexType W_l_hat;
123 <    int nBonds;
123 >    int nBonds, Nbonds;
124      SphericalHarmonic sphericalHarmonic;
125      int i, j;
126      // Make arrays for Wigner3jm
# Line 153 | Line 154 | namespace oopse {
154               rb = mol->nextRigidBody(rbIter)) {
155            rb->updateAtoms();
156          }        
157 <      }      
158 <      
157 >      }          
158 >            
159        // outer loop is over the selected StuntDoubles:
160  
161        for (sd = seleMan_.beginSelected(i); sd != NULL;
162             sd = seleMan_.nextSelected(i)) {
163  
164          myIndex = sd->getGlobalIndex();
164 –
165 –        // For this central atom, zero out nBonds and QBar_lm
166 –
165          nBonds = 0;
168 –      
166          for (int m = -lNumber_; m <= lNumber_; m++) {
167 <          QBar_lm[m] = 0.0;
167 >          q_lm[m] = 0.0;
168          }
169          
170          // inner loop is over all other atoms in the system:
# Line 197 | Line 194 | namespace oopse {
194                  
195                  for(int m = -lNumber_; m <= lNumber_; m++){
196                    sphericalHarmonic.setM(m);
197 <                  QBar_lm[m] += sphericalHarmonic.getValueAt(costheta,phi);
197 >                  q_lm[m] += sphericalHarmonic.getValueAt(costheta, phi);
198                  }
199                  nBonds++;
200                }  
201              }
202            }
203          }
204 <        
205 <        // Normalize Qbar2
206 <        for (int m = -lNumber_;m <= lNumber_; m++){
207 <          QBar_lm[m] /= nBonds;
208 <          std::cout << "m = " << m << " QBLM = " << QBar_lm[m] << "\n";
209 <        }
204 >        RealType ql = 0.0;
205 >        for(int m=-lNumber_; m<=lNumber_; m++) {          
206 >          ql += norm(QBar_lm[m]);
207 >        }        
208 >        ql *= 4.0*NumericConstant::PI/(RealType)(2*lNumber_+1);
209 >        collectHistogram(sqrt(ql)/(RealType)nBonds);
210  
211 <        // Find second order invariant Q_l
212 <
213 <        QSq_l = 0.0;
217 <        for (int m = -lNumber_; m <= lNumber_; m++){
218 <          QSq_l += norm(QBar_lm[m]);
211 >        Nbonds += nBonds;
212 >        for (int m=-lNumber_; m<=lNumber_; m++) {
213 >          QBar_lm[m] += q_lm[m];
214          }
220 –        std::cout << "qsq_l = " << QSq_l << "\n";
221 –        Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI /
222 –                   (2.0*(RealType)lNumber_ + 1.0));
223 –
224 –        // Find Third Order Invariant W_l
225 –        
226 –        W_l = 0.0;
227 –        for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
228 –          // Zero work array
229 –          for (int ii = 0; ii < mSize_; ii++){
230 –            THRCOF[i] = 0.0;
231 –          }
232 –          // Get Wigner coefficients
233 –          m1Pass = (double)m1;
234 –          Wigner3jm(&l_, &l_, &l_, &m1Pass, &m2Min, &m2Max, THRCOF, &mSize_, &error);
235 –          for (int m_index = 1; m_index < (int)(m2Max - m2Min-1.0); m_index++) {
236 –            m2 = floor(m2Min) + m_index - 1;
237 –            m3 = -m1-m2;
238 –            W_l += THRCOF[m_index]*QBar_lm[m1]*QBar_lm[m2]*QBar_lm[m3];
239 –          }
240 –        }
241 –        
242 –        W_l_hat = W_l / pow(QSq_l, 1.5);
243 –        
244 –        // accumulate histogram data for Q_l and W_l_hat:
245 –
246 –        std::cout << "Ql = " << Q_l << " Wl = " << W_l_hat << "\n";
247 –        collectHistogram(Q_l, real(W_l_hat));
248 –                
215        }
216      }
217 +
218 +    // Normalize Qbar2
219 +    for (int m = -lNumber_;m <= lNumber_; m++){
220 +      QBar_lm[m] /= Nbonds;
221 +    }
222      
223 <    writeOrderParameter();
223 >    // Find second order invariant Q_l
224      
225 +    QSq_l = 0.0;
226 +    for (int m = -lNumber_; m <= lNumber_; m++){
227 +      QSq_l += norm(QBar_lm[m]);
228 +    }
229 +    
230 +    std::cout << "qsl = " << QSq_l << "\n";
231 +    Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / (RealType)(2*lNumber_ + 1));
232 +    
233 +    // Find Third Order Invariant W_l
234 +    
235 +    W_l = 0.0;
236 +    for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
237 +      // Zero work array
238 +      for (int ii = 0; ii < mSize_; ii++){
239 +        THRCOF[ii] = 0.0;
240 +      }
241 +      // Get Wigner coefficients
242 +      m1Pass = (double)m1;
243 +      
244 +      Wigner3jm(&l_, &l_, &l_,
245 +                &m1Pass, &m2Min, &m2Max,
246 +                THRCOF, &mSize_, &error);
247 +      
248 +      for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
249 +        m2 = (int)floor(m2Min) + mmm;
250 +        m3 = -m1-m2;
251 +        W_l += THRCOF[mmm] * QBar_lm[m1] * QBar_lm[m2] * QBar_lm[m3];
252 +      }
253 +    }
254 +    
255 +    W_l_hat = W_l / pow(QSq_l, 1.5);              
256 +    
257 +    writeOrderParameter(Q_l, real(W_l_hat));    
258    }
259  
260 +  void BondOrderParameter::collectHistogram(RealType q_l) {
261  
262 <  void BondOrderParameter::collectHistogram(RealType Q_l, RealType W_l_hat) {
263 <
259 <    if (Q_l >= MinQ_ && Q_l < MaxQ_) {
260 <      int qbin = (Q_l - MinQ_) / deltaQ_;
262 >    if (q_l >= MinQ_ && q_l < MaxQ_) {
263 >      int qbin = (q_l - MinQ_) / deltaQ_;
264        Q_histogram_[qbin] += 1;
265        Qcount_++;
266 <      sumQ_ += Q_l;
267 <      sumQ2_ += Q_l * Q_l;
266 >      sumQ_ += q_l;
267 >      sumQ2_ += q_l * q_l;
268      } else {
269        sprintf( painCave.errMsg,
270 <               "Q_l value outside reasonable range\n");
270 >               "q_l value outside reasonable range\n");
271        painCave.severity = OOPSE_ERROR;
272        painCave.isFatal = 1;
273        simError();  
274      }
275  
273 –    if (W_l_hat >= MinW_ && W_l_hat < MaxW_) {
274 –      int wbin = (W_l_hat - MinW_) / deltaW_;
275 –      W_histogram_[wbin] += 1;
276 –      Wcount_++;
277 –      sumW_  += W_l_hat;
278 –      sumW2_ += W_l_hat*W_l_hat;
279 –    } else {
280 –      sprintf( painCave.errMsg,
281 –               "W_l_hat value outside reasonable range\n");
282 –      painCave.severity = OOPSE_ERROR;
283 –      painCave.isFatal = 1;
284 –      simError();  
285 –    }
276    }  
277  
288 –  void BondOrderParameter::writeOrderParameter() {
278  
279 <    std::ofstream osq((getOutputFileName() + "q").c_str());
279 >  void BondOrderParameter::writeOrderParameter(RealType ql, RealType Wlhat) {
280  
281 <    if (osq.is_open()) {
281 >    std::ofstream os(getOutputFileName().c_str());
282  
283 <      RealType qAvg = sumQ_ / (RealType) Qcount_;
295 <      RealType qStdDev = sumQ2_ / (RealType) Qcount_ - qAvg*qAvg;
283 >    if (os.is_open()) {
284        
285 <      osq << "# Bond Order Parameter Q_" << lNumber_ << "\n";
286 <      osq << "# selection: (" << selectionScript_ << ")\n";
287 <      osq << "# <Q_" << lNumber_ << ">: " << qAvg << "\n";
288 <      osq << "# std. dev.: " << qStdDev << "\n";
289 <      
285 >      os << "# Bond Order Parameters\n";
286 >      os << "# selection: (" << selectionScript_ << ")\n";
287 >      os << "# \n";
288 >      os << "# <Q_" << lNumber_ << ">: " << ql << "\n";
289 >      os << "# <W_" << lNumber_ << ">: " << Wlhat << "\n";
290        // Normalize by number of frames and write it out:
291        for (int i = 0; i < Q_histogram_.size(); ++i) {
292          RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;
293 <        osq << Qval << "\t" << Q_histogram_[i] / frameCounter_ << "\n";
293 >        os << Qval << "\t" << (RealType)Q_histogram_[i] / (RealType)Qcount_ << "\n";
294        }
307 –      
308 –      osq.close();
309 –    } else {
310 –      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
311 –              (getOutputFileName() + "q").c_str());
312 –      painCave.isFatal = 1;
313 –      simError();  
314 –    }
295  
296 <    std::ofstream osw((getOutputFileName() + "w").c_str());
296 >      os.close();
297  
318 –    if (osw.is_open()) {
319 –
320 –      RealType wAvg = sumW_ / (RealType) Wcount_;
321 –      RealType wStdDev = sumW2_ / (RealType) Wcount_ - wAvg*wAvg;
322 –      
323 –      osw << "# Bond Order Parameter W_" << lNumber_ << "\n";
324 –      osw << "# selection: (" << selectionScript_ << ")\n";
325 –      osw << "# <W_" << lNumber_ << ">: " << wAvg << "\n";
326 –      osw << "# std. dev.: " << wStdDev << "\n";
327 –      
328 –      // Normalize by number of frames and write it out:
329 –      for (int i = 0; i < W_histogram_.size(); ++i) {
330 –        RealType Wval = MinW_ + (i + 0.5) * deltaW_;
331 –        osw << Wval << "\t" << W_histogram_[i] / frameCounter_ << "\n";
332 –      }
333 –      
334 –      osw.close();
298      } else {
299        sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
300 <              (getOutputFileName() + "w").c_str());
300 >              getOutputFileName().c_str());
301        painCave.isFatal = 1;
302        simError();  
303      }

Diff Legend

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