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 1046 by gezelter, Thu Sep 21 20:43:17 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  
# Line 114 | Line 114 | namespace oopse {
114      RealType phi;
115      RealType r;
116      RealType dist;
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 <        }
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;
216 <        for (int m = -lNumber_; m <= lNumber_; m++){
217 <          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          }
219 –
220 –        Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI /
221 –                   (2.0*(RealType)lNumber_ + 1.0));
222 –
223 –        // Find Third Order Invariant W_l
224 –        
225 –        W_l = 0.0;
226 –        for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
227 –          // Zero work array
228 –          for (int ii = 0; ii < mSize_; ii++){
229 –            THRCOF[ii] = 0.0;
230 –          }
231 –          // Get Wigner coefficients
232 –          m1Pass = (double)m1;
233 –
234 –          Wigner3jm(&l_, &l_, &l_,
235 –                    &m1Pass, &m2Min, &m2Max,
236 –                    THRCOF, &mSize_, &error);
237 –
238 –          for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
239 –            m2 = (int)floor(m2Min) + mmm;
240 –            m3 = -m1-m2;
241 –            W_l += THRCOF[mmm] * QBar_lm[m1] * QBar_lm[m2] * QBar_lm[m3];
242 –          }
243 –        }
244 –        
245 –        W_l_hat = W_l / pow(QSq_l, 1.5);
246 –        
247 –        // accumulate histogram data for Q_l and W_l_hat:
248 –
249 –        std::cout << "Ql = " << Q_l << " Wl = " << W_l_hat << "\n";
250 –        collectHistogram(Q_l, real(W_l_hat));
251 –                
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 <
262 <    if (Q_l >= MinQ_ && Q_l < MaxQ_) {
263 <      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  
276 –    if (W_l_hat >= MinW_ && W_l_hat < MaxW_) {
277 –      int wbin = (W_l_hat - MinW_) / deltaW_;
278 –      W_histogram_[wbin] += 1;
279 –      Wcount_++;
280 –      sumW_  += W_l_hat;
281 –      sumW2_ += W_l_hat*W_l_hat;
282 –    } else {
283 –      sprintf( painCave.errMsg,
284 –               "W_l_hat value outside reasonable range\n");
285 –      painCave.severity = OOPSE_ERROR;
286 –      painCave.isFatal = 1;
287 –      simError();  
288 –    }
276    }  
277  
291 –  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_;
298 <      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" << (RealType)Q_histogram_[i] / (RealType)Qcount_ << "\n";
293 >        os << Qval << "\t" << (RealType)Q_histogram_[i] / (RealType)Qcount_ << "\n";
294        }
310 –      
311 –      osq.close();
312 –    } else {
313 –      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
314 –              (getOutputFileName() + "q").c_str());
315 –      painCave.isFatal = 1;
316 –      simError();  
317 –    }
295  
296 <    std::ofstream osw((getOutputFileName() + "w").c_str());
296 >      os.close();
297  
321 –    if (osw.is_open()) {
322 –
323 –      RealType wAvg = sumW_ / (RealType) Wcount_;
324 –      RealType wStdDev = sumW2_ / (RealType) Wcount_ - wAvg*wAvg;
325 –      
326 –      osw << "# Bond Order Parameter W_" << lNumber_ << "\n";
327 –      osw << "# selection: (" << selectionScript_ << ")\n";
328 –      osw << "# <W_" << lNumber_ << ">: " << wAvg << "\n";
329 –      osw << "# std. dev.: " << wStdDev << "\n";
330 –      
331 –      // Normalize by number of frames and write it out:
332 –      for (int i = 0; i < W_histogram_.size(); ++i) {
333 –        RealType Wval = MinW_ + (i + 0.5) * deltaW_;
334 –        osw << Wval << "\t" << (RealType)W_histogram_[i] / (RealType)Wcount_ << "\n";
335 –      }
336 –      
337 –      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