--- trunk/src/applications/staticProps/P2OrderParameter.cpp	2011/03/03 20:32:49	1542
+++ trunk/src/applications/staticProps/P2OrderParameter.cpp	2011/10/18 13:44:44	1657
@@ -61,7 +61,6 @@ namespace OpenMD {
     if (!evaluator1_.isDynamic()) {
       seleMan1_.setSelectionSet(evaluator1_.evaluate());
     }
-    
   }
 
   P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename, 
@@ -122,8 +121,7 @@ namespace OpenMD {
     SimInfo::MoleculeIterator mi;
     Molecule::RigidBodyIterator rbIter;
     StuntDouble* sd;
-    int i;
-
+    int i, ii;
   
     DumpReader reader(info_, dumpFilename_);    
     int nFrames = reader.getNFrames();
@@ -131,7 +129,7 @@ namespace OpenMD {
     for (int i = 0; i < nFrames; i += step_) {
       reader.readFrame(i);
       currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
-
+      
       for (mol = info_->beginMolecule(mi); mol != NULL; 
            mol = info_->nextMolecule(mi)) {
         //change the positions of atoms which belong to the rigidbodies
@@ -148,8 +146,8 @@ namespace OpenMD {
         if  (evaluator1_.isDynamic()) 
           seleMan1_.setSelectionSet(evaluator1_.evaluate());
         
-        for (sd = seleMan1_.beginSelected(i); sd != NULL; 
-             sd = seleMan1_.nextSelected(i)) {
+        for (sd = seleMan1_.beginSelected(ii); sd != NULL; 
+             sd = seleMan1_.nextSelected(ii)) {
           if (sd->isDirectional()) {
             Vector3d vec = sd->getA().getColumn(2);
             vec.normalize();
@@ -173,10 +171,12 @@ namespace OpenMD {
         orderTensor /= sdPairs_.size();
       }
       
+
       orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
       
       Vector3d eigenvalues;
       Mat3x3d eigenvectors;    
+
       Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors);
       
       int which;
@@ -196,35 +196,40 @@ namespace OpenMD {
       }   
 
       RealType angle = 0.0;
-
+      RealType p4 = 0.0;
       
       if (doVect_) {
-        for (sd = seleMan1_.beginSelected(i); sd != NULL; 
-             sd = seleMan1_.nextSelected(i)) {
+        for (sd = seleMan1_.beginSelected(ii); sd != NULL; 
+             sd = seleMan1_.nextSelected(ii)) {
           if (sd->isDirectional()) {
             Vector3d vec = sd->getA().getColumn(2);
             vec.normalize();
-            angle += acos(dot(vec, director));
+            RealType myDot = dot(vec, director);
+            angle += acos(myDot);
+            p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0;
           }
         }
         angle = angle/(seleMan1_.getSelectionCount()*NumericConstant::PI)*180.0;
-        
+        p4 /= (seleMan1_.getSelectionCount());
       } else {
         for (vector<pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
           Vector3d vec = j->first->getPos() - j->second->getPos();
           if (usePeriodicBoundaryConditions_)
             currentSnapshot_->wrapVector(vec);
-          vec.normalize();
-          
-          angle += acos(dot(vec, director)) ;
+          vec.normalize();          
+          RealType myDot = dot(vec, director);
+          angle += acos(myDot);
+          p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0;
         }
         angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
+        p4 /= (seleMan1_.getSelectionCount());
       }
 
       OrderParam param;
       param.p2 = p2;
       param.director = director;
       param.angle = angle;
+      param.p4 = p4;
 
       orderParams_.push_back(param);       
     
@@ -242,14 +247,15 @@ namespace OpenMD {
     if (!doVect_) {
       os << "selection2: (" << selectionScript2_ << ")\n";
     }
-    os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";    
+    os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\t<p4>\n";    
 
     for (size_t i = 0; i < orderParams_.size(); ++i) {
       os <<  orderParams_[i].p2 << "\t"
          <<  orderParams_[i].director[0] << "\t"
          <<  orderParams_[i].director[1] << "\t"
          <<  orderParams_[i].director[2] << "\t"
-         <<  orderParams_[i].angle << "\n";
+         <<  orderParams_[i].angle << "\t"
+         <<  orderParams_[i].p4 << "\n";
 
     }