--- trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/06/27 16:19:28	994
+++ trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/07/07 14:18:36	1001
@@ -39,20 +39,29 @@
  * such damages.
  */
 
-#include "applications/staticProps/P2OrderParameter.hpp"
+
+/* 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"
 #include "io/DumpReader.hpp"
 #include "primitives/Molecule.hpp"
 #include "utils/NumericConstant.hpp"
+#include "math/RealSphericalHarmonic.hpp"
 namespace oopse {
 
 
-P2OrderParameter::P2OrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
+BondOrderParameter::BondOrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1,
+ const std::string& sele2, double rCut, int lNumber)
   : StaticAnalyser(info, filename),
-    selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info), 
-    seleMan1_(info), seleMan2_(info){
+    selectionScript1_(sele1), evaluator1_(info), 
+    seleMan1_(info){
 
-    setOutputName(getPrefix(filename) + ".p2");
+    setOutputName(getPrefix(filename) + ".obo");
         
     evaluator1_.loadScriptString(sele1);
     evaluator2_.loadScriptString(sele2);
@@ -67,30 +76,16 @@ P2OrderParameter::P2OrderParameter(SimInfo* info, cons
         simError();  
     }
 
-    if (!evaluator2_.isDynamic()) {
-      seleMan2_.setSelectionSet(evaluator2_.evaluate());
-    }else {
-        sprintf( painCave.errMsg,
-                 "--sele2 must be static selection\n");
-        painCave.severity = OOPSE_ERROR;
-        painCave.isFatal = 1;
-        simError();  
-    }
+/* Set up cutoff radius and type of order parameter we are calcuating*/
+	lNumber_ = lNumber;
+	rCut_ = rCut;
+	mSize_ = 2*lNumber_+1;
 
-    if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
-        sprintf( painCave.errMsg,
-                 "The number of selected Stuntdoubles are not the same in --sele1 and sele2\n");
-        painCave.severity = OOPSE_ERROR;
-        painCave.isFatal = 1;
-        simError();  
-
-    }
-
   int i;
   int j;
   StuntDouble* sd1;
   StuntDouble* sd2;
-  for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
+  for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan1_.beginSelected(j);
      sd1 != NULL && sd2 != NULL;
      sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
 
@@ -100,19 +95,30 @@ P2OrderParameter::P2OrderParameter(SimInfo* info, cons
     
   }
 
-void P2OrderParameter::process() {
+void BondOrderParameter::process() {
   Molecule* mol;
   RigidBody* rb;
   SimInfo::MoleculeIterator mi;
   Molecule::RigidBodyIterator rbIter;
+  RealType theta;
+  RealType phi;
+  RealType r;
+  RealType dist;
+  RealType* QBar_lm;
+  int nBonds;
+  RealSphericalHarmonic sphericalHarmonic;
   
+  
   DumpReader reader(info_, dumpFilename_);    
   int nFrames = reader.getNFrames();
 
+   /*Set the l for the spherical harmonic, it doesn't change*/
+	sphericalHarmonic.setL(lNumber_);
+
   for (int i = 0; i < nFrames; i += step_) {
     reader.readFrame(i);
     currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
-
+	nBonds = 0;
     
     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
         //change the positions of atoms which belong to the rigidbodies
@@ -122,62 +128,108 @@ void P2OrderParameter::process() {
         
     }      
 
-      Mat3x3d orderTensor(0.0);
+      /* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
       for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
           Vector3d vec = j->first->getPos() - j->second->getPos();
           currentSnapshot_->wrapVector(vec);
-          vec.normalize();
-          orderTensor +=outProduct(vec, vec);
+ /* The spherical harmonics are wrt any arbitray coordiate sysetm, 
+  * 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 cuttoff*/
+  		if (r<rCut_){  		
+  	    theta = atan(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++;
+  		}  
       }
       
-      orderTensor /= sdPairs_.size();
-      orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
+     /*Normalize Qbar*/
+     for (int m = -lNumber_;m <= lNumber_; m++){
+      QBar_lm(m) = QBar_lm(m)/nBonds;	
+     }
       
-      Vector3d eigenvalues;
-      Mat3x3d eigenvectors;    
-      Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors);
-      
-      int which;
-      RealType maxEval = 0.0;
-      for(int k = 0; k< 3; k++){
-        if(fabs(eigenvalues[k]) > maxEval){
-          which = k;
-          maxEval = fabs(eigenvalues[k]);
-        }
-      }
-      RealType p2 = 1.5 * maxEval;
-      
-      //the eigen vector is already normalized in SquareMatrix3::diagonalize
-      Vector3d director = eigenvectors.getColumn(which);
-      if (director[0] < 0) {
-          director.negate();
-      }   
+    
+  }
+  
+  /*Normalize by number of frames*/
+    for (int m = -lNumber_;m <= lNumber_; m++){
+      QBar_lm(m) = QBar_lm(m)/nFrames;	
+     }
+     
+     
+     
+     /* Find second order invariant Q_l*/
+     
+       for (int m = -lNumber_;m <= lNumber_; m++){
+     	 QSq_l += pow(QBar_lm(m),2);	
+     	}
+     Q_l = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);
+     
+     /* Find Third Order Invariant W_l*/
+     for (int m = -lNumber_;m<= lNumber_;m++){
+     	
+     	
+     }
+     
+     
+  writeOrderParameter();
+  
+}
 
-      RealType angle = 0.0;
-      for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
-          Vector3d vec = j->first->getPos() - j->second->getPos();
-          currentSnapshot_->wrapVector(vec);
-          vec.normalize();
+void BondOrderParameter::initalizeHistogram() {
+    std::fill(histogram_.begin(), histogram_.end(), 0);
+  }
 
-          angle += acos(dot(vec, director)) ;
-      }
-      angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
+ void BondOrderParameter::processHistogram() {
 
-       OrderParam param;
-       param.p2 = p2;
-       param.director = director;
-       param.angle = angle;
+    int nPairs = getNPairs();
+    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
+    RealType pairDensity = nPairs /volume * 2.0;
+    RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
 
-        orderParams_.push_back(param);       
+    for(int i = 0 ; i < histogram_.size(); ++i){
+
+      RealType rLower = i * deltaR_;
+      RealType rUpper = rLower + deltaR_;
+      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
+      RealType nIdeal = volSlice * pairConstant;
+
+      avgGofr_[i] += histogram_[i] / nIdeal;    
+    }
+
+  }
+
+  void BondOrderParameter::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
+
+    if (sd1 == sd2) {
+      return;
+    }
     
+    Vector3d pos1 = sd1->getPos();
+    Vector3d pos2 = sd2->getPos();
+    Vector3d r12 = pos2 - pos1;
+    currentSnapshot_->wrapVector(r12);
+
+    RealType distance = r12.length();
+
+    if (distance < len_) {
+      int whichBin = distance / deltaR_;
+      histogram_[whichBin] += 2;
+    }
   }
 
-  writeP2();
-  
-}
 
-void P2OrderParameter::writeP2() {
 
+
+
+
+void BondOrderParameter::writeOrderParameter() {
+
     std::ofstream os(getOutputFileName().c_str());
     os << "#radial distribution function\n";
     os<< "#selection1: (" << selectionScript1_ << ")\t";