--- trunk/src/applications/staticProps/BondOrderParameter.cpp 2006/09/21 18:04:52 1044 +++ trunk/src/applications/staticProps/BondOrderParameter.cpp 2006/09/22 01:36:27 1048 @@ -71,12 +71,12 @@ namespace oopse { lNumber_ = lNumber; rCut_ = rCut; - mSize_ = 2*lNumber_+1; + mSize_ = 2*lNumber_+1; // Q can take values from 0 to 1 MinQ_ = 0.0; - MaxQ_ = 1.0; + MaxQ_ = 3.0; deltaQ_ = (MaxQ_ - MinQ_) / nbins; Q_histogram_.resize(nbins); @@ -114,12 +114,13 @@ namespace oopse { RealType phi; RealType r; RealType dist; - std::map QBar_lm; + std::map q_lm; + std::map QBar_lm; RealType QSq_l; RealType Q_l; ComplexType W_l; ComplexType W_l_hat; - int nBonds; + int nBonds, Nbonds; SphericalHarmonic sphericalHarmonic; int i, j; // Make arrays for Wigner3jm @@ -153,21 +154,17 @@ namespace oopse { rb = mol->nextRigidBody(rbIter)) { rb->updateAtoms(); } - } - + } + // outer loop is over the selected StuntDoubles: for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) { myIndex = sd->getGlobalIndex(); - - // For this central atom, zero out nBonds and QBar_lm - nBonds = 0; - for (int m = -lNumber_; m <= lNumber_; m++) { - QBar_lm[m] = 0.0; + q_lm[m] = 0.0; } // inner loop is over all other atoms in the system: @@ -197,144 +194,110 @@ namespace oopse { for(int m = -lNumber_; m <= lNumber_; m++){ sphericalHarmonic.setM(m); - QBar_lm[m] += sphericalHarmonic.getValueAt(costheta,phi); + q_lm[m] += sphericalHarmonic.getValueAt(costheta, phi); } nBonds++; } } } } - - // Normalize Qbar2 - for (int m = -lNumber_;m <= lNumber_; m++){ - QBar_lm[m] /= nBonds; - std::cout << "m = " << m << " QBLM = " << QBar_lm[m] << "\n"; - } + RealType ql = 0.0; + for(int m=-lNumber_; m<=lNumber_; m++) { + ql += norm(QBar_lm[m]); + } + ql *= 4.0*NumericConstant::PI/(RealType)(2*lNumber_+1); + collectHistogram(sqrt(ql)/(RealType)nBonds); - // Find second order invariant Q_l - - QSq_l = 0.0; - for (int m = -lNumber_; m <= lNumber_; m++){ - QSq_l += norm(QBar_lm[m]); + Nbonds += nBonds; + for (int m=-lNumber_; m<=lNumber_; m++) { + QBar_lm[m] += q_lm[m]; } - std::cout << "qsq_l = " << QSq_l << "\n"; - Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / - (2.0*(RealType)lNumber_ + 1.0)); - - // Find Third Order Invariant W_l - - W_l = 0.0; - for (int m1 = -lNumber_; m1 <= lNumber_; m1++) { - // Zero work array - for (int ii = 0; ii < mSize_; ii++){ - THRCOF[i] = 0.0; - } - // Get Wigner coefficients - m1Pass = (double)m1; - Wigner3jm(&l_, &l_, &l_, &m1Pass, &m2Min, &m2Max, THRCOF, &mSize_, &error); - for (int m_index = 1; m_index < (int)(m2Max - m2Min-1.0); m_index++) { - m2 = floor(m2Min) + m_index - 1; - m3 = -m1-m2; - W_l += THRCOF[m_index]*QBar_lm[m1]*QBar_lm[m2]*QBar_lm[m3]; - } - } - - W_l_hat = W_l / pow(QSq_l, 1.5); - - // accumulate histogram data for Q_l and W_l_hat: - - std::cout << "Ql = " << Q_l << " Wl = " << W_l_hat << "\n"; - collectHistogram(Q_l, real(W_l_hat)); - } } + + // Normalize Qbar2 + for (int m = -lNumber_;m <= lNumber_; m++){ + QBar_lm[m] /= Nbonds; + } - writeOrderParameter(); + // Find second order invariant Q_l + QSq_l = 0.0; + for (int m = -lNumber_; m <= lNumber_; m++){ + QSq_l += norm(QBar_lm[m]); + } + + std::cout << "qsl = " << QSq_l << "\n"; + Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / (RealType)(2*lNumber_ + 1)); + + // Find Third Order Invariant W_l + + W_l = 0.0; + for (int m1 = -lNumber_; m1 <= lNumber_; m1++) { + // Zero work array + for (int ii = 0; ii < mSize_; ii++){ + THRCOF[ii] = 0.0; + } + // Get Wigner coefficients + m1Pass = (double)m1; + + Wigner3jm(&l_, &l_, &l_, + &m1Pass, &m2Min, &m2Max, + THRCOF, &mSize_, &error); + + for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) { + m2 = (int)floor(m2Min) + mmm; + m3 = -m1-m2; + W_l += THRCOF[mmm] * QBar_lm[m1] * QBar_lm[m2] * QBar_lm[m3]; + } + } + + W_l_hat = W_l / pow(QSq_l, 1.5); + + writeOrderParameter(Q_l, real(W_l_hat)); } + void BondOrderParameter::collectHistogram(RealType q_l) { - void BondOrderParameter::collectHistogram(RealType Q_l, RealType W_l_hat) { - - if (Q_l >= MinQ_ && Q_l < MaxQ_) { - int qbin = (Q_l - MinQ_) / deltaQ_; + if (q_l >= MinQ_ && q_l < MaxQ_) { + int qbin = (q_l - MinQ_) / deltaQ_; Q_histogram_[qbin] += 1; Qcount_++; - sumQ_ += Q_l; - sumQ2_ += Q_l * Q_l; + sumQ_ += q_l; + sumQ2_ += q_l * q_l; } else { sprintf( painCave.errMsg, - "Q_l value outside reasonable range\n"); + "q_l value outside reasonable range\n"); painCave.severity = OOPSE_ERROR; painCave.isFatal = 1; simError(); } - if (W_l_hat >= MinW_ && W_l_hat < MaxW_) { - int wbin = (W_l_hat - MinW_) / deltaW_; - W_histogram_[wbin] += 1; - Wcount_++; - sumW_ += W_l_hat; - sumW2_ += W_l_hat*W_l_hat; - } else { - sprintf( painCave.errMsg, - "W_l_hat value outside reasonable range\n"); - painCave.severity = OOPSE_ERROR; - painCave.isFatal = 1; - simError(); - } } - void BondOrderParameter::writeOrderParameter() { - std::ofstream osq((getOutputFileName() + "q").c_str()); + void BondOrderParameter::writeOrderParameter(RealType ql, RealType Wlhat) { - if (osq.is_open()) { + std::ofstream os(getOutputFileName().c_str()); - RealType qAvg = sumQ_ / (RealType) Qcount_; - RealType qStdDev = sumQ2_ / (RealType) Qcount_ - qAvg*qAvg; + if (os.is_open()) { - osq << "# Bond Order Parameter Q_" << lNumber_ << "\n"; - osq << "# selection: (" << selectionScript_ << ")\n"; - osq << "# : " << qAvg << "\n"; - osq << "# std. dev.: " << qStdDev << "\n"; - + os << "# Bond Order Parameters\n"; + os << "# selection: (" << selectionScript_ << ")\n"; + os << "# \n"; + os << "# : " << ql << "\n"; + os << "# : " << Wlhat << "\n"; // Normalize by number of frames and write it out: for (int i = 0; i < Q_histogram_.size(); ++i) { RealType Qval = MinQ_ + (i + 0.5) * deltaQ_; - osq << Qval << "\t" << Q_histogram_[i] / frameCounter_ << "\n"; + osq << Qval << "\t" << (RealType)Q_histogram_[i] / (RealType)Qcount_ << "\n"; } - - osq.close(); - } else { - sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n", - (getOutputFileName() + "q").c_str()); - painCave.isFatal = 1; - simError(); - } - std::ofstream osw((getOutputFileName() + "w").c_str()); + os.close(); - if (osw.is_open()) { - - RealType wAvg = sumW_ / (RealType) Wcount_; - RealType wStdDev = sumW2_ / (RealType) Wcount_ - wAvg*wAvg; - - osw << "# Bond Order Parameter W_" << lNumber_ << "\n"; - osw << "# selection: (" << selectionScript_ << ")\n"; - osw << "# : " << wAvg << "\n"; - osw << "# std. dev.: " << wStdDev << "\n"; - - // Normalize by number of frames and write it out: - for (int i = 0; i < W_histogram_.size(); ++i) { - RealType Wval = MinW_ + (i + 0.5) * deltaW_; - osw << Wval << "\t" << W_histogram_[i] / frameCounter_ << "\n"; - } - - osw.close(); } else { sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n", - (getOutputFileName() + "w").c_str()); + getOutputFileName().c_str()); painCave.isFatal = 1; simError(); }