--- trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/09/26 14:33:56	1053
+++ trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/11/21 20:44:54	1094
@@ -37,14 +37,15 @@
  * arising out of the use of or inability to use software, even if the
  * University of Notre Dame has been advised of the possibility of
  * such damages.
+ *
+ *  BondOrderParameter.cpp
+ *  OOPSE-4
+ *
+ *  Created by J. Daniel Gezelter on 09/26/06.
+ *  @author  J. Daniel Gezelter
+ *  @version $Id: BondOrderParameter.cpp,v 1.21 2006-11-21 20:44:54 gezelter Exp $
+ *
  */
-
-
-/* Creates orientational bond order parameters as outlined by
- *     Bond-orientaional order in liquids and glasses, Steinhart,Nelson,Ronchetti
- *     Phys Rev B, 28,784,1983
- * 
- */
  
 #include "applications/staticProps/BondOrderParameter.hpp"
 #include "utils/simError.h"
@@ -118,13 +119,23 @@ namespace oopse {
         }
       }
     }
+    delete [] THRCOF;
+    THRCOF = NULL;
   }
   
   BondOrderParameter::~BondOrderParameter() {
     Q_histogram_.clear();
     W_histogram_.clear();
+    for (int l = 0; l <= lMax_; l++) {
+      for (int m = -l; m <= l; m++) {
+        w3j[std::make_pair(l,m)].clear();
+      }
+    }
+    w3j.clear();
+    m2Min.clear();
+    m2Max.clear();
   }
-
+  
   void BondOrderParameter::initalizeHistogram() {
     for (int bin = 0; bin < nBins_; bin++) {
       for (int l = 0; l <= lMax_; l++) {
@@ -175,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);
@@ -219,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))                
@@ -248,26 +262,15 @@ 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++) {
@@ -290,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)];
           }
         }
       }
@@ -357,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();  
@@ -385,7 +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";
       }
@@ -406,14 +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";
       }