feature: updating S(Q) and pair correlation functions

Author: agdelma

include/action.h

@@ -135,6 +135,7 @@ class ActionBase {
         beadLocator bead2,bead3;        // Bead indexers
         dVec sep,sep2;                  // The spatial separation between beads.
         double dSep;                    // The discretization for the separation histogram
+        double dPerSep;                 // The PBC discretization for the separation histogram 
 
         /* Update the separation histogram */
         void updateSepHist(const dVec &);   

src/action.cpp

@@ -43,7 +43,8 @@ ActionBase::ActionBase (const Path &_path, LookupTable &_lookup,
     sepHist = 0;
     cylSepHist.resize(NPCFSEP);
     cylSepHist = 0;
-    dSep = 0.5*sqrt(sum(path.boxPtr->periodic))*path.boxPtr->side[NDIM-1] / (1.0*NPCFSEP);
+    dSep = 0.5*sqrt(NDIM)*path.boxPtr->side[NDIM-1] / (1.0*NPCFSEP);
+    dPerSep = 0.5*sqrt(sum(path.boxPtr->periodic))*path.boxPtr->side[NDIM-1] / (1.0*NPCFSEP);
 
     /* Needed for canonical ensemble weighting */
     canonical = constants()->canonical();
@@ -70,9 +71,11 @@ ActionBase::~ActionBase() {
 ******************************************************************************/
 inline void ActionBase::updateSepHist(const dVec &_sep) {
     dVec psep;
-    psep = path.boxPtr->periodic*_sep;
+    /* Temporarly computing full radial g(r) */
+    /* psep = path.boxPtr->periodic*_sep; */
+    psep = _sep;
     int nR = int(sqrt(dot(psep,psep))/dSep);
-    if (nR < NPCFSEP)
+    if (nR >= 0 && nR < NPCFSEP)
         ++sepHist(nR);
 }
 
@@ -658,7 +661,7 @@ double LocalAction::V(const int slice, const double maxR) {
 
                 sep = path.getSeparation(bead2,bead1);
                 if (r2sq < maxR*maxR) {
-                    int nR = int(abs(sep[2])/dSep);
+                    int nR = int(abs(sep[2])/dPerSep);
                     if (nR < NPCFSEP)
                         ++cylSepHist(nR);
                 }

src/estimator.cpp

@@ -2574,7 +2574,7 @@ PairCorrelationEstimator::PairCorrelationEstimator (const Path &_path,
 
     /* The normalization factor for the pair correlation function depends 
      * on the dimensionality, and container type */
-    if (path.boxPtr->name == "Cylinder") {
+    if (path.boxPtr->name == "XXX") {
         for (int n = 0; n < NPCFSEP; n++)
             norm(n) = 0.5*path.boxPtr->side[NDIM-1] / dR;
     }