--- trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/09/26 01:30:32	1052
+++ trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/09/26 16:08:44	1054
@@ -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.18 2006-09-26 16:08:44 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"
@@ -85,13 +86,56 @@ namespace oopse {
     MinW_ = -0.25;
     MaxW_ = 0.25;
     deltaW_ = (MaxW_ - MinW_) / nbins;
-  }
 
+    // Make arrays for Wigner3jm
+    double* THRCOF = new double[2*lMax_+1];
+    // Variables for Wigner routine
+    double lPass, m1Pass, m2m, m2M;
+    int error, mSize;
+    mSize = 2*lMax_+1;
+
+    for (int l = 0; l <= lMax_; l++) {
+      lPass = (double)l;
+      for (int m1 = -l; m1 <= l; m1++) {
+        m1Pass = (double)m1;
+
+        std::pair<int,int> lm = std::make_pair(l, m1);
+        
+        // Zero work array
+        for (int ii = 0; ii < 2*l + 1; ii++){
+          THRCOF[ii] = 0.0;
+        }
+            
+        // Get Wigner coefficients
+        Wigner3jm(&lPass, &lPass, &lPass, 
+                  &m1Pass, &m2m, &m2M, 
+                  THRCOF, &mSize, &error);
+        
+        m2Min[lm] = (int)floor(m2m);
+        m2Max[lm] = (int)floor(m2M);
+        
+        for (int mmm = 0; mmm < (int)(m2M - m2m); mmm++) {
+          w3j[lm].push_back(THRCOF[mmm]);
+        }
+      }
+    }
+    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++) {
@@ -128,12 +172,6 @@ namespace oopse {
     int nBonds, Nbonds;
     SphericalHarmonic sphericalHarmonic;
     int i, j;
-    // Make arrays for Wigner3jm
-    double* THRCOF = new double[2*lMax_+1];
-    // Variables for Wigner routine
-    double lPass, m1Pass, m2Min, m2Max;
-    int error, m1, m2, m3, mSize;
-    mSize = 2*lMax_+1;
 
     DumpReader reader(info_, dumpFilename_);    
     int nFrames = reader.getNFrames();
@@ -245,27 +283,14 @@ namespace oopse {
     
         for (int l = 0; l <= lMax_; l++) {
           w[l] = 0.0;
-          lPass = (double)l;
           for (int m1 = -l; m1 <= l; m1++) {
-            // Zero work array
-            for (int ii = 0; ii < 2*l + 1; ii++){
-              THRCOF[ii] = 0.0;
+            std::pair<int,int> lm = std::make_pair(l, m1);
+            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] * 
+                q[std::make_pair(l,m2)] *  q[std::make_pair(l,m3)];
             }
-            // Get Wigner coefficients
-            m1Pass = (double)m1;
-            
-            Wigner3jm(&lPass, &lPass, &lPass, 
-                      &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] * 
-                q[std::make_pair(l,m1)] * 
-                q[std::make_pair(l,m2)] * 
-                q[std::make_pair(l,m3)];
-            }
           }
           
           w_hat[l] = w[l] / pow(q2[l], 1.5);
@@ -303,27 +328,14 @@ namespace oopse {
     
     for (int l = 0; l <= lMax_; l++) {
       W[l] = 0.0;
-      lPass = (double)l;
       for (int m1 = -l; m1 <= l; m1++) {
-        // Zero work array
-        for (int ii = 0; ii < 2*l + 1; ii++){
-          THRCOF[ii] = 0.0;
+        std::pair<int,int> lm = std::make_pair(l, m1);
+        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] * 
+            QBar[std::make_pair(l,m2)] * QBar[std::make_pair(l,m3)];
         }
-        // Get Wigner coefficients
-        m1Pass = (double)m1;
-        
-        Wigner3jm(&lPass, &lPass, &lPass, 
-                  &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[std::make_pair(l,m1)] * 
-            QBar[std::make_pair(l,m2)] * 
-            QBar[std::make_pair(l,m3)];
-        }
       }
       
       W_hat[l] = W[l] / pow(Q2[l], 1.5);
@@ -384,7 +396,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];
         }
         osq << "\n";
       }
@@ -412,7 +425,8 @@ namespace oopse {
         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];
         }
         osw << "\n";
       }