| 6 |
|
* redistribute this software in source and binary code form, provided |
| 7 |
|
* that the following conditions are met: |
| 8 |
|
* |
| 9 |
< |
* 1. Acknowledgement of the program authors must be made in any |
| 10 |
< |
* publication of scientific results based in part on use of the |
| 11 |
< |
* program. An acceptable form of acknowledgement is citation of |
| 12 |
< |
* the article in which the program was described (Matthew |
| 13 |
< |
* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
| 14 |
< |
* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
| 15 |
< |
* Parallel Simulation Engine for Molecular Dynamics," |
| 16 |
< |
* J. Comput. Chem. 26, pp. 252-271 (2005)) |
| 17 |
< |
* |
| 18 |
< |
* 2. Redistributions of source code must retain the above copyright |
| 9 |
> |
* 1. Redistributions of source code must retain the above copyright |
| 10 |
|
* notice, this list of conditions and the following disclaimer. |
| 11 |
|
* |
| 12 |
< |
* 3. Redistributions in binary form must reproduce the above copyright |
| 12 |
> |
* 2. Redistributions in binary form must reproduce the above copyright |
| 13 |
|
* notice, this list of conditions and the following disclaimer in the |
| 14 |
|
* documentation and/or other materials provided with the |
| 15 |
|
* distribution. |
| 28 |
|
* arising out of the use of or inability to use software, even if the |
| 29 |
|
* University of Notre Dame has been advised of the possibility of |
| 30 |
|
* such damages. |
| 31 |
+ |
* |
| 32 |
+ |
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
| 33 |
+ |
* research, please cite the appropriate papers when you publish your |
| 34 |
+ |
* work. Good starting points are: |
| 35 |
+ |
* |
| 36 |
+ |
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
| 37 |
+ |
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
| 38 |
+ |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
| 39 |
+ |
* [4] Vardeman & Gezelter, in progress (2009). |
| 40 |
|
*/ |
| 41 |
|
|
| 42 |
|
#include <cmath> |
| 46 |
|
#include "math/Polynomial.hpp" |
| 47 |
|
#include "primitives/Molecule.hpp" |
| 48 |
|
#include "primitives/StuntDouble.hpp" |
| 49 |
< |
#include "utils/OOPSEConstant.hpp" |
| 49 |
> |
#include "utils/PhysicalConstants.hpp" |
| 50 |
|
#include "utils/Tuple.hpp" |
| 51 |
|
|
| 52 |
|
#ifndef IS_MPI |
| 57 |
|
|
| 58 |
|
#define HONKING_LARGE_VALUE 1.0e10 |
| 59 |
|
|
| 60 |
< |
namespace oopse { |
| 60 |
> |
namespace OpenMD { |
| 61 |
|
|
| 62 |
|
RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { |
| 63 |
|
|
| 111 |
|
std::cerr << "WARNING! RNEMD Type Unknown!\n"; |
| 112 |
|
} |
| 113 |
|
|
| 114 |
+ |
output3DTemp_ = false; |
| 115 |
+ |
if (simParams->haveRNEMD_outputDimensionalTemperature()) { |
| 116 |
+ |
output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature(); |
| 117 |
+ |
} |
| 118 |
+ |
|
| 119 |
|
#ifdef IS_MPI |
| 120 |
|
if (worldRank == 0) { |
| 121 |
|
#endif |
| 122 |
|
|
| 123 |
|
std::string rnemdFileName; |
| 119 |
– |
std::string xTempFileName; |
| 120 |
– |
std::string yTempFileName; |
| 121 |
– |
std::string zTempFileName; |
| 124 |
|
switch(rnemdType_) { |
| 125 |
|
case rnemdKineticSwap : |
| 126 |
|
case rnemdKineticScale : |
| 131 |
|
case rnemdPy : |
| 132 |
|
case rnemdPyScale : |
| 133 |
|
rnemdFileName = "momemtum.log"; |
| 132 |
– |
xTempFileName = "temperatureX.log"; |
| 133 |
– |
yTempFileName = "temperatureY.log"; |
| 134 |
– |
zTempFileName = "temperatureZ.log"; |
| 135 |
– |
xTempLog_.open(xTempFileName.c_str()); |
| 136 |
– |
yTempLog_.open(yTempFileName.c_str()); |
| 137 |
– |
zTempLog_.open(zTempFileName.c_str()); |
| 134 |
|
break; |
| 135 |
|
case rnemdPz : |
| 136 |
|
case rnemdPzScale : |
| 141 |
|
} |
| 142 |
|
rnemdLog_.open(rnemdFileName.c_str()); |
| 143 |
|
|
| 144 |
+ |
std::string xTempFileName; |
| 145 |
+ |
std::string yTempFileName; |
| 146 |
+ |
std::string zTempFileName; |
| 147 |
+ |
if (output3DTemp_) { |
| 148 |
+ |
xTempFileName = "temperatureX.log"; |
| 149 |
+ |
yTempFileName = "temperatureY.log"; |
| 150 |
+ |
zTempFileName = "temperatureZ.log"; |
| 151 |
+ |
xTempLog_.open(xTempFileName.c_str()); |
| 152 |
+ |
yTempLog_.open(yTempFileName.c_str()); |
| 153 |
+ |
zTempLog_.open(zTempFileName.c_str()); |
| 154 |
+ |
} |
| 155 |
+ |
|
| 156 |
|
#ifdef IS_MPI |
| 157 |
|
} |
| 158 |
|
#endif |
| 160 |
|
set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime()); |
| 161 |
|
set_RNEMD_nBins(simParams->getRNEMD_nBins()); |
| 162 |
|
midBin_ = nBins_ / 2; |
| 163 |
+ |
if (simParams->haveRNEMD_binShift()) { |
| 164 |
+ |
if (simParams->getRNEMD_binShift()) { |
| 165 |
+ |
zShift_ = 0.5 / (RealType)(nBins_); |
| 166 |
+ |
} else { |
| 167 |
+ |
zShift_ = 0.0; |
| 168 |
+ |
} |
| 169 |
+ |
} else { |
| 170 |
+ |
zShift_ = 0.0; |
| 171 |
+ |
} |
| 172 |
+ |
//std::cerr << "we have zShift_ = " << zShift_ << "\n"; |
| 173 |
+ |
//shift slabs by half slab width, might be useful in heterogeneous systems |
| 174 |
+ |
//set to 0.0 if not using it; can NOT be used in status output yet |
| 175 |
|
if (simParams->haveRNEMD_logWidth()) { |
| 176 |
< |
rnemdLogWidth_ = simParams->getRNEMD_logWidth(); |
| 177 |
< |
if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) { |
| 176 |
> |
set_RNEMD_logWidth(simParams->getRNEMD_logWidth()); |
| 177 |
> |
/*arbitary rnemdLogWidth_ no checking |
| 178 |
> |
if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) { |
| 179 |
|
std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; |
| 180 |
|
std::cerr << "Automaically set back to default.\n"; |
| 181 |
|
rnemdLogWidth_ = nBins_; |
| 182 |
< |
} |
| 182 |
> |
}*/ |
| 183 |
|
} else { |
| 184 |
< |
rnemdLogWidth_ = nBins_; |
| 184 |
> |
set_RNEMD_logWidth(nBins_); |
| 185 |
|
} |
| 186 |
|
valueHist_.resize(rnemdLogWidth_, 0.0); |
| 187 |
|
valueCount_.resize(rnemdLogWidth_, 0); |
| 188 |
|
xTempHist_.resize(rnemdLogWidth_, 0.0); |
| 189 |
|
yTempHist_.resize(rnemdLogWidth_, 0.0); |
| 190 |
|
zTempHist_.resize(rnemdLogWidth_, 0.0); |
| 191 |
+ |
xyzTempCount_.resize(rnemdLogWidth_, 0); |
| 192 |
|
|
| 193 |
|
set_RNEMD_exchange_total(0.0); |
| 194 |
|
if (simParams->haveRNEMD_targetFlux()) { |
| 216 |
|
|
| 217 |
|
RNEMD::~RNEMD() { |
| 218 |
|
delete randNumGen_; |
| 219 |
< |
|
| 198 |
< |
std::cerr << "total fail trials: " << failTrialCount_ << "\n"; |
| 219 |
> |
|
| 220 |
|
#ifdef IS_MPI |
| 221 |
|
if (worldRank == 0) { |
| 222 |
|
#endif |
| 223 |
+ |
std::cerr << "total fail trials: " << failTrialCount_ << "\n"; |
| 224 |
|
rnemdLog_.close(); |
| 225 |
|
if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) |
| 226 |
|
std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n"; |
| 227 |
< |
if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) { |
| 227 |
> |
if (output3DTemp_) { |
| 228 |
|
xTempLog_.close(); |
| 229 |
|
yTempLog_.close(); |
| 230 |
|
zTempLog_.close(); |
| 268 |
|
// which bin is this stuntdouble in? |
| 269 |
|
// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
| 270 |
|
|
| 271 |
< |
int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
| 271 |
> |
int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_; |
| 272 |
|
|
| 273 |
|
|
| 274 |
|
// if we're in bin 0 or the middleBin |
| 283 |
|
|
| 284 |
|
value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + |
| 285 |
|
vel[2]*vel[2]); |
| 286 |
+ |
/* |
| 287 |
|
if (sd->isDirectional()) { |
| 288 |
|
Vector3d angMom = sd->getJ(); |
| 289 |
|
Mat3x3d I = sd->getI(); |
| 299 |
|
+ angMom[1]*angMom[1]/I(1, 1) |
| 300 |
|
+ angMom[2]*angMom[2]/I(2, 2); |
| 301 |
|
} |
| 302 |
< |
} |
| 302 |
> |
} no exchange of angular momenta |
| 303 |
> |
*/ |
| 304 |
|
//make exchangeSum_ comparable between swap & scale |
| 305 |
|
//temporarily without using energyConvert |
| 306 |
< |
//value = value * 0.5 / OOPSEConstant::energyConvert; |
| 306 |
> |
//value = value * 0.5 / PhysicalConstants::energyConvert; |
| 307 |
|
value *= 0.5; |
| 308 |
|
break; |
| 309 |
|
case rnemdPx : |
| 354 |
|
bool my_min_found = min_found; |
| 355 |
|
bool my_max_found = max_found; |
| 356 |
|
|
| 357 |
< |
// Even if we didn't find a minimum, did someone else? |
| 358 |
< |
MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, |
| 359 |
< |
1, MPI::BOOL, MPI::LAND); |
| 360 |
< |
|
| 357 |
> |
// Even if we didn't find a minimum, did someone else? debugging... |
| 358 |
> |
//MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, |
| 359 |
> |
// 1, MPI::BOOL, MPI::LAND); |
| 360 |
> |
MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR); |
| 361 |
|
// Even if we didn't find a maximum, did someone else? |
| 362 |
< |
MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, |
| 363 |
< |
1, MPI::BOOL, MPI::LAND); |
| 364 |
< |
|
| 362 |
> |
//MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, |
| 363 |
> |
// 1, MPI::BOOL, MPI::LAND); |
| 364 |
> |
MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR); |
| 365 |
|
struct { |
| 366 |
|
RealType val; |
| 367 |
|
int rank; |
| 397 |
|
#endif |
| 398 |
|
|
| 399 |
|
if (max_found && min_found) { |
| 400 |
< |
if (min_val< max_val) { |
| 400 |
> |
if (min_val < max_val) { |
| 401 |
|
|
| 402 |
|
#ifdef IS_MPI |
| 403 |
|
if (max_vals.rank == worldRank && min_vals.rank == worldRank) { |
| 404 |
|
// I have both maximum and minimum, so proceed like a single |
| 405 |
|
// processor version: |
| 406 |
|
#endif |
| 407 |
< |
// objects to be swapped: velocity & angular velocity |
| 407 |
> |
// objects to be swapped: velocity ONLY |
| 408 |
|
Vector3d min_vel = min_sd->getVel(); |
| 409 |
|
Vector3d max_vel = max_sd->getVel(); |
| 410 |
|
RealType temp_vel; |
| 413 |
|
case rnemdKineticSwap : |
| 414 |
|
min_sd->setVel(max_vel); |
| 415 |
|
max_sd->setVel(min_vel); |
| 416 |
+ |
/* |
| 417 |
|
if (min_sd->isDirectional() && max_sd->isDirectional()) { |
| 418 |
|
Vector3d min_angMom = min_sd->getJ(); |
| 419 |
|
Vector3d max_angMom = max_sd->getJ(); |
| 420 |
|
min_sd->setJ(max_angMom); |
| 421 |
|
max_sd->setJ(min_angMom); |
| 422 |
< |
} |
| 422 |
> |
} no angular momentum exchange |
| 423 |
> |
*/ |
| 424 |
|
break; |
| 425 |
|
case rnemdPx : |
| 426 |
|
temp_vel = min_vel.x(); |
| 464 |
|
switch(rnemdType_) { |
| 465 |
|
case rnemdKineticSwap : |
| 466 |
|
max_sd->setVel(min_vel); |
| 467 |
< |
|
| 467 |
> |
/* |
| 468 |
|
if (max_sd->isDirectional()) { |
| 469 |
|
Vector3d min_angMom; |
| 470 |
|
Vector3d max_angMom = max_sd->getJ(); |
| 477 |
|
status); |
| 478 |
|
|
| 479 |
|
max_sd->setJ(min_angMom); |
| 480 |
< |
} |
| 480 |
> |
} no angular momentum exchange |
| 481 |
> |
*/ |
| 482 |
|
break; |
| 483 |
|
case rnemdPx : |
| 484 |
|
max_vel.x() = min_vel.x(); |
| 511 |
|
switch(rnemdType_) { |
| 512 |
|
case rnemdKineticSwap : |
| 513 |
|
min_sd->setVel(max_vel); |
| 514 |
< |
|
| 514 |
> |
/* |
| 515 |
|
if (min_sd->isDirectional()) { |
| 516 |
|
Vector3d min_angMom = min_sd->getJ(); |
| 517 |
|
Vector3d max_angMom; |
| 524 |
|
status); |
| 525 |
|
|
| 526 |
|
min_sd->setJ(max_angMom); |
| 527 |
< |
} |
| 527 |
> |
} no angular momentum exchange |
| 528 |
> |
*/ |
| 529 |
|
break; |
| 530 |
|
case rnemdPx : |
| 531 |
|
min_vel.x() = max_vel.x(); |
| 598 |
|
// which bin is this stuntdouble in? |
| 599 |
|
// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
| 600 |
|
|
| 601 |
< |
int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
| 601 |
> |
int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_; |
| 602 |
|
|
| 603 |
|
// if we're in bin 0 or the middleBin |
| 604 |
|
if (binNo == 0 || binNo == midBin_) { |
| 670 |
|
a001 = Kcx * px * px; |
| 671 |
|
*/ |
| 672 |
|
|
| 673 |
< |
//scale all three dimensions, let x = y |
| 673 |
> |
//scale all three dimensions, let c_x = c_y |
| 674 |
|
a000 = Kcx + Kcy; |
| 675 |
|
a110 = Kcz; |
| 676 |
|
c0 = targetFlux_ - Kcx - Kcy - Kcz; |
| 706 |
|
+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); |
| 707 |
|
break; |
| 708 |
|
case rnemdPzScale ://we don't really do this, do we? |
| 709 |
+ |
c = 1 - targetFlux_ / Pcz; |
| 710 |
+ |
a000 = Kcx; |
| 711 |
+ |
a110 = Kcy; |
| 712 |
+ |
c0 = Kcz * c * c - Kcx - Kcy - Kcz; |
| 713 |
+ |
a001 = px * px * Khx; |
| 714 |
+ |
a111 = py * py * Khy; |
| 715 |
+ |
b01 = -2.0 * Khx * px * (1.0 + px); |
| 716 |
+ |
b11 = -2.0 * Khy * py * (1.0 + py); |
| 717 |
+ |
c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) |
| 718 |
+ |
+ Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0); |
| 719 |
+ |
break; |
| 720 |
|
default : |
| 721 |
|
break; |
| 722 |
|
} |
| 760 |
|
r2 = *ri; |
| 761 |
|
//check if FindRealRoots() give the right answer |
| 762 |
|
if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) { |
| 763 |
< |
std::cerr << "WARNING! eq solvers might have mistakes!\n"; |
| 763 |
> |
sprintf(painCave.errMsg, |
| 764 |
> |
"RNEMD Warning: polynomial solve seems to have an error!"); |
| 765 |
> |
painCave.isFatal = 0; |
| 766 |
> |
simError(); |
| 767 |
|
failRootCount_++; |
| 768 |
|
} |
| 769 |
|
//might not be useful w/o rescaling coefficients |
| 804 |
|
+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2); |
| 805 |
|
break; |
| 806 |
|
case rnemdPzScale : |
| 807 |
+ |
diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) |
| 808 |
+ |
+ fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2); |
| 809 |
|
default : |
| 810 |
|
break; |
| 811 |
|
} |
| 841 |
|
z = bestPair.second; |
| 842 |
|
break; |
| 843 |
|
case rnemdPzScale : |
| 844 |
+ |
x = bestPair.first; |
| 845 |
+ |
y = bestPair.second; |
| 846 |
+ |
z = c; |
| 847 |
+ |
break; |
| 848 |
|
default : |
| 849 |
|
break; |
| 850 |
|
} |
| 871 |
|
exchangeSum_ += targetFlux_; |
| 872 |
|
//we may want to check whether the exchange has been successful |
| 873 |
|
} else { |
| 874 |
< |
std::cerr << "exchange NOT performed!\n"; |
| 874 |
> |
std::cerr << "exchange NOT performed!\n";//MPI incompatible |
| 875 |
|
failTrialCount_++; |
| 876 |
|
} |
| 877 |
|
|
| 909 |
|
StuntDouble* sd; |
| 910 |
|
int idx; |
| 911 |
|
|
| 912 |
+ |
// alternative approach, track all molecules instead of only those selected for scaling/swapping |
| 913 |
+ |
//SimInfo::MoleculeIterator miter; |
| 914 |
+ |
//std::vector<StuntDouble*>::iterator iiter; |
| 915 |
+ |
//Molecule* mol; |
| 916 |
+ |
//StuntDouble* integrableObject; |
| 917 |
+ |
//for (mol = info_->beginMolecule(miter); mol != NULL; |
| 918 |
+ |
// mol = info_->nextMolecule(miter)) |
| 919 |
+ |
// integrableObject is essentially sd |
| 920 |
+ |
//for (integrableObject = mol->beginIntegrableObject(iiter); |
| 921 |
+ |
// integrableObject != NULL; |
| 922 |
+ |
// integrableObject = mol->nextIntegrableObject(iiter)) |
| 923 |
|
for (sd = seleMan_.beginSelected(selei); sd != NULL; |
| 924 |
|
sd = seleMan_.nextSelected(selei)) { |
| 925 |
|
|
| 935 |
|
// which bin is this stuntdouble in? |
| 936 |
|
// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] |
| 937 |
|
|
| 938 |
< |
int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; |
| 939 |
< |
|
| 938 |
> |
int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) % |
| 939 |
> |
rnemdLogWidth_; |
| 940 |
> |
/* no symmetrization allowed due to arbitary rnemdLogWidth_ value |
| 941 |
|
if (rnemdLogWidth_ == midBin_ + 1) |
| 942 |
|
if (binNo > midBin_) |
| 943 |
|
binNo = nBins_ - binNo; |
| 944 |
< |
|
| 944 |
> |
*/ |
| 945 |
|
RealType mass = sd->getMass(); |
| 946 |
|
Vector3d vel = sd->getVel(); |
| 947 |
|
RealType value; |
| 975 |
|
valueCount_[binNo] +=3; |
| 976 |
|
} |
| 977 |
|
} |
| 978 |
< |
value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb; |
| 978 |
> |
value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb; |
| 979 |
|
|
| 980 |
|
break; |
| 981 |
|
case rnemdPx : |
| 982 |
|
case rnemdPxScale : |
| 983 |
|
value = mass * vel[0]; |
| 984 |
|
valueCount_[binNo]++; |
| 925 |
– |
xVal = mass * vel.x() * vel.x() / OOPSEConstant::energyConvert |
| 926 |
– |
/ OOPSEConstant::kb; |
| 927 |
– |
yVal = mass * vel.y() * vel.y() / OOPSEConstant::energyConvert |
| 928 |
– |
/ OOPSEConstant::kb; |
| 929 |
– |
zVal = mass * vel.z() * vel.z() / OOPSEConstant::energyConvert |
| 930 |
– |
/ OOPSEConstant::kb; |
| 931 |
– |
xTempHist_[binNo] += xVal; |
| 932 |
– |
yTempHist_[binNo] += yVal; |
| 933 |
– |
zTempHist_[binNo] += zVal; |
| 985 |
|
break; |
| 986 |
|
case rnemdPy : |
| 987 |
|
case rnemdPyScale : |
| 990 |
|
break; |
| 991 |
|
case rnemdPz : |
| 992 |
|
case rnemdPzScale : |
| 993 |
< |
value = mass * vel[2]; |
| 993 |
> |
value = pos.z(); //temporarily for homogeneous systems ONLY |
| 994 |
|
valueCount_[binNo]++; |
| 995 |
|
break; |
| 996 |
|
case rnemdUnknown : |
| 997 |
|
default : |
| 998 |
+ |
value = 1.0; |
| 999 |
+ |
valueCount_[binNo]++; |
| 1000 |
|
break; |
| 1001 |
|
} |
| 1002 |
|
valueHist_[binNo] += value; |
| 950 |
– |
} |
| 1003 |
|
|
| 1004 |
+ |
if (output3DTemp_) { |
| 1005 |
+ |
xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert |
| 1006 |
+ |
/ PhysicalConstants::kb; |
| 1007 |
+ |
yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert |
| 1008 |
+ |
/ PhysicalConstants::kb; |
| 1009 |
+ |
zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert |
| 1010 |
+ |
/ PhysicalConstants::kb; |
| 1011 |
+ |
xTempHist_[binNo] += xVal; |
| 1012 |
+ |
yTempHist_[binNo] += yVal; |
| 1013 |
+ |
zTempHist_[binNo] += zVal; |
| 1014 |
+ |
xyzTempCount_[binNo]++; |
| 1015 |
+ |
} |
| 1016 |
+ |
} |
| 1017 |
|
} |
| 1018 |
|
|
| 1019 |
|
void RNEMD::getStarted() { |
| 1020 |
+ |
collectData(); |
| 1021 |
+ |
/* now should be able to output profile in step 0, but might not be useful |
| 1022 |
|
Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 1023 |
|
Stats& stat = currentSnap_->statData; |
| 1024 |
|
stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; |
| 1025 |
+ |
*/ |
| 1026 |
+ |
getStatus(); |
| 1027 |
|
} |
| 1028 |
|
|
| 1029 |
|
void RNEMD::getStatus() { |
| 1034 |
|
|
| 1035 |
|
stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; |
| 1036 |
|
//or to be more meaningful, define another item as exchangeSum_ / time |
| 1037 |
+ |
int j; |
| 1038 |
|
|
| 969 |
– |
|
| 1039 |
|
#ifdef IS_MPI |
| 1040 |
|
|
| 1041 |
|
// all processors have the same number of bins, and STL vectors pack their |
| 1045 |
|
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 1046 |
|
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0], |
| 1047 |
|
rnemdLogWidth_, MPI::INT, MPI::SUM); |
| 1048 |
< |
|
| 1048 |
> |
if (output3DTemp_) { |
| 1049 |
> |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0], |
| 1050 |
> |
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 1051 |
> |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0], |
| 1052 |
> |
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 1053 |
> |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0], |
| 1054 |
> |
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 1055 |
> |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0], |
| 1056 |
> |
rnemdLogWidth_, MPI::INT, MPI::SUM); |
| 1057 |
> |
} |
| 1058 |
|
// If we're the root node, should we print out the results |
| 1059 |
|
int worldRank = MPI::COMM_WORLD.Get_rank(); |
| 1060 |
|
if (worldRank == 0) { |
| 1061 |
|
#endif |
| 984 |
– |
int j; |
| 1062 |
|
rnemdLog_ << time; |
| 1063 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1064 |
|
rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j]; |
| 988 |
– |
valueHist_[j] = 0.0; |
| 1065 |
|
} |
| 1066 |
|
rnemdLog_ << "\n"; |
| 1067 |
< |
if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) { |
| 1067 |
> |
if (output3DTemp_) { |
| 1068 |
|
xTempLog_ << time; |
| 1069 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1070 |
< |
xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j]; |
| 995 |
< |
xTempHist_[j] = 0.0; |
| 1070 |
> |
xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j]; |
| 1071 |
|
} |
| 1072 |
|
xTempLog_ << "\n"; |
| 1073 |
|
yTempLog_ << time; |
| 1074 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1075 |
< |
yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j]; |
| 1001 |
< |
yTempHist_[j] = 0.0; |
| 1075 |
> |
yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j]; |
| 1076 |
|
} |
| 1077 |
|
yTempLog_ << "\n"; |
| 1078 |
|
zTempLog_ << time; |
| 1079 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1080 |
< |
zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j]; |
| 1007 |
< |
zTempHist_[j] = 0.0; |
| 1080 |
> |
zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j]; |
| 1081 |
|
} |
| 1082 |
|
zTempLog_ << "\n"; |
| 1083 |
|
} |
| 1011 |
– |
for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0; |
| 1084 |
|
#ifdef IS_MPI |
| 1085 |
|
} |
| 1086 |
|
#endif |
| 1087 |
+ |
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1088 |
+ |
valueCount_[j] = 0; |
| 1089 |
+ |
valueHist_[j] = 0.0; |
| 1090 |
+ |
} |
| 1091 |
+ |
if (output3DTemp_) |
| 1092 |
+ |
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1093 |
+ |
xTempHist_[j] = 0.0; |
| 1094 |
+ |
yTempHist_[j] = 0.0; |
| 1095 |
+ |
zTempHist_[j] = 0.0; |
| 1096 |
+ |
xyzTempCount_[j] = 0; |
| 1097 |
+ |
} |
| 1098 |
|
} |
| 1016 |
– |
|
| 1099 |
|
} |