ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/rnemd/RNEMD.cpp
(Generate patch)

Comparing trunk/src/integrators/RNEMD.cpp (file contents):
Revision 1338 by skuang, Thu Apr 23 18:22:30 2009 UTC vs.
Revision 1339 by gezelter, Thu Apr 23 18:31:05 2009 UTC

# Line 158 | Line 158 | namespace oopse {
158        // which bin is this stuntdouble in?
159        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
160  
161 <      int binNo = midBin + int(nBins_ * (pos.z()) / hmat(2,2));
161 >      int binNo = int((nBins_-1) * (pos.z() + 0.5*hmat(2,2)) / hmat(2,2));
162  
163        //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
164  
# Line 319 | Line 319 | namespace oopse {
319      bool max_found = false;
320      StuntDouble* max_sd;
321      */
322 <    std::vector<RealType> valueHist;
323 <    std::vector<int> valueCount;
322 >    std::vector<RealType> valueHist;  // keeps track of what's being averaged
323 >    std::vector<int> valueCount; // keeps track of the number of degrees of
324 >                                 // freedom being averaged
325      valueHist.resize(nBins_);
326      valueCount.resize(nBins_);
327      //do they initialize themselves to zero automatically?
# Line 339 | Line 340 | namespace oopse {
340        // which bin is this stuntdouble in?
341        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
342        
343 <      int binNo = midBin + int(nBins_ * (pos.z()) / hmat(2,2));
343 >      int binNo = int((nBins_-1) * (pos.z()+0.5*hmat(2,2)) / hmat(2,2));
344        
345        //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
346        
# Line 354 | Line 355 | namespace oopse {
355          value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
356                          vel[2]*vel[2]);
357          
358 +        valueCount[binNo] += 3;
359 +
360          if (sd->isDirectional()) {
361            Vector3d angMom = sd->getJ();
362            Mat3x3d I = sd->getI();
# Line 364 | Line 367 | namespace oopse {
367              int k = (i + 2) % 3;
368              value += angMom[j] * angMom[j] / I(j, j) +
369                angMom[k] * angMom[k] / I(k, k);
370 +
371 +            valueCount[binNo] +=2;
372 +
373            } else {                        
374              value += angMom[0]*angMom[0]/I(0, 0)
375                + angMom[1]*angMom[1]/I(1, 1)
376                + angMom[2]*angMom[2]/I(2, 2);
377 +            valueCount[binNo] +=3;
378 +
379            }
380          }
381          //std::cerr <<"this value = " << value << "\n";
382 <        value = value * 0.5 / OOPSEConstant::energyConvert;
382 >        value *= 0.5 / OOPSEConstant::energyConvert;  // get it in kcal / mol
383 >        value *= 2.0 / OOPSEConstant::kb;             // convert to temperature
384          //std::cerr <<"this value = " << value << "\n";
385          break;
386        case rnemdPx :
387          value = mass * vel[0];
388 +        valueCount[binNo]++;
389          break;
390        case rnemdPy :
391          value = mass * vel[1];
392 +        valueCount[binNo]++;
393          break;
394        case rnemdPz :
395          value = mass * vel[2];
396 +        valueCount[binNo]++;
397          break;
398        case rnemdUnknown :
399        default :
# Line 389 | Line 401 | namespace oopse {
401        }
402        //std::cerr << "bin = " << binNo << " value = " << value ;
403        valueHist[binNo] += value;
392 –      valueCount[binNo]++;
404        //std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n";
405      }
406      

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines