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 1043 by gezelter, Thu Sep 21 14:45:08 2006 UTC vs.
Revision 1047 by gezelter, Thu Sep 21 21:47:17 2006 UTC

# Line 72 | Line 72 | namespace oopse {
72      lNumber_ = lNumber;
73      rCut_ = rCut;
74      mSize_ = 2*lNumber_+1;    
75 –
76 –    // Q can take values from 0 to 1
77 –
78 –    MinQ_ = 0.0;
79 –    MaxQ_ = 1.0;
80 –    deltaQ_ = (MaxQ_ - MinQ_) / nbins;
81 –    Q_histogram_.resize(nbins);
82 –
83 –    // W_6 for icosahedral clusters is 11 / sqrt(4199) = 0.169754, so we'll
84 –    // use values for MinW_ and MaxW_ that are slightly larger than this:
85 –
86 –    MinW_ = -0.18;
87 –    MaxW_ = 0.18;
88 –    deltaW_ = (MaxW_ - MinW_) / nbins;
89 –    W_histogram_.resize(nbins);
90 –
75    }
76  
77    BondOrderParameter::~BondOrderParameter() {
94 –    Q_histogram_.clear();
95 –    W_histogram_.clear();
78    }
79  
98 –  void BondOrderParameter::initalizeHistogram() {
99 –    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
100 –    std::fill(W_histogram_.begin(), W_histogram_.end(), 0);
101 –  }
102 –
80    void BondOrderParameter::process() {
81      Molecule* mol;
82      Atom* atom;
# Line 114 | Line 91 | namespace oopse {
91      RealType phi;
92      RealType r;
93      RealType dist;
94 <    std::map<int, ComplexType> QBar_lm;
94 >    std::map<int,ComplexType> QBar_lm;
95      RealType QSq_l;
96      RealType Q_l;
97      ComplexType W_l;
# Line 155 | Line 132 | namespace oopse {
132          }        
133        }      
134        
135 +      nBonds = 0;
136 +      
137 +      for (int m = -lNumber_; m <= lNumber_; m++) {
138 +        QBar_lm[m] = 0.0;
139 +      }
140 +      
141        // outer loop is over the selected StuntDoubles:
142  
143        for (sd = seleMan_.beginSelected(i); sd != NULL;
144             sd = seleMan_.nextSelected(i)) {
145  
146          myIndex = sd->getGlobalIndex();
164 –
165 –        // For this central atom, zero out nBonds and QBar_lm
166 –
167 –        nBonds = 0;
168 –      
169 –        for (int m = -lNumber_; m <= lNumber_; m++) {
170 –          QBar_lm[m] = 0.0;
171 –        }
147          
148          // inner loop is over all other atoms in the system:
149          
# Line 197 | Line 172 | namespace oopse {
172                  
173                  for(int m = -lNumber_; m <= lNumber_; m++){
174                    sphericalHarmonic.setM(m);
175 <                  QBar_lm[m] += sphericalHarmonic.getValueAt(costheta,phi);
175 >                  QBar_lm[m] += sphericalHarmonic.getValueAt(costheta, phi);
176                  }
177                  nBonds++;
178                }  
179              }
180            }
181          }
182 <        
183 <        // Normalize Qbar2
209 <        for (int m = -lNumber_;m <= lNumber_; m++){
210 <          QBar_lm[m] /= nBonds;
211 <          std::cout << "m = " << m << " QBLM = " << QBar_lm[m] << "\n";
212 <        }
182 >      }      
183 >    }
184  
185 <        // Find second order invariant Q_l
186 <
187 <        QSq_l = 0.0;
217 <        for (int m = -lNumber_; m <= lNumber_; m++){
218 <          QSq_l += norm(QBar_lm[m]);
219 <        }
220 <        std::cout << "qsq_l = " << QSq_l << "\n";
221 <        Q_l = sqrt(QSq_l*(4.0 * NumericConstant::PI / (2.0*(RealType)lNumber_ + 1)));
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[i] = 0.0;
230 <          }
231 <          // Get Wigner coefficients
232 <          m1Pass = (double)m1;
233 <          Wigner3jm(&l_, &l_, &l_, &m1Pass, &m2Min, &m2Max, THRCOF, &mSize_, &error);
234 <          for (int m_index = 1; m_index < (int)(m2Max - m2Min-1.0); m_index++) {
235 <            m2 = floor(m2Min) + m_index - 1;
236 <            m3 = -m1-m2;
237 <            W_l += THRCOF[m_index]*QBar_lm[m1+lNumber_]*QBar_lm[m2+lNumber_]*QBar_lm[m3+lNumber_];
238 <          }
239 <        }
240 <        
241 <        W_l_hat = W_l / pow(QSq_l, 1.5);
242 <        
243 <        // accumulate histogram data for Q_l and W_l_hat:
244 <
245 <        std::cout << "Ql = " << Q_l << " Wl = " << W_l_hat << "\n";
246 <        collectHistogram(Q_l, real(W_l_hat));
247 <                
248 <      }
185 >    // Normalize Qbar2
186 >    for (int m = -lNumber_;m <= lNumber_; m++){
187 >      QBar_lm[m] /= nBonds;
188      }
189      
190 <    writeOrderParameter();
190 >    // Find second order invariant Q_l
191      
192 <  }
193 <
194 <
256 <  void BondOrderParameter::collectHistogram(RealType Q_l, RealType W_l_hat) {
257 <
258 <    if (Q_l >= MinQ_ && Q_l < MaxQ_) {
259 <      int qbin = (Q_l - MinQ_) / deltaQ_;
260 <      Q_histogram_[qbin] += 1;
261 <      Qcount_++;
262 <      sumQ_ += Q_l;
263 <      sumQ2_ += Q_l * Q_l;
264 <    } else {
265 <      sprintf( painCave.errMsg,
266 <               "Q_l value outside reasonable range\n");
267 <      painCave.severity = OOPSE_ERROR;
268 <      painCave.isFatal = 1;
269 <      simError();  
192 >    QSq_l = 0.0;
193 >    for (int m = -lNumber_; m <= lNumber_; m++){
194 >      QSq_l += norm(QBar_lm[m]);
195      }
196 <
197 <    if (W_l_hat >= MinW_ && W_l_hat < MaxW_) {
198 <      int wbin = (W_l_hat - MinW_) / deltaW_;
199 <      W_histogram_[wbin] += 1;
200 <      Wcount_++;
201 <      sumW_  += W_l_hat;
202 <      sumW2_ += W_l_hat*W_l_hat;
203 <    } else {
204 <      sprintf( painCave.errMsg,
205 <               "W_l_hat value outside reasonable range\n");
206 <      painCave.severity = OOPSE_ERROR;
207 <      painCave.isFatal = 1;
208 <      simError();  
209 <    }
285 <  }  
286 <
287 <  void BondOrderParameter::writeOrderParameter() {
288 <
289 <    std::ofstream osq((getOutputFileName() + "q").c_str());
290 <
291 <    if (osq.is_open()) {
292 <
293 <      RealType qAvg = sumQ_ / (RealType) Qcount_;
294 <      RealType qStdDev = sumQ2_ / (RealType) Qcount_ - qAvg*qAvg;
196 >    
197 >    std::cout << "qsl = " << QSq_l << "\n";
198 >    Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / (RealType)(2*lNumber_ + 1));
199 >    
200 >    // Find Third Order Invariant W_l
201 >    
202 >    W_l = 0.0;
203 >    for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
204 >      // Zero work array
205 >      for (int ii = 0; ii < mSize_; ii++){
206 >        THRCOF[ii] = 0.0;
207 >      }
208 >      // Get Wigner coefficients
209 >      m1Pass = (double)m1;
210        
211 <      osq << "# Bond Order Parameter Q_" << lNumber_ << "\n";
212 <      osq << "# selection: (" << selectionScript_ << ")\n";
213 <      osq << "# <Q_" << lNumber_ << ">: " << qAvg << "\n";
299 <      osq << "# std. dev.: " << qStdDev << "\n";
211 >      Wigner3jm(&l_, &l_, &l_,
212 >                &m1Pass, &m2Min, &m2Max,
213 >                THRCOF, &mSize_, &error);
214        
215 <      // Normalize by number of frames and write it out:
216 <      for (int i = 0; i < Q_histogram_.size(); ++i) {
217 <        RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;
218 <        osq << Qval << "\t" << Q_histogram_[i] / frameCounter_ << "\n";
215 >      for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
216 >        m2 = (int)floor(m2Min) + mmm;
217 >        m3 = -m1-m2;
218 >        W_l += THRCOF[mmm] * QBar_lm[m1] * QBar_lm[m2] * QBar_lm[m3];
219        }
306 –      
307 –      osq.close();
308 –    } else {
309 –      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
310 –              (getOutputFileName() + "q").c_str());
311 –      painCave.isFatal = 1;
312 –      simError();  
220      }
221 +    
222 +    W_l_hat = W_l / pow(QSq_l, 1.5);              
223 +    
224 +    writeOrderParameter(Q_l, real(W_l_hat));    
225 +  }
226  
315 –    std::ofstream osw((getOutputFileName() + "w").c_str());
227  
228 <    if (osw.is_open()) {
228 >  void BondOrderParameter::writeOrderParameter(RealType ql, RealType Wlhat) {
229  
230 <      RealType wAvg = sumW_ / (RealType) Wcount_;
231 <      RealType wStdDev = sumW2_ / (RealType) Wcount_ - wAvg*wAvg;
230 >    std::ofstream os(getOutputFileName().c_str());
231 >
232 >    if (os.is_open()) {
233        
234 <      osw << "# Bond Order Parameter W_" << lNumber_ << "\n";
235 <      osw << "# selection: (" << selectionScript_ << ")\n";
236 <      osw << "# <W_" << lNumber_ << ">: " << wAvg << "\n";
237 <      osw << "# std. dev.: " << wStdDev << "\n";
238 <      
239 <      // Normalize by number of frames and write it out:
240 <      for (int i = 0; i < W_histogram_.size(); ++i) {
329 <        RealType Wval = MinW_ + (i + 0.5) * deltaW_;
330 <        osw << Wval << "\t" << W_histogram_[i] / frameCounter_ << "\n";
331 <      }
332 <      
333 <      osw.close();
234 >      os << "# Bond Order Parameters\n";
235 >      os << "# selection: (" << selectionScript_ << ")\n";
236 >      os << "# \n";
237 >      os << "# <Q_" << lNumber_ << ">: " << ql << "\n";
238 >      os << "# <W_" << lNumber_ << ">: " << Wlhat << "\n";
239 >      os.close();
240 >
241      } else {
242        sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
243 <              (getOutputFileName() + "w").c_str());
243 >              getOutputFileName().c_str());
244        painCave.isFatal = 1;
245        simError();  
246      }

Diff Legend

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