--- trunk/src/integrators/RNEMD.cpp 2009/10/19 13:39:04 1368 +++ branches/development/src/integrators/RNEMD.cpp 2011/09/13 22:05:04 1627 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,15 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ #include @@ -46,18 +46,19 @@ #include "math/Polynomial.hpp" #include "primitives/Molecule.hpp" #include "primitives/StuntDouble.hpp" -#include "utils/OOPSEConstant.hpp" +#include "utils/PhysicalConstants.hpp" #include "utils/Tuple.hpp" #ifndef IS_MPI #include "math/SeqRandNumGen.hpp" #else +#include #include "math/ParallelRandNumGen.hpp" #endif #define HONKING_LARGE_VALUE 1.0e10 -namespace oopse { +namespace OpenMD { RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { @@ -154,7 +155,7 @@ namespace oopse { midBin_ = nBins_ / 2; if (simParams->haveRNEMD_logWidth()) { rnemdLogWidth_ = simParams->getRNEMD_logWidth(); - if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) { + if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) { std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; std::cerr << "Automaically set back to default.\n"; rnemdLogWidth_ = nBins_; @@ -194,11 +195,11 @@ namespace oopse { RNEMD::~RNEMD() { delete randNumGen_; - - std::cerr << "total fail trials: " << failTrialCount_ << "\n"; + #ifdef IS_MPI if (worldRank == 0) { #endif + std::cerr << "total fail trials: " << failTrialCount_ << "\n"; rnemdLog_.close(); if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n"; @@ -279,7 +280,7 @@ namespace oopse { } //make exchangeSum_ comparable between swap & scale //temporarily without using energyConvert - //value = value * 0.5 / OOPSEConstant::energyConvert; + //value = value * 0.5 / PhysicalConstants::energyConvert; value *= 0.5; break; case rnemdPx : @@ -642,7 +643,7 @@ namespace oopse { a001 = Kcx * px * px; */ - //scale all three dimensions, let x = y + //scale all three dimensions, let c_x = c_y a000 = Kcx + Kcy; a110 = Kcz; c0 = targetFlux_ - Kcx - Kcy - Kcz; @@ -678,6 +679,17 @@ namespace oopse { + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); break; case rnemdPzScale ://we don't really do this, do we? + c = 1 - targetFlux_ / Pcz; + a000 = Kcx; + a110 = Kcy; + c0 = Kcz * c * c - Kcx - Kcy - Kcz; + a001 = px * px * Khx; + a111 = py * py * Khy; + b01 = -2.0 * Khx * px * (1.0 + px); + b11 = -2.0 * Khy * py * (1.0 + py); + c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) + + Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0); + break; default : break; } @@ -721,7 +733,10 @@ namespace oopse { r2 = *ri; //check if FindRealRoots() give the right answer if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) { - std::cerr << "WARNING! eq solvers might have mistakes!\n"; + sprintf(painCave.errMsg, + "RNEMD Warning: polynomial solve seems to have an error!"); + painCave.isFatal = 0; + simError(); failRootCount_++; } //might not be useful w/o rescaling coefficients @@ -797,6 +812,10 @@ namespace oopse { z = bestPair.second; break; case rnemdPzScale : + x = bestPair.first; + y = bestPair.second; + z = c; + break; default : break; } @@ -823,7 +842,7 @@ namespace oopse { exchangeSum_ += targetFlux_; //we may want to check whether the exchange has been successful } else { - std::cerr << "exchange NOT performed!\n"; + std::cerr << "exchange NOT performed!\n";//MPI incompatible failTrialCount_++; } @@ -915,19 +934,19 @@ namespace oopse { valueCount_[binNo] +=3; } } - value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb; + value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb; break; case rnemdPx : case rnemdPxScale : value = mass * vel[0]; valueCount_[binNo]++; - xVal = mass * vel.x() * vel.x() / OOPSEConstant::energyConvert - / OOPSEConstant::kb; - yVal = mass * vel.y() * vel.y() / OOPSEConstant::energyConvert - / OOPSEConstant::kb; - zVal = mass * vel.z() * vel.z() / OOPSEConstant::energyConvert - / OOPSEConstant::kb; + xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert + / PhysicalConstants::kb; + yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert + / PhysicalConstants::kb; + zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert + / PhysicalConstants::kb; xTempHist_[binNo] += xVal; yTempHist_[binNo] += yVal; zTempHist_[binNo] += zVal; @@ -965,8 +984,8 @@ namespace oopse { stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; //or to be more meaningful, define another item as exchangeSum_ / time + int j; - #ifdef IS_MPI // all processors have the same number of bins, and STL vectors pack their @@ -976,42 +995,52 @@ namespace oopse { rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0], rnemdLogWidth_, MPI::INT, MPI::SUM); - + if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) { + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0], + rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0], + rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0], + rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); + } // If we're the root node, should we print out the results int worldRank = MPI::COMM_WORLD.Get_rank(); if (worldRank == 0) { #endif - int j; rnemdLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j]; - valueHist_[j] = 0.0; } rnemdLog_ << "\n"; if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) { xTempLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j]; - xTempHist_[j] = 0.0; } xTempLog_ << "\n"; yTempLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j]; - yTempHist_[j] = 0.0; } yTempLog_ << "\n"; zTempLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j]; - zTempHist_[j] = 0.0; } zTempLog_ << "\n"; } - for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0; #ifdef IS_MPI } #endif + for (j = 0; j < rnemdLogWidth_; j++) { + valueCount_[j] = 0; + valueHist_[j] = 0.0; + } + if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) + for (j = 0; j < rnemdLogWidth_; j++) { + xTempHist_[j] = 0.0; + yTempHist_[j] = 0.0; + zTempHist_[j] = 0.0; + } } - }