--- trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/09/18 21:31:23	1038
+++ trunk/src/applications/staticProps/BondOrderParameter.cpp	2006/09/21 18:04:52	1044
@@ -45,189 +45,298 @@
  *     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"
+#include "math/SphericalHarmonic.hpp"
+
 namespace oopse {
 
+  BondOrderParameter::BondOrderParameter(SimInfo* info, 
+                                         const std::string& filename, 
+                                         const std::string& sele,
+                                         double rCut, int lNumber, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
+    
+    setOutputName(getPrefix(filename) + ".bo");
 
-  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), evaluator1_(info),
-      seleMan1_(info){
-
-    setOutputName(getPrefix(filename) + ".obo");
-
-    evaluator1_.loadScriptString(sele1);
-    evaluator2_.loadScriptString(sele2);
-
-    if (!evaluator1_.isDynamic()) {
-      seleMan1_.setSelectionSet(evaluator1_.evaluate());
-    }else {
-      sprintf( painCave.errMsg,
-	       "--sele1 must be static selection\n");
-      painCave.severity = OOPSE_ERROR;
-      painCave.isFatal = 1;
-      simError();
+    evaluator_.loadScriptString(sele);
+    if (!evaluator_.isDynamic()) {
+      seleMan_.setSelectionSet(evaluator_.evaluate());
     }
 
-    /* Set up cutoff radius and type of order parameter we are calcuating*/
+    // Set up cutoff radius and order of the Legendre Polynomial:
+
     lNumber_ = lNumber;
     rCut_ = rCut;
-    mSize_ = 2*lNumber_+1;
+    mSize_ = 2*lNumber_+1;    
 
-    int i;
-    int j;
-    StuntDouble* sd1;
-    StuntDouble* sd2;
-    for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan1_.beginSelected(j);
-	 sd1 != NULL && sd2 != NULL;
-	 sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
-      for (sd2 = seleMan1_.beginSelected(j),sd2
-	     sdPairs_.push_back(std::make_pair(sd1, sd2));
-	   }
+    // Q can take values from 0 to 1
 
+    MinQ_ = 0.0;
+    MaxQ_ = 1.0;
+    deltaQ_ = (MaxQ_ - MinQ_) / nbins;
+    Q_histogram_.resize(nbins);
 
-    }
+    // 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:
 
-    void BondOrderParameter::process
-      () {
-      Molecule* mol;
-      RigidBody* rb;
-      SimInfo::MoleculeIterator mi;
-      Molecule::RigidBodyIterator rbIter;
-      RealType theta;
-      RealType phi;
-      RealType r;
-      RealType dist;
-      RealType* QBar_lm;
-      RealType QSq_l;
-      int nBonds;
-      int m, m_index;
-      RealSphericalHarmonic sphericalHarmonic;
+    MinW_ = -0.18;
+    MaxW_ = 0.18;
+    deltaW_ = (MaxW_ - MinW_) / nbins;
+    W_histogram_.resize(nbins);
 
+  }
 
-      DumpReader reader(info_, dumpFilename_);
-      int nFrames = reader.getNFrames();
+  BondOrderParameter::~BondOrderParameter() {
+    Q_histogram_.clear();
+    W_histogram_.clear();
+  }
 
-      /*Set the l for the spherical harmonic, it doesn't change*/
-      sphericalHarmonic.setL(lNumber_);
+  void BondOrderParameter::initalizeHistogram() {
+    std::fill(Q_histogram_.begin(), Q_histogram_.end(), 0);
+    std::fill(W_histogram_.begin(), W_histogram_.end(), 0);
+  }
 
-      for (int i = 0; i < nFrames; i += step_) {
-	reader.readFrame(i);
-	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
-	nBonds = 0;
+  void BondOrderParameter::process() {
+    Molecule* mol;
+    Atom* atom;
+    RigidBody* rb;
+    int myIndex;
+    SimInfo::MoleculeIterator mi;
+    Molecule::RigidBodyIterator rbIter;
+    Molecule::AtomIterator ai;
+    StuntDouble* sd;
+    Vector3d vec;
+    RealType costheta;
+    RealType phi;
+    RealType r;
+    RealType dist;
+    std::map<int, ComplexType> QBar_lm;
+    RealType QSq_l;
+    RealType Q_l;
+    ComplexType W_l;
+    ComplexType W_l_hat;
+    int nBonds;
+    SphericalHarmonic sphericalHarmonic;
+    int i, j;
+    // Make arrays for Wigner3jm
+    double* THRCOF = new double[mSize_];
+    // Variables for Wigner routine
+    double l_ = (double)lNumber_;
+    double m1Pass, m2Min, m2Max;
+    int error, m1, m2, m3;
 
-	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
-	  //change the positions of atoms which belong to the rigidbodies
-	  for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
-	    rb->updateAtoms();
-	  }
+    // Set the l for the spherical harmonic, it doesn't change
+    sphericalHarmonic.setL(lNumber_);
 
-	}
+    DumpReader reader(info_, dumpFilename_);    
+    int nFrames = reader.getNFrames();
+    frameCounter_ = 0;
 
+    for (int istep = 0; istep < nFrames; istep += step_) {
+      reader.readFrame(istep);
+      frameCounter_++;
+      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
+      
+      if (evaluator_.isDynamic()) {
+        seleMan_.setSelectionSet(evaluator_.evaluate());
+      }
 
-	/* Setup QBar */
-	QBar_lm = new double[mSize_];
+      // update the positions of atoms which belong to the rigidbodies
 
-	/* 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);
-	  /* 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));
+      for (mol = info_->beginMolecule(mi); mol != NULL; 
+           mol = info_->nextMolecule(mi)) {
+        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
+             rb = mol->nextRigidBody(rbIter)) {
+          rb->updateAtoms();
+        }        
+      }      
+      
+      // outer loop is over the selected StuntDoubles:
 
-	  /* 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_index = 0; m_index < mSize_; m_index++){
-	      sphericalHarmonic.setM(m_index-lNumber_);
-	      QBar_lm(m_index) += sphericalHarmonic.getValueAt(theta,phi);
-	    }
-	    nBonds++;
-	  }
-	}
+      for (sd = seleMan_.beginSelected(i); sd != NULL; 
+           sd = seleMan_.nextSelected(i)) {
 
-	/*Normalize Qbar by number of Bonds*/
-	for ( int m_index = 0;m_index < mSize_; m_index++){
-	  QBar_lm(m_index) = QBar_lm(m_index)/nBonds;
-	}
+        myIndex = sd->getGlobalIndex();
 
+        // For this central atom, zero out nBonds and QBar_lm
 
-      }
+        nBonds = 0;
+       
+        for (int m = -lNumber_; m <= lNumber_; m++) {
+          QBar_lm[m] = 0.0;
+        }
+        
+        // inner loop is over all other atoms in the system:
+        
+        for (mol = info_->beginMolecule(mi); mol != NULL; 
+             mol = info_->nextMolecule(mi)) {
+          for (atom = mol->beginAtom(ai); atom != NULL; 
+               atom = mol->nextAtom(ai)) {
 
-      /*Normalize by number of frames*/
-      for ( int m_index = 0;m_index < mSize_; m_index++){
-	QBar_lm(m_index) = QBar_lm(m_index)/nFrames;
-      }
+            if (atom->getGlobalIndex() != myIndex) {
 
+              vec = sd->getPos() - atom->getPos();       
+              currentSnapshot_->wrapVector(vec);
+              
+              // Calculate "bonds" and build Q_lm(r) where 
+              //      Q_lm = Y_lm(theta(r),phi(r))                
+              // The spherical harmonics are wrt any arbitrary coordinate
+              // system, we choose standard spherical coordinates 
+              
+              r = vec.length();
+              
+              // Check to see if neighbor is in bond cutoff 
+              
+              if (r < rCut_) { 
+                costheta = vec.z() / r; 
+                phi = atan2(vec.y(), vec.x());
+                
+                for(int m = -lNumber_; m <= lNumber_; m++){
+                  sphericalHarmonic.setM(m);
+                  QBar_lm[m] += sphericalHarmonic.getValueAt(costheta,phi);
+                }
+                nBonds++;
+              }  
+            }
+          }
+        }
+        
+        // Normalize Qbar2
+        for (int m = -lNumber_;m <= lNumber_; m++){
+          QBar_lm[m] /= nBonds;	
+          std::cout << "m = " << m << " QBLM = " << QBar_lm[m] << "\n";
+        }
 
+        // Find second order invariant Q_l
 
-      /* Find second order invariant Q_l*/
+        QSq_l = 0.0;
+        for (int m = -lNumber_; m <= lNumber_; m++){
+          QSq_l += norm(QBar_lm[m]);
+        }
+        std::cout << "qsq_l = " << QSq_l << "\n";
+        Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / 
+                   (2.0*(RealType)lNumber_ + 1.0));
 
-      for (int m_index = 0 ;m_index <= sizeM_; m++){
-	QSq_l += pow(QBar_lm(m),2);
-      }
-      Q_l_ = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);
+        // Find Third Order Invariant W_l
+        
+        W_l = 0.0;
+        for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
+          // Zero work array
+          for (int ii = 0; ii < mSize_; ii++){
+            THRCOF[i] = 0.0;
+          }
+          // Get Wigner coefficients
+          m1Pass = (double)m1;
+          Wigner3jm(&l_, &l_, &l_, &m1Pass, &m2Min, &m2Max, THRCOF, &mSize_, &error);
+          for (int m_index = 1; m_index < (int)(m2Max - m2Min-1.0); m_index++) {
+            m2 = floor(m2Min) + m_index - 1;
+            m3 = -m1-m2;
+            W_l += THRCOF[m_index]*QBar_lm[m1]*QBar_lm[m2]*QBar_lm[m3];
+          }
+        }
+        
+        W_l_hat = W_l / pow(QSq_l, 1.5);
+        
+        // accumulate histogram data for Q_l and W_l_hat:
 
