--- trunk/src/applications/staticProps/BondOrderParameter.cpp 2006/09/20 20:13:40 1041 +++ trunk/src/applications/staticProps/BondOrderParameter.cpp 2006/09/21 20:43:17 1046 @@ -51,7 +51,7 @@ #include "io/DumpReader.hpp" #include "primitives/Molecule.hpp" #include "utils/NumericConstant.hpp" -#include "math/RealSphericalHarmonic.hpp" +#include "math/SphericalHarmonic.hpp" namespace oopse { @@ -76,7 +76,7 @@ namespace oopse { // 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); @@ -104,25 +104,34 @@ namespace oopse { Molecule* mol; Atom* atom; RigidBody* rb; + int myIndex; SimInfo::MoleculeIterator mi; Molecule::RigidBodyIterator rbIter; Molecule::AtomIterator ai; StuntDouble* sd; - RealType theta; + Vector3d vec; + RealType costheta; RealType phi; RealType r; RealType dist; - std::map QBar_lm; + std::map QBar_lm; RealType QSq_l; RealType Q_l; + ComplexType W_l; + ComplexType W_l_hat; int nBonds; - RealSphericalHarmonic sphericalHarmonic; + SphericalHarmonic sphericalHarmonic; int i, j; - + // Make arrays for Wigner3jm + double* THRCOF = new double[mSize_]; + // Variables for Wigner routine + double l_ = (double)lNumber_; + double m1Pass, m2Min, m2Max; + int error, m1, m2, m3; + // Set the l for the spherical harmonic, it doesn't change sphericalHarmonic.setL(lNumber_); - DumpReader reader(info_, dumpFilename_); int nFrames = reader.getNFrames(); frameCounter_ = 0; @@ -151,6 +160,8 @@ namespace oopse { 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; @@ -166,32 +177,35 @@ namespace oopse { for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + if (atom->getGlobalIndex() != myIndex) { - Vector3d vec = sd->getPos() - atom->getPos(); - currentSnapshot_->wrapVector(vec); - - // Calculate "bonds" and build Q_lm(r) where - // Q_lm = Y_lm(theta(r),phi(r)) - // The spherical harmonics are wrt any arbitrary coordinate - // system, we choose standard spherical coordinates - - r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2)); - - // Check to see if neighbor is in bond cutoff - - if (r < rCut_) { - theta = atan2(vec.y(), vec.x()); - phi = acos(vec.z()/r); - for(int m = -lNumber_; m <= lNumber_; m++){ - sphericalHarmonic.setM(m); - QBar_lm[m] += sphericalHarmonic.getValueAt(theta,phi); - } - nBonds++; - } + vec = sd->getPos() - atom->getPos(); + currentSnapshot_->wrapVector(vec); + + // Calculate "bonds" and build Q_lm(r) where + // Q_lm = Y_lm(theta(r),phi(r)) + // The spherical harmonics are wrt any arbitrary coordinate + // system, we choose standard spherical coordinates + + r = vec.length(); + + // Check to see if neighbor is in bond cutoff + + if (r < rCut_) { + costheta = vec.z() / r; + phi = atan2(vec.y(), vec.x()); + + for(int m = -lNumber_; m <= lNumber_; m++){ + sphericalHarmonic.setM(m); + QBar_lm[m] += sphericalHarmonic.getValueAt(costheta,phi); + } + nBonds++; + } + } } } - // Normalize Qbar + // Normalize Qbar2 for (int m = -lNumber_;m <= lNumber_; m++){ QBar_lm[m] /= nBonds; } @@ -200,42 +214,40 @@ namespace oopse { QSq_l = 0.0; for (int m = -lNumber_; m <= lNumber_; m++){ - QSq_l += pow(QBar_lm[m], 2); + QSq_l += norm(QBar_lm[m]); } - Q_l = sqrt(QSq_l*(4.0 * NumericConstant::PI / (2.0*(RealType)lNumber_ + 1))); - - // Find Third Order Invariant W_l - // Make arrays for Wigner3jm - double* THRCOF = new double[mSize_]; - // Variables for Wigner routine - double l_ = (double)lNumber_; - double m1Pass, m2Min, m2Max; - int error, m1, m2, m3; + Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / + (2.0*(RealType)lNumber_ + 1.0)); - RealType W_l; - RealType W_l_hat; + // 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; + THRCOF[ii] = 0.0; } // Get Wigner coefficients m1Pass = (double)m1; - Wigner3jm(&l_, &l_, &l_, &m1Pass, &m2Min, &m2Max, THRCOF, &mSize_, &error); - for (int m_index = 1; i < (int)(m2Max - m2Min-1.0); m_index++) { - m2 = floor(m2Min) + m_index - 1; + + 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[m_index]*QBar_lm[m1+lNumber_]*QBar_lm[m2+lNumber_]*QBar_lm[m3+lNumber_]; + W_l += THRCOF[mmm] * 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: - collectHistogram(Q_l, W_l_hat); + std::cout << "Ql = " << Q_l << " Wl = " << W_l_hat << "\n"; + collectHistogram(Q_l, real(W_l_hat)); } } @@ -293,7 +305,7 @@ namespace oopse { // 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(); @@ -319,7 +331,7 @@ namespace oopse { // 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 << Wval << "\t" << (RealType)W_histogram_[i] / (RealType)Wcount_ << "\n"; } osw.close();