--- trunk/src/applications/staticProps/RNEMDStats.cpp 2013/06/18 16:07:27 1881 +++ trunk/src/applications/staticProps/RNEMDStats.cpp 2013/06/18 16:16:40 1883 @@ -96,8 +96,8 @@ namespace OpenMD { vector binMass(nBins_, 0.0); vector binVel(nBins_, V3Zero); vector binKE(nBins_, 0.0); - vector binDof(nBins_, 0); - vector binCount(nBins_, 0); + vector binDof(nBins_, 0); + vector binCount(nBins_, 0); for (mol = info_->beginMolecule(mi); mol != NULL; @@ -122,14 +122,14 @@ namespace OpenMD { // figure out where that object is: Vector3d pos = sd->getPos(); - currentSnapshot_->wrapVector(pos); + Vector3d vel = sd->getVel(); + RealType m = sd->getMass(); + currentSnapshot_->wrapVector(pos); int bin = getBin(pos); - binCount[bin]++; + binCount[bin] += 1; - RealType m = sd->getMass(); binMass[bin] += m; - Vector3d vel = sd->getVel(); binVel[bin] += vel; binKE[bin] += 0.5 * (m * vel.lengthSquare()); binDof[bin] += 3; @@ -153,16 +153,18 @@ namespace OpenMD { } } - for (int i = 0; i < nBins_; i++) { - RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb * - PhysicalConstants::energyConvert); - RealType den = binMass[i] * nBins_ * PhysicalConstants::densityConvert - / volume_; - Vector3d vel = binVel[i] / RealType(binCount[i]); - dynamic_cast(temperature->accumulator[i])->add(temp); - dynamic_cast(velocity->accumulator[i])->add(vel); - dynamic_cast(density->accumulator[i])->add(den); - dynamic_cast(counts_->accumulator[i])->add(1); + for (unsigned int i = 0; i < nBins_; i++) { + if (binDof[i] > 0) { + RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb * + PhysicalConstants::energyConvert); + RealType den = binMass[i] * nBins_ * PhysicalConstants::densityConvert + / volume_; + Vector3d vel = binVel[i] / RealType(binCount[i]); + dynamic_cast(temperature->accumulator[i])->add(temp); + dynamic_cast(velocity->accumulator[i])->add(vel); + dynamic_cast(density->accumulator[i])->add(den); + dynamic_cast(counts_->accumulator[i])->add(1); + } } }