| 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 |
|
|
| 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? |
| 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 |
|
|
| 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(); |
| 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 : |
| 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 |
|
|
