--- trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/09/26 16:08:44	1054
+++ trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/12/20 21:03:11	1100
@@ -43,7 +43,7 @@
  *
  *  Created by J. Daniel Gezelter on 09/26/06.
  *  @author  J. Daniel Gezelter
- *  @version $Id: BondOrderParameter.cpp,v 1.18 2006-09-26 16:08:44 gezelter Exp $
+ *  @version $Id: BondOrderParameter.cpp,v 1.22 2006-12-20 21:03:11 gezelter Exp $
  *
  */
  
@@ -83,8 +83,8 @@ namespace oopse {
     // W_6 for icosahedral clusters is 11 / sqrt(4199) = 0.169754, so we'll
     // use values for MinW_ and MaxW_ that are slightly larger than this:
 
-    MinW_ = -0.25;
-    MaxW_ = 0.25;
+    MinW_ = -1.1;
+    MaxW_ = 1.1;
     deltaW_ = (MaxW_ - MinW_) / nbins;
 
     // Make arrays for Wigner3jm
@@ -105,7 +105,7 @@ namespace oopse {
         for (int ii = 0; ii < 2*l + 1; ii++){
           THRCOF[ii] = 0.0;
         }
-            
+
         // Get Wigner coefficients
         Wigner3jm(&lPass, &lPass, &lPass, 
                   &m1Pass, &m2m, &m2M, 
@@ -114,7 +114,7 @@ namespace oopse {
         m2Min[lm] = (int)floor(m2m);
         m2Max[lm] = (int)floor(m2M);
         
-        for (int mmm = 0; mmm < (int)(m2M - m2m); mmm++) {
+        for (int mmm = 0; mmm <= (int)(m2M - m2m); mmm++) {
           w3j[lm].push_back(THRCOF[mmm]);
         }
       }
@@ -186,6 +186,7 @@ namespace oopse {
     Q.resize(lMax_+1);
     W.resize(lMax_+1);
     W_hat.resize(lMax_+1);
+    Nbonds = 0;
 
     for (int istep = 0; istep < nFrames; istep += step_) {
       reader.readFrame(istep);
@@ -230,7 +231,9 @@ namespace oopse {
             if (atom->getGlobalIndex() != myIndex) {
 
               vec = sd->getPos() - atom->getPos();       
-              currentSnapshot_->wrapVector(vec);
+
+              if (usePeriodicBoundaryConditions_) 
+                currentSnapshot_->wrapVector(vec);
               
               // Calculate "bonds" and build Q_lm(r) where 
               //      Q_lm = Y_lm(theta(r),phi(r))                
@@ -258,34 +261,23 @@ namespace oopse {
           }
         }
         
-        
-        for (int l = 0; l <= lMax_; l++) {         
-          q_l[l] = 0.0;
-          for(int m = -l; m <= l; m++) { 
-            q_l[l] += norm(q[std::make_pair(l,m)]);
-          }     
-          q_l[l] *= 4.0*NumericConstant::PI/(RealType)(2*l + 1);
-          q_l[l] = sqrt(q_l[l])/(RealType)nBonds;
-        }
-
-        // Find second order invariant Q_l
         
         for (int l = 0; l <= lMax_; l++) {
           q2[l] = 0.0;
           for (int m = -l; m <= l; m++){
+            q[std::make_pair(l,m)] /= (RealType)nBonds;            
             q2[l] += norm(q[std::make_pair(l,m)]);
           }
-          q_l[l] = sqrt(q2[l] * 4.0 * NumericConstant::PI / 
-                        (RealType)(2*l + 1))/(RealType)nBonds;
+          q_l[l] = sqrt(q2[l] * 4.0 * NumericConstant::PI / (RealType)(2*l + 1));
         }
-
+        
         // Find Third Order Invariant W_l
     
         for (int l = 0; l <= lMax_; l++) {
           w[l] = 0.0;
           for (int m1 = -l; m1 <= l; m1++) {
             std::pair<int,int> lm = std::make_pair(l, m1);
-            for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) {
+            for (int mmm = 0; mmm <= (m2Max[lm] - m2Min[lm]); mmm++) {
               int m2 = m2Min[lm] + mmm;
               int m3 = -m1-m2;
               w[l] += w3j[lm][mmm] * q[lm] * 
@@ -301,7 +293,7 @@ namespace oopse {
         Nbonds += nBonds;
         for (int l = 0; l <= lMax_;  l++) {
           for (int m = -l; m <= l; m++) {
-            QBar[std::make_pair(l,m)] += q[std::make_pair(l,m)];
+            QBar[std::make_pair(l,m)] += (RealType)nBonds*q[std::make_pair(l,m)];
           }
         }
       }
@@ -330,7 +322,7 @@ namespace oopse {
       W[l] = 0.0;
       for (int m1 = -l; m1 <= l; m1++) {
         std::pair<int,int> lm = std::make_pair(l, m1);
-        for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) {
+        for (int mmm = 0; mmm <= (m2Max[lm] - m2Min[lm]); mmm++) {
           int m2 = m2Min[lm] + mmm;
           int m3 = -m1-m2;
           W[l] += w3j[lm][mmm] * QBar[lm] * 
@@ -368,7 +360,7 @@ namespace oopse {
         Wcount_[l]++;      
       } else {
         sprintf( painCave.errMsg,
-                 "Re[w_hat] value outside reasonable range\n");
+                 "Re[w_hat] value (%lf) outside reasonable range\n", real(what[l]));
         painCave.severity = OOPSE_ERROR;
         painCave.isFatal = 1;
         simError();  
@@ -396,8 +388,8 @@ namespace oopse {
         RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;               
         osq << Qval;
         for (int l = 0; l <= lMax_; l++) {
-          osq << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)] / 
-            (RealType)Qcount_[l];
+
+          osq << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)]/(RealType)Qcount_[l]/deltaQ_;
         }
         osq << "\n";
       }
@@ -418,15 +410,15 @@ namespace oopse {
       osw << "# selection: (" << selectionScript_ << ")\n";
       osw << "# \n";
       for (int l = 0; l <= lMax_; l++) {
-        osw << "# <W_" << l << ">: " << real(What[l]) << "\n";
+        osw << "# <W_" << l << ">: " << real(What[l]) << "\t" << imag(What[l]) << "\n";
       }
       // Normalize by number of frames and write it out:
       for (int i = 0; i < nBins_; ++i) {
         RealType Wval = MinW_ + (i + 0.5) * deltaW_;               
         osw << Wval;
         for (int l = 0; l <= lMax_; l++) {
-          osw << "\t" << (RealType)W_histogram_[std::make_pair(i,l)] / 
-            (RealType)Wcount_[l];
+
+          osw << "\t" << (RealType)W_histogram_[std::make_pair(i,l)]/(RealType)Wcount_[l]/deltaW_;
         }
         osw << "\n";
       }