-      /* Find Third Order Invariant W_l*/
-
-      /* Make arrays for Wigner3jm */
-      double* THRCOF = new double[mSize_];
-      /* Variables for Wigner routine */
-      double l_ = (double)lNumber_;
-      double m2Min;
-      double m2Max;
-      int error;
-      int m1;
-      int m2;
-      int m3;
-
-      for (int m1 = -lNumber_;m <= lNumber_;m1++){
-	/* Zero work array */
-	for (i=0; i<mSize_;i++){
-	  THRCOF[i] = 0.0;      
-	}
-	/* Get wigner coefficients */
-	Wigner3jm(&l_,&l_,&l_,&(double)m1,&m2Min,&m2Max&,THRCOF,&mSize_,&error);
-	for (m_index=1; i<(m2Max-M2Min-1.0);m_index++){
-	  m2 = floor(m2Min) + m_index - 1;
-	  m3 = -m1-m2;
-	  W_l_ += THRCOF(m_index)*QBar_lm(m1+lNumber_)*QBar_lm(m2+lNumber_)*QBar_lm(m3+lNumber_);
-	}
+        std::cout << "Ql = " << Q_l << " Wl = " << W_l_hat << "\n";
+        collectHistogram(Q_l, real(W_l_hat));
+                
       }
+    }
+    
+    writeOrderParameter();
+    
+  }
 
 
-      writeOrderParameter();
+  void BondOrderParameter::collectHistogram(RealType Q_l, RealType W_l_hat) {
 
+    if (Q_l >= MinQ_ && Q_l < MaxQ_) {
+      int qbin = (Q_l - MinQ_) / deltaQ_;
+      Q_histogram_[qbin] += 1;
+      Qcount_++;
+      sumQ_ += Q_l;
+      sumQ2_ += Q_l * Q_l;
+    } else {
+      sprintf( painCave.errMsg,
+               "Q_l value outside reasonable range\n");
+      painCave.severity = OOPSE_ERROR;
+      painCave.isFatal = 1;
+      simError();  
     }
 
+    if (W_l_hat >= MinW_ && W_l_hat < MaxW_) {
+      int wbin = (W_l_hat - MinW_) / deltaW_;
+      W_histogram_[wbin] += 1;
+      Wcount_++;
+      sumW_  += W_l_hat;
+      sumW2_ += W_l_hat*W_l_hat;
+    } else {
+      sprintf( painCave.errMsg,
+               "W_l_hat value outside reasonable range\n");
+      painCave.severity = OOPSE_ERROR;
+      painCave.isFatal = 1;
+      simError();  
+    }
+  }  
 
