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 1042 by gezelter, Wed Sep 20 22:16:23 2006 UTC vs.
Revision 1048 by gezelter, Fri Sep 22 01:36:27 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 104 | Line 104 | namespace oopse {
104      Molecule* mol;
105      Atom* atom;
106      RigidBody* rb;
107 +    int myIndex;
108      SimInfo::MoleculeIterator mi;
109      Molecule::RigidBodyIterator rbIter;
110      Molecule::AtomIterator ai;
111      StuntDouble* sd;
112 +    Vector3d vec;
113      RealType costheta;
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 <    int nBonds;
121 >    ComplexType W_l;
122 >    ComplexType W_l_hat;
123 >    int nBonds, Nbonds;
124      SphericalHarmonic sphericalHarmonic;
125      int i, j;
126 <  
126 >    // Make arrays for Wigner3jm
127 >    double* THRCOF = new double[mSize_];
128 >    // Variables for Wigner routine
129 >    double l_ = (double)lNumber_;
130 >    double m1Pass, m2Min, m2Max;
131 >    int error, m1, m2, m3;
132 >
133      // Set the l for the spherical harmonic, it doesn't change
134      sphericalHarmonic.setL(lNumber_);
135  
125 –
136      DumpReader reader(info_, dumpFilename_);    
137      int nFrames = reader.getNFrames();
138      frameCounter_ = 0;
# Line 144 | 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 <        // For this central atom, zero out nBonds and QBar_lm
155 <
164 >        myIndex = sd->getGlobalIndex();
165          nBonds = 0;
157 –      
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 166 | Line 174 | namespace oopse {
174            for (atom = mol->beginAtom(ai); atom != NULL;
175                 atom = mol->nextAtom(ai)) {
176  
177 +            if (atom->getGlobalIndex() != myIndex) {
178  
179 <            Vector3d vec = sd->getPos() - atom->getPos();      
180 <            currentSnapshot_->wrapVector(vec);
172 <            
173 <            // Calculate "bonds" and build Q_lm(r) where
174 <            //      Q_lm = Y_lm(theta(r),phi(r))                
175 <            // The spherical harmonics are wrt any arbitrary coordinate
176 <            // system, we choose standard spherical coordinates
177 <            
178 <            r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
179 <            
180 <            // Check to see if neighbor is in bond cutoff
181 <            
182 <            if (r < rCut_) {            
183 <              costheta = vec.z() / r;
184 <              phi = atan2(vec.y(), vec.x());
179 >              vec = sd->getPos() - atom->getPos();      
180 >              currentSnapshot_->wrapVector(vec);
181                
182 <              for(int m = -lNumber_; m <= lNumber_; m++){
183 <                sphericalHarmonic.setM(m);
184 <                QBar_lm[m] += sphericalHarmonic.getValueAt(costheta,phi);
185 <              }
186 <              nBonds++;
187 <            }  
182 >              // Calculate "bonds" and build Q_lm(r) where
183 >              //      Q_lm = Y_lm(theta(r),phi(r))                
184 >              // The spherical harmonics are wrt any arbitrary coordinate
185 >              // system, we choose standard spherical coordinates
186 >              
187 >              r = vec.length();
188 >              
189 >              // Check to see if neighbor is in bond cutoff
190 >              
191 >              if (r < rCut_) {
192 >                costheta = vec.z() / r;
193 >                phi = atan2(vec.y(), vec.x());
194 >                
195 >                for(int m = -lNumber_; m <= lNumber_; m++){
196 >                  sphericalHarmonic.setM(m);
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;
203 <        for (int m = -lNumber_; m <= lNumber_; m++){
204 <          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          }
215 <        Q_l = sqrt(QSq_l*(4.0 * NumericConstant::PI / (2.0*(RealType)lNumber_ + 1)));
216 <    
208 <        // Find Third Order Invariant W_l
215 >      }
216 >    }
217  
218 <        // Make arrays for Wigner3jm
219 <        double* THRCOF = new double[mSize_];
220 <        // Variables for Wigner routine
221 <        double l_ = (double)lNumber_;
222 <        double m1Pass, m2Min, m2Max;
223 <        int error, m1, m2, m3;
224 <
225 <        ComplexType W_l;
226 <        ComplexType W_l_hat;
227 <        W_l = 0.0;
228 <        for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
229 <          // Zero work array
230 <          for (int ii = 0; ii < mSize_; ii++){
231 <            THRCOF[i] = 0.0;
232 <          }
233 <          // Get Wigner coefficients
234 <          m1Pass = (double)m1;
235 <          Wigner3jm(&l_, &l_, &l_, &m1Pass, &m2Min, &m2Max, THRCOF, &mSize_, &error);
236 <          for (int m_index = 1; i < (int)(m2Max - m2Min-1.0); m_index++) {
237 <            m2 = floor(m2Min) + m_index - 1;
238 <            m3 = -m1-m2;
239 <            W_l += THRCOF[m_index]*QBar_lm[m1+lNumber_]*QBar_lm[m2+lNumber_]*QBar_lm[m3+lNumber_];
232 <          }
233 <        }
234 <
235 <        W_l_hat = W_l / pow(QSq_l, 1.5);
236 <
237 <        // accumulate histogram data for Q_l and W_l_hat:
238 <
239 <        collectHistogram(Q_l, real(W_l_hat));
240 <                
218 >    // Normalize Qbar2
219 >    for (int m = -lNumber_;m <= lNumber_; m++){
220 >      QBar_lm[m] /= Nbonds;
221 >    }
222 >    
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 <    writeOrderParameter();
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 <
251 <    if (Q_l >= MinQ_ && Q_l < MaxQ_) {
252 <      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  
265 –    if (W_l_hat >= MinW_ && W_l_hat < MaxW_) {
266 –      int wbin = (W_l_hat - MinW_) / deltaW_;
267 –      W_histogram_[wbin] += 1;
268 –      Wcount_++;
269 –      sumW_  += W_l_hat;
270 –      sumW2_ += W_l_hat*W_l_hat;
271 –    } else {
272 –      sprintf( painCave.errMsg,
273 –               "W_l_hat value outside reasonable range\n");
274 –      painCave.severity = OOPSE_ERROR;
275 –      painCave.isFatal = 1;
276 –      simError();  
277 –    }
276    }  
277  
280 –  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_;
287 <      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 >        osq << Qval << "\t" << (RealType)Q_histogram_[i] / (RealType)Qcount_ << "\n";
294        }
299 –      
300 –      osq.close();
301 –    } else {
302 –      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
303 –              (getOutputFileName() + "q").c_str());
304 –      painCave.isFatal = 1;
305 –      simError();  
306 –    }
295  
296 <    std::ofstream osw((getOutputFileName() + "w").c_str());
296 >      os.close();
297  
310 –    if (osw.is_open()) {
311 –
312 –      RealType wAvg = sumW_ / (RealType) Wcount_;
313 –      RealType wStdDev = sumW2_ / (RealType) Wcount_ - wAvg*wAvg;
314 –      
315 –      osw << "# Bond Order Parameter W_" << lNumber_ << "\n";
316 –      osw << "# selection: (" << selectionScript_ << ")\n";
317 –      osw << "# <W_" << lNumber_ << ">: " << wAvg << "\n";
318 –      osw << "# std. dev.: " << wStdDev << "\n";
319 –      
320 –      // Normalize by number of frames and write it out:
321 –      for (int i = 0; i < W_histogram_.size(); ++i) {
322 –        RealType Wval = MinW_ + (i + 0.5) * deltaW_;
323 –        osw << Wval << "\t" << W_histogram_[i] / frameCounter_ << "\n";
324 –      }
325 –      
326 –      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