| 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 |  | } |