--- trunk/src/applications/staticProps/GofAngle2.cpp	2014/09/26 22:22:28	2022
+++ trunk/src/applications/staticProps/GofAngle2.cpp	2014/10/02 14:35:14	2023
@@ -49,24 +49,118 @@ namespace OpenMD {
 
 namespace OpenMD {
 
-  GofAngle2::GofAngle2(SimInfo* info, const std::string& filename, const std::string& sele1, 
+  GofAngle2::GofAngle2(SimInfo* info, const std::string& filename, 
+                       const std::string& sele1, 
 		       const std::string& sele2, int nangleBins)
-    : RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins) {
+    : RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins), 
+      evaluator3_(info), 
+      seleMan3_(info), doSele3_(false) {
+    
+    setOutputName(getPrefix(filename) + ".gto");
+    
+    deltaCosAngle_ = 2.0 / nAngleBins_;
+    
+    histogram_.resize(nAngleBins_);
+    avgGofr_.resize(nAngleBins_);
+    for (int i = 0 ; i < nAngleBins_; ++i) {
+      histogram_[i].resize(nAngleBins_);
+      avgGofr_[i].resize(nAngleBins_);
+    }   
+  }
 
-      setOutputName(getPrefix(filename) + ".gto");
+  GofAngle2::GofAngle2(SimInfo* info, const std::string& filename, 
+                       const std::string& sele1, 
+		       const std::string& sele2, 
+                       const std::string& sele3, int nangleBins)
+    : RadialDistrFunc(info, filename, sele1, sele2), nAngleBins_(nangleBins), 
+      evaluator3_(info), selectionScript3_(sele3),
+      seleMan3_(info), doSele3_(true) {
+    
+    setOutputName(getPrefix(filename) + ".gto");
 
-      deltaCosAngle_ = 2.0 / nAngleBins_;
+    deltaCosAngle_ = 2.0 / nAngleBins_;
+    
+    histogram_.resize(nAngleBins_);
+    avgGofr_.resize(nAngleBins_);
+    for (int i = 0 ; i < nAngleBins_; ++i) {
+      histogram_[i].resize(nAngleBins_);
+      avgGofr_[i].resize(nAngleBins_);
+    }    
+    evaluator3_.loadScriptString(sele3);      
+    if (!evaluator3_.isDynamic()) {
+      seleMan3_.setSelectionSet(evaluator3_.evaluate());
+    }
+  }
 
-      histogram_.resize(nAngleBins_);
-      avgGofr_.resize(nAngleBins_);
-      for (int i = 0 ; i < nAngleBins_; ++i) {
-        histogram_[i].resize(nAngleBins_);
-        avgGofr_[i].resize(nAngleBins_);
-      }    
+  void GofAngle2::processNonOverlapping( SelectionManager& sman1, 
+                                         SelectionManager& sman2) {
+    StuntDouble* sd1;
+    StuntDouble* sd2;
+    StuntDouble* sd3;
+    int i;    
+    int j;
+    int k;
+    
+    // This is the same as a non-overlapping pairwise loop structure:
+    // for (int i = 0;  i < ni ; ++i ) {
+    //   for (int j = 0; j < nj; ++j) {} 
+    // }
 
+    if (doSele3_) {
+      if  (evaluator3_.isDynamic()) {
+        seleMan3_.setSelectionSet(evaluator3_.evaluate());
+      }
+      if (sman1.getSelectionCount() != seleMan3_.getSelectionCount() ) {
+        RadialDistrFunc::processNonOverlapping( sman1, sman2 );
+      }
+
+      for (sd1 = sman1.beginSelected(i), sd3 = seleMan3_.beginSelected(k); 
+           sd1 != NULL && sd3 != NULL; 
+           sd1 = sman1.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
+        for (sd2 = sman2.beginSelected(j); sd2 != NULL; 
+             sd2 = sman2.nextSelected(j)) {
+          collectHistogram(sd1, sd2, sd3);
+        }
+      }
+    } else {
+      RadialDistrFunc::processNonOverlapping( sman1, sman2 );
     }
+  }
 
+  void GofAngle2::processOverlapping( SelectionManager& sman) {
+    StuntDouble* sd1;
+    StuntDouble* sd2;
+    StuntDouble* sd3;
+    int i;    
+    int j;
+    int k;
 
+    // This is the same as a pairwise loop structure:
+    // for (int i = 0;  i < n-1 ; ++i ) {
+    //   for (int j = i + 1; j < n; ++j) {} 
+    // }
+    
+    if (doSele3_) {
+      if  (evaluator3_.isDynamic()) {
+        seleMan3_.setSelectionSet(evaluator3_.evaluate());
+      }
+      if (sman.getSelectionCount() != seleMan3_.getSelectionCount() ) {
+        RadialDistrFunc::processOverlapping( sman);
+      }
+      for (sd1 = sman.beginSelected(i), sd3 = seleMan3_.beginSelected(k); 
+           sd1 != NULL && sd3 != NULL; 
+           sd1 = sman.nextSelected(i), sd3 = seleMan3_.nextSelected(k)) {
+        for (j  = i, sd2 = sman.nextSelected(j); sd2 != NULL; 
+             sd2 = sman.nextSelected(j)) {
+          collectHistogram(sd1, sd2, sd3);
+        }            
+      }
+    } else {
+      RadialDistrFunc::processOverlapping( sman);
+    }    
+  }
+
+
   void GofAngle2::preProcess() {
 
     for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
@@ -104,6 +198,22 @@ namespace OpenMD {
     MultipoleAdapter ma1 = MultipoleAdapter(atype1);
     MultipoleAdapter ma2 = MultipoleAdapter(atype2);
 
+    if (!sd1->isDirectional()) {
+      sprintf(painCave.errMsg, 
+              "GofAngle2: attempted to use a non-directional object: %s\n", 
+              sd1->getType().c_str());
+      painCave.isFatal = 1;
+      simError();  
+    }
+
+    if (!sd2->isDirectional()) {
+      sprintf(painCave.errMsg, 
+              "GofAngle2: attempted to use a non-directional object: %s\n", 
+              sd2->getType().c_str());
+      painCave.isFatal = 1;
+      simError();  
+    }
+
     Vector3d dipole1, dipole2;
     if (ma1.isDipole())         
         dipole1 = sd1->getDipole();
@@ -131,13 +241,72 @@ namespace OpenMD {
     ++npairs_;
   }
 
+  void GofAngle2::collectHistogram(StuntDouble* sd1, StuntDouble* sd2, 
+                                   StuntDouble* sd3) {
+
+    if (sd1 == sd2) {
+      return;
+    }
+
+    Vector3d p1 = sd1->getPos();
+    Vector3d p3 = sd3->getPos();
+
+    Vector3d c = 0.5 * (p1 + p3);
+    Vector3d r13 = p3 - p1;
+
+    Vector3d r12 = sd2->getPos() - c;
+  
+    if (usePeriodicBoundaryConditions_) {
+      currentSnapshot_->wrapVector(r12);
+      currentSnapshot_->wrapVector(r13);
+    }
+    r12.normalize();
+    r13.normalize();
+
+    if (!sd2->isDirectional()) {
+      sprintf(painCave.errMsg, 
+              "GofAngle2: attempted to use a non-directional object: %s\n", 
+              sd2->getType().c_str());
+      painCave.isFatal = 1;
+      simError();  
+    }
+
+    AtomType* atype2 = static_cast<Atom*>(sd2)->getAtomType();
+    MultipoleAdapter ma2 = MultipoleAdapter(atype2);
+
+    Vector3d dipole2;
+    if (ma2.isDipole())         
+        dipole2 = sd2->getDipole();
+    else
+        dipole2 = sd2->getA().transpose() * V3Z;
+    
+    dipole2.normalize();    
+
+    RealType cosAngle1 = dot(r12, r13);
+    RealType cosAngle2 = dot(r13, dipole2);
+
+    RealType halfBin = (nAngleBins_ - 1) * 0.5;
+    int angleBin1 = int(halfBin * (cosAngle1 + 1.0));
+    int angleBin2 = int(halfBin * (cosAngle2 + 1.0));
+
+    ++histogram_[angleBin1][angleBin2];    
+    ++npairs_;
+
+  }
+
   void GofAngle2::writeRdf() {
     std::ofstream rdfStream(outputFilename_.c_str());
     if (rdfStream.is_open()) {
       rdfStream << "#radial distribution function\n";
       rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
-      rdfStream << "selection2: (" << selectionScript2_ << ")\n";
-      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = " << deltaCosAngle_ << "\n";
+      rdfStream << "selection2: (" << selectionScript2_ << ")";
+      if (doSele3_) {
+        rdfStream << "\tselection3: (" << selectionScript3_ << ")\n";
+      } else {
+        rdfStream << "\n";
+      }
+      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = "
+                << deltaCosAngle_ << "\n";
       for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
 	// RealType cosAngle1 = -1.0 + (i + 0.5)*deltaCosAngle_;
         
@@ -145,13 +314,13 @@ namespace OpenMD {
 	  // RealType cosAngle2 = -1.0 + (j + 0.5)*deltaCosAngle_;
 	  rdfStream <<avgGofr_[i][j]/nProcessed_ << "\t";
 	}
-
 	rdfStream << "\n";
       }
         
     } else {
 
-      sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n", outputFilename_.c_str());
+      sprintf(painCave.errMsg, "GofAngle2: unable to open %s\n", 
+              outputFilename_.c_str());
       painCave.isFatal = 1;
       simError();  
     }