positionRescaleFactor2Tag,
// End multigrator
tcoupleCoefficientTag,
+ stochRescaleCoefficientTag,
minimizeCoefficientTag,
momentumCorrectionTag,
#if USE_BARRIER
SimpleBroadcastObject<Tensor> positionRescaleFactor2;
// End multigrator
SimpleBroadcastObject<BigReal> tcoupleCoefficient;
+ SimpleBroadcastObject<BigReal> stochRescaleCoefficient;
SimpleBroadcastObject<BigReal> minimizeCoefficient;
SimpleBroadcastObject<Vector> momentumCorrection;
#if USE_BARRIER
SimpleBroadcastObject<int> scriptBarrier;
SimpleBroadcastObject<int> traceBarrier;
SimpleBroadcastObject<Vector> accelMDRescaleFactor;
- SimpleBroadcastObject<BigReal> adaptTemperature;
+ SimpleBroadcastObject<BigReal> adaptTemperature;
#ifdef MEASURE_NAMD_WITH_PAPI
SimpleBroadcastObject<int> papiMeasureBarrier;
#endif
- ControllerBroadcasts(const LDObjHandle *ldObjPtr = 0) :
+ ControllerBroadcasts(const LDObjHandle *ldObjPtr = 0) :
velocityRescaleFactor(velocityRescaleFactorTag, ldObjPtr),
positionRescaleFactor(positionRescaleFactorTag, ldObjPtr),
// For multigrator
positionRescaleFactor2(positionRescaleFactor2Tag, ldObjPtr),
// End multigrator
tcoupleCoefficient(tcoupleCoefficientTag, ldObjPtr),
+ stochRescaleCoefficient(stochRescaleCoefficientTag, ldObjPtr),
minimizeCoefficient(minimizeCoefficientTag, ldObjPtr),
momentumCorrection(momentumCorrectionTag, ldObjPtr),
#if USE_BARRIER
cycleBarrier(cycleBarrierTag, ldObjPtr),
#endif
accelMDRescaleFactor(accelMDRescaleFactorTag, ldObjPtr),
- adaptTemperature(adaptTemperatureTag, ldObjPtr),
+ adaptTemperature(adaptTemperatureTag, ldObjPtr),
scriptBarrier(scriptBarrierTag, ldObjPtr),
#ifdef MEASURE_NAMD_WITH_PAPI
papiMeasureBarrier(papiMeasureTag, ldObjPtr),
random = new Random(simParams->randomSeed);
random->split(0,PatchMap::Object()->numPatches()+1);
- rescaleVelocities_sumTemps = 0; rescaleVelocities_numTemps = 0;
+ rescaleVelocities_sumTemps = 0;
+ rescaleVelocities_numTemps = 0;
+ if (simParams->stochRescaleOn) {
+ stochRescale_count = 0;
+ stochRescaleTimefactor =\
+ exp(-simParams->stochRescaleFreq*simParams->dt*0.001/simParams->stochRescalePeriod);
+ }
berendsenPressure_avg = 0; berendsenPressure_count = 0;
// strainRate tensor is symmetric to avoid rotation
langevinPiston_strainRate =
adaptTempUpdate(step);
rescaleVelocities(step);
tcoupleVelocities(step);
+ stochRescaleVelocities(step);
berendsenPressure(step);
langevinPiston1(step);
rescaleaccelMD(step);
}
}
+/*
+* Generate random numbers and broadcast the velocity rescaling factor.
+*/
+void Controller::stochRescaleVelocities(int step)
+{
+ if ( simParams->stochRescaleOn )
+ {
+ ++stochRescale_count;
+ if ( stochRescale_count == simParams->stochRescaleFreq )
+ {
+ const BigReal stochRescaleTemp = simParams->stochRescaleTemp;
+
+ BigReal coefficient = 1.;
+ if ( temperature > 0. )
+ {
+ BigReal R1 = random->gaussian();
+ BigReal gammaShape = 0.5*(numDegFreedom - 1);
+ // R2sum is the sum of (numDegFreedom - 1) squared normal variables, which is
+ // chi-squared distributed. This is in turn a special case of the Gamma
+ // distribution, which converges to a normal distribution in the limit of a
+ // large shape parameter.
+ BigReal R2sum = 2*(gammaShape + sqrt(gammaShape)*random->gaussian()) + R1*R1;
+ BigReal tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
+
+ coefficient = sqrt(stochRescaleTimefactor + (1 - stochRescaleTimefactor)*tempfactor*R2sum
+ + 2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*stochRescaleTimefactor));
+ }
+ broadcast->stochRescaleCoefficient.publish(step,coefficient);
+ stochRescale_count = 0;
+ }
+ }
+}
+
static char *FORMAT(BigReal X)
{
static char tmp_string[25];
int rescaleVelocities_numTemps;
void reassignVelocities(int);
void tcoupleVelocities(int);
+ void stochRescaleVelocities(int);
+ int stochRescale_count;
+ BigReal stochRescaleTimefactor;
void berendsenPressure(int);
// Tensor berendsenPressure_avg;
// int berendsenPressure_count;
random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
rescaleVelocities_numTemps = 0;
+ stochRescale_count = 0;
berendsenPressure_count = 0;
// patch->write_tip4_props();
}
rescaleVelocities(step);
tcoupleVelocities(timestep,step);
+ stochRescaleVelocities(timestep,step);
berendsenPressure(step);
if ( ! commOnly ) {
}
}
+/*
+* Rescale velocities with the scale factor sent from the corresponding Controller method.
+*/
+void Sequencer::stochRescaleVelocities(BigReal dt_fs, int step)
+{
+ if ( simParams->stochRescaleOn )
+ {
+ FullAtom *a = patch->atom.begin();
+ int numAtoms = patch->numAtoms;
+ ++stochRescale_count;
+ if ( stochRescale_count == simParams->stochRescaleFreq )
+ {
+ // Blocking receive for the temperature coupling coefficient.
+ BigReal coefficient = broadcast->stochRescaleCoefficient.get(step);
+ DebugM(4, "stochastically rescaling velocities at step " << step << " by " << coefficient << "\n");
+ for ( int i = 0; i < numAtoms; ++i )
+ {
+ a[i].velocity *= coefficient;
+ }
+ stochRescale_count = 0;
+ }
+ }
+}
+
void Sequencer::saveForce(const int ftag)
{
patch->saveForce(ftag);
void reinitVelocities(void);
void rescaleVelocitiesByFactor(BigReal);
void tcoupleVelocities(BigReal,int);
+ void stochRescaleVelocities(BigReal,int);
+ int stochRescale_count;
void berendsenPressure(int);
int berendsenPressure_count;
int checkpoint_berendsenPressure_count;
"containing the temperature coupling term B(i);\n"
"default is 'O'", PARSE_STRING);
+ opts.optionalB("main", "stochRescale",
+ "Should stochastic velocity rescaling be performed?",
+ &stochRescaleOn, FALSE);
+ opts.require("stochRescale", "stochRescaleTemp",
+ "Temperature for stochastic velocity rescaling",
+ &stochRescaleTemp);
+ opts.range("stochRescaleTemp", NOT_NEGATIVE);
+ opts.units("stochRescaleTemp", N_KELVIN);
+ opts.require("stochRescale", "stochRescalePeriod",
+ "Time scale for stochastic velocity rescaling (ps)",
+ &stochRescalePeriod);
+ opts.range("stochRescalePeriod", POSITIVE);
+ opts.require("stochRescale", "stochRescaleFreq",
+ "Number of steps between stochastic rescalings",
+ &stochRescaleFreq);
+ opts.range("stochRescaleFreq", POSITIVE);
+
opts.optional("main", "rescaleFreq", "Number of steps between "
"velocity rescaling", &rescaleFreq);
opts.range("rescaleFreq", POSITIVE);
// active
BigReal tCoupleTemp; // Temperature for temp coupling
+ Bool stochRescaleOn; // Flag TRUE-> stochastic velocity rescaling
+ BigReal stochRescaleTemp; // Temperature for stochastic velocity rescaling
+ BigReal stochRescalePeriod; // Timescale (ps) for stochastic velocity rescaling
+ int stochRescaleFreq; // How frequently stochastic velocity rescaling occurs
+
int rescaleFreq; // Velocity rescale frequency
BigReal rescaleTemp; // Temperature to rescale to