-    void BondOrderParameter::writeOrderParameter() {
+  void BondOrderParameter::writeOrderParameter() {
 
-      std::ofstream os(getOutputFileName().c_str());
-      os << "#radial distribution function\n";
-      os<< "#selection1: (" << selectionScript1_ << ")\t";
-      os << "selection2: (" << selectionScript2_ << ")\n";
-      os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";
+    std::ofstream osq((getOutputFileName() + "q").c_str());
 
-      for (std::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";
+    if (osq.is_open()) {
 
+      RealType qAvg = sumQ_ / (RealType) Qcount_;
+      RealType qStdDev = sumQ2_ / (RealType) Qcount_ - qAvg*qAvg;
+      
+      osq << "# Bond Order Parameter Q_" << lNumber_ << "\n";
+      osq << "# selection: (" << selectionScript_ << ")\n";
+      osq << "# <Q_" << lNumber_ << ">: " << qAvg << "\n";
+      osq << "# std. dev.: " << qStdDev << "\n";
+      
+      // Normalize by number of frames and write it out:
+      for (int i = 0; i < Q_histogram_.size(); ++i) {
+        RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;
+        osq << Qval << "\t" << Q_histogram_[i] / frameCounter_ << "\n";
       }
-
+      
+      osq.close();
+    } else {
+      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n", 
+              (getOutputFileName() + "q").c_str());
+      painCave.isFatal = 1;
+      simError();  
     }
 
-  }
+    std::ofstream osw((getOutputFileName() + "w").c_str());
 
+    if (osw.is_open()) {
+
+      RealType wAvg = sumW_ / (RealType) Wcount_;
+      RealType wStdDev = sumW2_ / (RealType) Wcount_ - wAvg*wAvg;
+      
+      osw << "# Bond Order Parameter W_" << lNumber_ << "\n";
+      osw << "# selection: (" << selectionScript_ << ")\n";
+      osw << "# <W_" << lNumber_ << ">: " << wAvg << "\n";
+      osw << "# std. dev.: " << wStdDev << "\n";
+      
+      // Normalize by number of frames and write it out:
+      for (int i = 0; i < W_histogram_.size(); ++i) {
+        RealType Wval = MinW_ + (i + 0.5) * deltaW_;
+        osw << Wval << "\t" << W_histogram_[i] / frameCounter_ << "\n";
+      }
+      
+      osw.close();
+    } else {
+      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n", 
+              (getOutputFileName() + "w").c_str());
+      painCave.isFatal = 1;
+      simError();  
+    }
+  }
+}