| 52 |
|
#ifndef IS_MPI |
| 53 |
|
#include "math/SeqRandNumGen.hpp" |
| 54 |
|
#else |
| 55 |
+ |
#include <mpi.h> |
| 56 |
|
#include "math/ParallelRandNumGen.hpp" |
| 57 |
|
#endif |
| 58 |
|
|
| 155 |
|
midBin_ = nBins_ / 2; |
| 156 |
|
if (simParams->haveRNEMD_logWidth()) { |
| 157 |
|
rnemdLogWidth_ = simParams->getRNEMD_logWidth(); |
| 158 |
< |
if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) { |
| 158 |
> |
if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) { |
| 159 |
|
std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; |
| 160 |
|
std::cerr << "Automaically set back to default.\n"; |
| 161 |
|
rnemdLogWidth_ = nBins_; |
| 195 |
|
|
| 196 |
|
RNEMD::~RNEMD() { |
| 197 |
|
delete randNumGen_; |
| 198 |
< |
|
| 198 |
< |
std::cerr << "total fail trials: " << failTrialCount_ << "\n"; |
| 198 |
> |
|
| 199 |
|
#ifdef IS_MPI |
| 200 |
|
if (worldRank == 0) { |
| 201 |
|
#endif |
| 202 |
+ |
std::cerr << "total fail trials: " << failTrialCount_ << "\n"; |
| 203 |
|
rnemdLog_.close(); |
| 204 |
|
if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) |
| 205 |
|
std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n"; |
| 643 |
|
a001 = Kcx * px * px; |
| 644 |
|
*/ |
| 645 |
|
|
| 646 |
< |
//scale all three dimensions, let x = y |
| 646 |
> |
//scale all three dimensions, let c_x = c_y |
| 647 |
|
a000 = Kcx + Kcy; |
| 648 |
|
a110 = Kcz; |
| 649 |
|
c0 = targetFlux_ - Kcx - Kcy - Kcz; |
| 679 |
|
+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); |
| 680 |
|
break; |
| 681 |
|
case rnemdPzScale ://we don't really do this, do we? |
| 682 |
+ |
c = 1 - targetFlux_ / Pcz; |
| 683 |
+ |
a000 = Kcx; |
| 684 |
+ |
a110 = Kcy; |
| 685 |
+ |
c0 = Kcz * c * c - Kcx - Kcy - Kcz; |
| 686 |
+ |
a001 = px * px * Khx; |
| 687 |
+ |
a111 = py * py * Khy; |
| 688 |
+ |
b01 = -2.0 * Khx * px * (1.0 + px); |
| 689 |
+ |
b11 = -2.0 * Khy * py * (1.0 + py); |
| 690 |
+ |
c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) |
| 691 |
+ |
+ Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0); |
| 692 |
+ |
break; |
| 693 |
|
default : |
| 694 |
|
break; |
| 695 |
|
} |
| 733 |
|
r2 = *ri; |
| 734 |
|
//check if FindRealRoots() give the right answer |
| 735 |
|
if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) { |
| 736 |
< |
std::cerr << "WARNING! eq solvers might have mistakes!\n"; |
| 736 |
> |
sprintf(painCave.errMsg, |
| 737 |
> |
"RNEMD Warning: polynomial solve seems to have an error!"); |
| 738 |
> |
painCave.isFatal = 0; |
| 739 |
> |
simError(); |
| 740 |
|
failRootCount_++; |
| 741 |
|
} |
| 742 |
|
//might not be useful w/o rescaling coefficients |
| 812 |
|
z = bestPair.second; |
| 813 |
|
break; |
| 814 |
|
case rnemdPzScale : |
| 815 |
+ |
x = bestPair.first; |
| 816 |
+ |
y = bestPair.second; |
| 817 |
+ |
z = c; |
| 818 |
+ |
break; |
| 819 |
|
default : |
| 820 |
|
break; |
| 821 |
|
} |
| 842 |
|
exchangeSum_ += targetFlux_; |
| 843 |
|
//we may want to check whether the exchange has been successful |
| 844 |
|
} else { |
| 845 |
< |
std::cerr << "exchange NOT performed!\n"; |
| 845 |
> |
std::cerr << "exchange NOT performed!\n";//MPI incompatible |
| 846 |
|
failTrialCount_++; |
| 847 |
|
} |
| 848 |
|
|
| 984 |
|
|
| 985 |
|
stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; |
| 986 |
|
//or to be more meaningful, define another item as exchangeSum_ / time |
| 987 |
+ |
int j; |
| 988 |
|
|
| 969 |
– |
|
| 989 |
|
#ifdef IS_MPI |
| 990 |
|
|
| 991 |
|
// all processors have the same number of bins, and STL vectors pack their |
| 995 |
|
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 996 |
|
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0], |
| 997 |
|
rnemdLogWidth_, MPI::INT, MPI::SUM); |
| 998 |
< |
|
| 998 |
> |
if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) { |
| 999 |
> |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0], |
| 1000 |
> |
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 1001 |
> |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0], |
| 1002 |
> |
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 1003 |
> |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0], |
| 1004 |
> |
rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); |
| 1005 |
> |
} |
| 1006 |
|
// If we're the root node, should we print out the results |
| 1007 |
|
int worldRank = MPI::COMM_WORLD.Get_rank(); |
| 1008 |
|
if (worldRank == 0) { |
| 1009 |
|
#endif |
| 984 |
– |
int j; |
| 1010 |
|
rnemdLog_ << time; |
| 1011 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1012 |
|
rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j]; |
| 988 |
– |
valueHist_[j] = 0.0; |
| 1013 |
|
} |
| 1014 |
|
rnemdLog_ << "\n"; |
| 1015 |
|
if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) { |
| 1016 |
|
xTempLog_ << time; |
| 1017 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1018 |
|
xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j]; |
| 995 |
– |
xTempHist_[j] = 0.0; |
| 1019 |
|
} |
| 1020 |
|
xTempLog_ << "\n"; |
| 1021 |
|
yTempLog_ << time; |
| 1022 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1023 |
|
yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j]; |
| 1001 |
– |
yTempHist_[j] = 0.0; |
| 1024 |
|
} |
| 1025 |
|
yTempLog_ << "\n"; |
| 1026 |
|
zTempLog_ << time; |
| 1027 |
|
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1028 |
|
zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j]; |
| 1007 |
– |
zTempHist_[j] = 0.0; |
| 1029 |
|
} |
| 1030 |
|
zTempLog_ << "\n"; |
| 1031 |
|
} |
| 1011 |
– |
for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0; |
| 1032 |
|
#ifdef IS_MPI |
| 1033 |
< |
} |
| 1033 |
> |
} |
| 1034 |
|
#endif |
| 1035 |
< |
|
| 1036 |
< |
|
| 1035 |
> |
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1036 |
> |
valueCount_[j] = 0; |
| 1037 |
> |
valueHist_[j] = 0.0; |
| 1038 |
> |
} |
| 1039 |
> |
if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) |
| 1040 |
> |
for (j = 0; j < rnemdLogWidth_; j++) { |
| 1041 |
> |
xTempHist_[j] = 0.0; |
| 1042 |
> |
yTempHist_[j] = 0.0; |
| 1043 |
> |
zTempHist_[j] = 0.0; |
| 1044 |
> |
} |
| 1045 |
|
} |
| 1018 |
– |
|
| 1046 |
|
} |
