| 57 |
|
|
| 58 |
|
#define HONKING_LARGE_VALUE 1.0e10 |
| 59 |
|
|
| 60 |
+ |
using namespace std; |
| 61 |
|
namespace OpenMD { |
| 62 |
|
|
| 63 |
< |
RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { |
| 63 |
> |
RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), |
| 64 |
> |
usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { |
| 65 |
|
|
| 66 |
|
failTrialCount_ = 0; |
| 67 |
|
failRootCount_ = 0; |
| 90 |
|
|
| 91 |
|
if (selectionCount > nIntegrable) { |
| 92 |
|
sprintf(painCave.errMsg, |
| 93 |
< |
"RNEMD warning: The current RNEMD_objectSelection,\n" |
| 93 |
> |
"RNEMD: The current RNEMD_objectSelection,\n" |
| 94 |
|
"\t\t%s\n" |
| 95 |
|
"\thas resulted in %d selected objects. However,\n" |
| 96 |
|
"\tthe total number of integrable objects in the system\n" |
| 100 |
|
rnemdObjectSelection_.c_str(), |
| 101 |
|
selectionCount, nIntegrable); |
| 102 |
|
painCave.isFatal = 0; |
| 103 |
+ |
painCave.severity = OPENMD_WARNING; |
| 104 |
|
simError(); |
| 102 |
– |
|
| 105 |
|
} |
| 106 |
|
|
| 107 |
< |
const std::string st = simParams->getRNEMD_exchangeType(); |
| 107 |
> |
const string st = simParams->getRNEMD_exchangeType(); |
| 108 |
|
|
| 109 |
< |
std::map<std::string, RNEMDTypeEnum>::iterator i; |
| 109 |
> |
map<string, RNEMDTypeEnum>::iterator i; |
| 110 |
|
i = stringToEnumMap_.find(st); |
| 111 |
|
rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second; |
| 112 |
|
if (rnemdType_ == rnemdUnknown) { |
| 113 |
< |
std::cerr << "WARNING! RNEMD Type Unknown!\n"; |
| 113 |
> |
sprintf(painCave.errMsg, |
| 114 |
> |
"RNEMD: The current RNEMD_exchangeType,\n" |
| 115 |
> |
"\t\t%s\n" |
| 116 |
> |
"\tis not one of the recognized exchange types.\n", |
| 117 |
> |
st.c_str()); |
| 118 |
> |
painCave.isFatal = 1; |
| 119 |
> |
painCave.severity = OPENMD_ERROR; |
| 120 |
> |
simError(); |
| 121 |
|
} |
| 122 |
< |
|
| 122 |
> |
|
| 123 |
|
output3DTemp_ = false; |
| 124 |
|
if (simParams->haveRNEMD_outputDimensionalTemperature()) { |
| 125 |
|
output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature(); |
| 129 |
|
if (worldRank == 0) { |
| 130 |
|
#endif |
| 131 |
|
|
| 132 |
< |
std::string rnemdFileName; |
| 132 |
> |
string rnemdFileName; |
| 133 |
|
switch(rnemdType_) { |
| 134 |
|
case rnemdKineticSwap : |
| 135 |
|
case rnemdKineticScale : |
| 150 |
|
} |
| 151 |
|
rnemdLog_.open(rnemdFileName.c_str()); |
| 152 |
|
|
| 153 |
< |
std::string xTempFileName; |
| 154 |
< |
std::string yTempFileName; |
| 155 |
< |
std::string zTempFileName; |
| 153 |
> |
string xTempFileName; |
| 154 |
> |
string yTempFileName; |
| 155 |
> |
string zTempFileName; |
| 156 |
|
if (output3DTemp_) { |
| 157 |
|
xTempFileName = "temperatureX.log"; |
| 158 |
|
yTempFileName = "temperatureY.log"; |
| 170 |
|
set_RNEMD_nBins(simParams->getRNEMD_nBins()); |
| 171 |
|
midBin_ = nBins_ / 2; |
| 172 |
|
if (simParams->haveRNEMD_binShift()) { |
| 173 |
< |
if (simParams->getRNEMD_binShift()) { |
| 174 |
< |
zShift_ = 0.5 / (RealType)(nBins_); |
| 175 |
< |
} else { |
| 176 |
< |
zShift_ = 0.0; |
| 177 |
< |
} |
| 173 |
> |
if (simParams->getRNEMD_binShift()) { |
| 174 |
> |
zShift_ = 0.5 / (RealType)(nBins_); |
| 175 |
> |
} else { |
| 176 |
> |
zShift_ = 0.0; |
| 177 |
> |
} |
| 178 |
|
} else { |
| 179 |
|
zShift_ = 0.0; |
| 180 |
|
} |
| 181 |
< |
//std::cerr << "we have zShift_ = " << zShift_ << "\n"; |
| 181 |
> |
//cerr << "we have zShift_ = " << zShift_ << "\n"; |
| 182 |
|
//shift slabs by half slab width, might be useful in heterogeneous systems |
| 183 |
|
//set to 0.0 if not using it; can NOT be used in status output yet |
| 184 |
|
if (simParams->haveRNEMD_logWidth()) { |
| 185 |
|
set_RNEMD_logWidth(simParams->getRNEMD_logWidth()); |
| 186 |
|
/*arbitary rnemdLogWidth_ no checking |
| 187 |
< |
if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) { |
| 188 |
< |
std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; |
| 189 |
< |
std::cerr << "Automaically set back to default.\n"; |
| 187 |
> |
if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) { |
| 188 |
> |
cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; |
| 189 |
> |
cerr << "Automaically set back to default.\n"; |
| 190 |
|
rnemdLogWidth_ = nBins_; |
| 191 |
|
}*/ |
| 192 |
|
} else { |
| 229 |
|
#ifdef IS_MPI |
| 230 |
|
if (worldRank == 0) { |
| 231 |
|
#endif |
| 232 |
< |
std::cerr << "total fail trials: " << failTrialCount_ << "\n"; |
| 232 |
> |
|
| 233 |
> |
sprintf(painCave.errMsg, |
| 234 |
> |
"RNEMD: total failed trials: %d\n", |
| 235 |
> |
failTrialCount_); |
| 236 |
> |
painCave.isFatal = 0; |
| 237 |
> |
painCave.severity = OPENMD_INFO; |
| 238 |
> |
simError(); |
| 239 |
> |
|
| 240 |
|
rnemdLog_.close(); |
| 241 |
< |
if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) |
| 242 |
< |
std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n"; |
| 241 |
> |
if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) { |
| 242 |
> |
sprintf(painCave.errMsg, |
| 243 |
> |
"RNEMD: total root-checking warnings: %d\n", |
| 244 |
> |
failRootCount_); |
| 245 |
> |
painCave.isFatal = 0; |
| 246 |
> |
painCave.severity = OPENMD_INFO; |
| 247 |
> |
simError(); |
| 248 |
> |
} |
| 249 |
|
if (output3DTemp_) { |
| 250 |
|
xTempLog_.close(); |
| 251 |
|
yTempLog_.close(); |
| 306 |
|
value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + |
| 307 |
|
vel[2]*vel[2]); |
| 308 |
|
/* |
| 309 |
< |
if (sd->isDirectional()) { |
| 309 |
> |
if (sd->isDirectional()) { |
| 310 |
|
Vector3d angMom = sd->getJ(); |
| 311 |
|
Mat3x3d I = sd->getI(); |
| 312 |
|
|
| 313 |
|
if (sd->isLinear()) { |
| 314 |
< |
int i = sd->linearAxis(); |
| 315 |
< |
int j = (i + 1) % 3; |
| 316 |
< |
int k = (i + 2) % 3; |
| 317 |
< |
value += angMom[j] * angMom[j] / I(j, j) + |
| 318 |
< |
angMom[k] * angMom[k] / I(k, k); |
| 314 |
> |
int i = sd->linearAxis(); |
| 315 |
> |
int j = (i + 1) % 3; |
| 316 |
> |
int k = (i + 2) % 3; |
| 317 |
> |
value += angMom[j] * angMom[j] / I(j, j) + |
| 318 |
> |
angMom[k] * angMom[k] / I(k, k); |
| 319 |
|
} else { |
| 320 |
< |
value += angMom[0]*angMom[0]/I(0, 0) |
| 321 |
< |
+ angMom[1]*angMom[1]/I(1, 1) |
| 322 |
< |
+ angMom[2]*angMom[2]/I(2, 2); |
| 320 |
> |
value += angMom[0]*angMom[0]/I(0, 0) |
| 321 |
> |
+ angMom[1]*angMom[1]/I(1, 1) |
| 322 |
> |
+ angMom[2]*angMom[2]/I(2, 2); |
| 323 |
|
} |
| 324 |
< |
} no exchange of angular momenta |
| 324 |
> |
} no exchange of angular momenta |
| 325 |
|
*/ |
| 326 |
|
//make exchangeSum_ comparable between swap & scale |
| 327 |
|
//temporarily without using energyConvert |
| 376 |
|
bool my_min_found = min_found; |
| 377 |
|
bool my_max_found = max_found; |
| 378 |
|
|
| 379 |
< |
// 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); |
| 379 |
> |
// Even if we didn't find a minimum, did someone else? |
| 380 |
|
MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR); |
| 381 |
|
// 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); |
| 382 |
|
MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR); |
| 383 |
|
struct { |
| 384 |
|
RealType val; |
| 432 |
|
min_sd->setVel(max_vel); |
| 433 |
|
max_sd->setVel(min_vel); |
| 434 |
|
/* |
| 435 |
< |
if (min_sd->isDirectional() && max_sd->isDirectional()) { |
| 435 |
> |
if (min_sd->isDirectional() && max_sd->isDirectional()) { |
| 436 |
|
Vector3d min_angMom = min_sd->getJ(); |
| 437 |
|
Vector3d max_angMom = max_sd->getJ(); |
| 438 |
|
min_sd->setJ(max_angMom); |
| 439 |
|
max_sd->setJ(min_angMom); |
| 440 |
< |
} no angular momentum exchange |
| 440 |
> |
} no angular momentum exchange |
| 441 |
|
*/ |
| 442 |
|
break; |
| 443 |
|
case rnemdPx : |
| 482 |
|
switch(rnemdType_) { |
| 483 |
|
case rnemdKineticSwap : |
| 484 |
|
max_sd->setVel(min_vel); |
| 485 |
< |
/* |
| 485 |
> |
//no angular momentum exchange for now |
| 486 |
> |
/* |
| 487 |
|
if (max_sd->isDirectional()) { |
| 488 |
|
Vector3d min_angMom; |
| 489 |
|
Vector3d max_angMom = max_sd->getJ(); |
| 490 |
< |
|
| 490 |
> |
|
| 491 |
|
// point-to-point swap of the angular momentum vector |
| 492 |
|
MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3, |
| 493 |
|
MPI::REALTYPE, min_vals.rank, 1, |
| 494 |
|
min_angMom.getArrayPointer(), 3, |
| 495 |
|
MPI::REALTYPE, min_vals.rank, 1, |
| 496 |
|
status); |
| 497 |
< |
|
| 497 |
> |
|
| 498 |
|
max_sd->setJ(min_angMom); |
| 499 |
< |
} no angular momentum exchange |
| 500 |
< |
*/ |
| 499 |
> |
} |
| 500 |
> |
*/ |
| 501 |
|
break; |
| 502 |
|
case rnemdPx : |
| 503 |
|
max_vel.x() = min_vel.x(); |
| 530 |
|
switch(rnemdType_) { |
| 531 |
|
case rnemdKineticSwap : |
| 532 |
|
min_sd->setVel(max_vel); |
| 533 |
< |
/* |
| 533 |
> |
// no angular momentum exchange for now |
| 534 |
> |
/* |
| 535 |
|
if (min_sd->isDirectional()) { |
| 536 |
|
Vector3d min_angMom = min_sd->getJ(); |
| 537 |
|
Vector3d max_angMom; |
| 538 |
< |
|
| 538 |
> |
|
| 539 |
|
// point-to-point swap of the angular momentum vector |
| 540 |
|
MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3, |
| 541 |
|
MPI::REALTYPE, max_vals.rank, 1, |
| 542 |
|
max_angMom.getArrayPointer(), 3, |
| 543 |
|
MPI::REALTYPE, max_vals.rank, 1, |
| 544 |
|
status); |
| 545 |
< |
|
| 545 |
> |
|
| 546 |
|
min_sd->setJ(max_angMom); |
| 547 |
< |
} no angular momentum exchange |
| 547 |
> |
} |
| 548 |
|
*/ |
| 549 |
|
break; |
| 550 |
|
case rnemdPx : |
| 565 |
|
} |
| 566 |
|
#endif |
| 567 |
|
exchangeSum_ += max_val - min_val; |
| 568 |
< |
} else { |
| 569 |
< |
std::cerr << "exchange NOT performed!\nmin_val > max_val.\n"; |
| 568 |
> |
} else { |
| 569 |
> |
sprintf(painCave.errMsg, |
| 570 |
> |
"RNEMD: exchange NOT performed because min_val > max_val\n"); |
| 571 |
> |
painCave.isFatal = 0; |
| 572 |
> |
painCave.severity = OPENMD_INFO; |
| 573 |
> |
simError(); |
| 574 |
|
failTrialCount_++; |
| 575 |
|
} |
| 576 |
|
} else { |
| 577 |
< |
std::cerr << "exchange NOT performed!\n"; |
| 578 |
< |
std::cerr << "at least one of the two slabs empty.\n"; |
| 577 |
> |
sprintf(painCave.errMsg, |
| 578 |
> |
"RNEMD: exchange NOT performed because at least one\n" |
| 579 |
> |
"\tof the two slabs is empty\n"); |
| 580 |
> |
painCave.isFatal = 0; |
| 581 |
> |
painCave.severity = OPENMD_INFO; |
| 582 |
> |
simError(); |
| 583 |
|
failTrialCount_++; |
| 584 |
|
} |
| 585 |
|
|
| 596 |
|
StuntDouble* sd; |
| 597 |
|
int idx; |
| 598 |
|
|
| 599 |
< |
std::vector<StuntDouble*> hotBin, coldBin; |
| 599 |
> |
vector<StuntDouble*> hotBin, coldBin; |
| 600 |
|
|
| 601 |
|
RealType Phx = 0.0; |
| 602 |
|
RealType Phy = 0.0; |
| 685 |
|
RealType a000, a110, c0, a001, a111, b01, b11, c1, c; |
| 686 |
|
switch(rnemdType_) { |
| 687 |
|
case rnemdKineticScale : |
| 688 |
< |
/*used hotBin coeff's & only scale x & y dimensions |
| 688 |
> |
// used hotBin coeff's & only scale x & y dimensions |
| 689 |
> |
/* |
| 690 |
|
RealType px = Phx / Pcx; |
| 691 |
|
RealType py = Phy / Pcy; |
| 692 |
|
a110 = Khy; |
| 693 |
|
c0 = - Khx - Khy - targetFlux_; |
| 694 |
|
a000 = Khx; |
| 695 |
< |
a111 = Kcy * py * py |
| 695 |
> |
a111 = Kcy * py * py; |
| 696 |
|
b11 = -2.0 * Kcy * py * (1.0 + py); |
| 697 |
|
c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_; |
| 698 |
|
b01 = -2.0 * Kcx * px * (1.0 + px); |
| 699 |
|
a001 = Kcx * px * px; |
| 700 |
< |
*/ |
| 672 |
< |
|
| 700 |
> |
*/ |
| 701 |
|
//scale all three dimensions, let c_x = c_y |
| 702 |
|
a000 = Kcx + Kcy; |
| 703 |
|
a110 = Kcz; |
| 707 |
|
b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py)); |
| 708 |
|
b11 = -2.0 * Khz * pz * (1.0 + pz); |
| 709 |
|
c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) |
| 710 |
< |
+ Khz * pz * (2.0 + pz) - targetFlux_; |
| 710 |
> |
+ Khz * pz * (2.0 + pz) - targetFlux_; |
| 711 |
|
break; |
| 712 |
|
case rnemdPxScale : |
| 713 |
|
c = 1 - targetFlux_ / Pcx; |
| 719 |
|
b01 = -2.0 * Khy * py * (1.0 + py); |
| 720 |
|
b11 = -2.0 * Khz * pz * (1.0 + pz); |
| 721 |
|
c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz) |
| 722 |
< |
+ Khx * (fastpow(c * px - px - 1.0, 2) - 1.0); |
| 722 |
> |
+ Khx * (fastpow(c * px - px - 1.0, 2) - 1.0); |
| 723 |
|
break; |
| 724 |
|
case rnemdPyScale : |
| 725 |
|
c = 1 - targetFlux_ / Pcy; |
| 731 |
|
b01 = -2.0 * Khx * px * (1.0 + px); |
| 732 |
|
b11 = -2.0 * Khz * pz * (1.0 + pz); |
| 733 |
|
c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz) |
| 734 |
< |
+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); |
| 734 |
> |
+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); |
| 735 |
|
break; |
| 736 |
|
case rnemdPzScale ://we don't really do this, do we? |
| 737 |
|
c = 1 - targetFlux_ / Pcz; |
| 779 |
|
poly.setCoefficient(2, u2); |
| 780 |
|
poly.setCoefficient(1, u1); |
| 781 |
|
poly.setCoefficient(0, u0); |
| 782 |
< |
std::vector<RealType> realRoots = poly.FindRealRoots(); |
| 782 |
> |
vector<RealType> realRoots = poly.FindRealRoots(); |
| 783 |
|
|
| 784 |
< |
std::vector<RealType>::iterator ri; |
| 784 |
> |
vector<RealType>::iterator ri; |
| 785 |
|
RealType r1, r2, alpha0; |
| 786 |
< |
std::vector<std::pair<RealType,RealType> > rps; |
| 786 |
> |
vector<pair<RealType,RealType> > rps; |
| 787 |
|
for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) { |
| 788 |
|
r2 = *ri; |
| 789 |
|
//check if FindRealRoots() give the right answer |
| 799 |
|
if (alpha0 >= 0.0) { |
| 800 |
|
r1 = sqrt(alpha0 / a000); |
| 801 |
|
if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6) |
| 802 |
< |
{ rps.push_back(std::make_pair(r1, r2)); } |
| 802 |
> |
{ rps.push_back(make_pair(r1, r2)); } |
| 803 |
|
if (r1 > 1e-6) { //r1 non-negative |
| 804 |
|
r1 = -r1; |
| 805 |
|
if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6) |
| 806 |
< |
{ rps.push_back(std::make_pair(r1, r2)); } |
| 806 |
> |
{ rps.push_back(make_pair(r1, r2)); } |
| 807 |
|
} |
| 808 |
|
} |
| 809 |
|
} |
| 810 |
< |
// Consider combininig together the solving pair part w/ the searching |
| 810 |
> |
// Consider combining together the solving pair part w/ the searching |
| 811 |
|
// best solution part so that we don't need the pairs vector |
| 812 |
|
if (!rps.empty()) { |
| 813 |
|
RealType smallestDiff = HONKING_LARGE_VALUE; |
| 814 |
|
RealType diff; |
| 815 |
< |
std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0); |
| 816 |
< |
std::vector<std::pair<RealType,RealType> >::iterator rpi; |
| 815 |
> |
pair<RealType,RealType> bestPair = make_pair(1.0, 1.0); |
| 816 |
> |
vector<pair<RealType,RealType> >::iterator rpi; |
| 817 |
|
for (rpi = rps.begin(); rpi != rps.end(); rpi++) { |
| 818 |
|
r1 = (*rpi).first; |
| 819 |
|
r2 = (*rpi).second; |
| 845 |
|
#ifdef IS_MPI |
| 846 |
|
if (worldRank == 0) { |
| 847 |
|
#endif |
| 848 |
< |
std::cerr << "we choose r1 = " << bestPair.first |
| 849 |
< |
<< " and r2 = " << bestPair.second << "\n"; |
| 848 |
> |
sprintf(painCave.errMsg, |
| 849 |
> |
"RNEMD: roots r1= %lf\tr2 = %lf\n", |
| 850 |
> |
bestPair.first, bestPair.second); |
| 851 |
> |
painCave.isFatal = 0; |
| 852 |
> |
painCave.severity = OPENMD_INFO; |
| 853 |
> |
simError(); |
| 854 |
|
#ifdef IS_MPI |
| 855 |
|
} |
| 856 |
|
#endif |
| 857 |
< |
|
| 857 |
> |
|
| 858 |
|
RealType x, y, z; |
| 859 |
< |
switch(rnemdType_) { |
| 860 |
< |
case rnemdKineticScale : |
| 861 |
< |
x = bestPair.first; |
| 862 |
< |
y = bestPair.first; |
| 863 |
< |
z = bestPair.second; |
| 864 |
< |
break; |
| 865 |
< |
case rnemdPxScale : |
| 866 |
< |
x = c; |
| 867 |
< |
y = bestPair.first; |
| 868 |
< |
z = bestPair.second; |
| 869 |
< |
break; |
| 870 |
< |
case rnemdPyScale : |
| 871 |
< |
x = bestPair.first; |
| 872 |
< |
y = c; |
| 873 |
< |
z = bestPair.second; |
| 874 |
< |
break; |
| 875 |
< |
case rnemdPzScale : |
| 876 |
< |
x = bestPair.first; |
| 877 |
< |
y = bestPair.second; |
| 878 |
< |
z = c; |
| 879 |
< |
break; |
| 880 |
< |
default : |
| 881 |
< |
break; |
| 882 |
< |
} |
| 883 |
< |
std::vector<StuntDouble*>::iterator sdi; |
| 859 |
> |
switch(rnemdType_) { |
| 860 |
> |
case rnemdKineticScale : |
| 861 |
> |
x = bestPair.first; |
| 862 |
> |
y = bestPair.first; |
| 863 |
> |
z = bestPair.second; |
| 864 |
> |
break; |
| 865 |
> |
case rnemdPxScale : |
| 866 |
> |
x = c; |
| 867 |
> |
y = bestPair.first; |
| 868 |
> |
z = bestPair.second; |
| 869 |
> |
break; |
| 870 |
> |
case rnemdPyScale : |
| 871 |
> |
x = bestPair.first; |
| 872 |
> |
y = c; |
| 873 |
> |
z = bestPair.second; |
| 874 |
> |
break; |
| 875 |
> |
case rnemdPzScale : |
| 876 |
> |
x = bestPair.first; |
| 877 |
> |
y = bestPair.second; |
| 878 |
> |
z = c; |
| 879 |
> |
break; |
| 880 |
> |
default : |
| 881 |
> |
break; |
| 882 |
> |
} |
| 883 |
> |
vector<StuntDouble*>::iterator sdi; |
| 884 |
|
Vector3d vel; |
| 885 |
|
for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { |
| 886 |
|
vel = (*sdi)->getVel(); |
| 903 |
|
exchangeSum_ += targetFlux_; |
| 904 |
|
//we may want to check whether the exchange has been successful |
| 905 |
|
} else { |
| 906 |
< |
std::cerr << "exchange NOT performed!\n";//MPI incompatible |
| 906 |
> |
sprintf(painCave.errMsg, |
| 907 |
> |
"RNEMD: exchange NOT performed!\n"); |
| 908 |
> |
painCave.isFatal = 0; |
| 909 |
> |
painCave.severity = OPENMD_INFO; |
| 910 |
> |
simError(); |
| 911 |
|
failTrialCount_++; |
| 912 |
|
} |
| 913 |
|
|
| 945 |
|
StuntDouble* sd; |
| 946 |
|
int idx; |
| 947 |
|
|
| 948 |
< |
// alternative approach, track all molecules instead of only those selected for scaling/swapping |
| 949 |
< |
//SimInfo::MoleculeIterator miter; |
| 950 |
< |
//std::vector<StuntDouble*>::iterator iiter; |
| 951 |
< |
//Molecule* mol; |
| 952 |
< |
//StuntDouble* integrableObject; |
| 953 |
< |
//for (mol = info_->beginMolecule(miter); mol != NULL; |
| 954 |
< |
// mol = info_->nextMolecule(miter)) |
| 955 |
< |
// integrableObject is essentially sd |
| 956 |
< |
//for (integrableObject = mol->beginIntegrableObject(iiter); |
| 957 |
< |
// integrableObject != NULL; |
| 958 |
< |
// integrableObject = mol->nextIntegrableObject(iiter)) |
| 948 |
> |
// alternative approach, track all molecules instead of only those |
| 949 |
> |
// selected for scaling/swapping: |
| 950 |
> |
/* |
| 951 |
> |
SimInfo::MoleculeIterator miter; |
| 952 |
> |
vector<StuntDouble*>::iterator iiter; |
| 953 |
> |
Molecule* mol; |
| 954 |
> |
StuntDouble* integrableObject; |
| 955 |
> |
for (mol = info_->beginMolecule(miter); mol != NULL; |
| 956 |
> |
mol = info_->nextMolecule(miter)) |
| 957 |
> |
integrableObject is essentially sd |
| 958 |
> |
for (integrableObject = mol->beginIntegrableObject(iiter); |
| 959 |
> |
integrableObject != NULL; |
| 960 |
> |
integrableObject = mol->nextIntegrableObject(iiter)) |
| 961 |
> |
*/ |
| 962 |
|
for (sd = seleMan_.beginSelected(selei); sd != NULL; |
| 963 |
|
sd = seleMan_.nextSelected(selei)) { |
| 964 |
|
|
| 976 |
|
|
| 977 |
|
int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) % |
| 978 |
|
rnemdLogWidth_; |
| 979 |
< |
/* no symmetrization allowed due to arbitary rnemdLogWidth_ value |
| 979 |
> |
// no symmetrization allowed due to arbitary rnemdLogWidth_ value |
| 980 |
> |
/* |
| 981 |
|
if (rnemdLogWidth_ == midBin_ + 1) |
| 982 |
|
if (binNo > midBin_) |
| 983 |
|
binNo = nBins_ - binNo; |
| 991 |
|
case rnemdKineticSwap : |
| 992 |
|
case rnemdKineticScale : |
| 993 |
|
|
| 994 |
< |
value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + |
| 955 |
< |
vel[2]*vel[2]); |
| 994 |
> |
value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); |
| 995 |
|
|
| 996 |
|
valueCount_[binNo] += 3; |
| 997 |
|
if (sd->isDirectional()) { |
| 1004 |
|
int k = (i + 2) % 3; |
| 1005 |
|
value += angMom[j] * angMom[j] / I(j, j) + |
| 1006 |
|
angMom[k] * angMom[k] / I(k, k); |
| 1007 |
< |
|
| 1007 |
> |
|
| 1008 |
|
valueCount_[binNo] +=2; |
| 1009 |
< |
|
| 1009 |
> |
|
| 1010 |
|
} else { |
| 1011 |
|
value += angMom[0]*angMom[0]/I(0, 0) |
| 1012 |
|
+ angMom[1]*angMom[1]/I(1, 1) |
| 1015 |
|
} |
| 1016 |
|
} |
| 1017 |
|
value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb; |
| 1018 |
< |
|
| 1018 |
> |
|
| 1019 |
|
break; |
| 1020 |
|
case rnemdPx : |
| 1021 |
|
case rnemdPxScale : |
| 1058 |
|
void RNEMD::getStarted() { |
| 1059 |
|
collectData(); |
| 1060 |
|
/* now should be able to output profile in step 0, but might not be useful |
| 1061 |
< |
Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 1062 |
< |
Stats& stat = currentSnap_->statData; |
| 1063 |
< |
stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; |
| 1061 |
> |
Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 1062 |
> |
Stats& stat = currentSnap_->statData; |
| 1063 |
> |
stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; |
| 1064 |
|
*/ |
| 1065 |
|
getStatus(); |
| 1066 |
|
} |