#include <Controller.h>
|
|
Definition at line 123 of file Controller.C. References avg_count, AVGXY, berendsenPressure_avg, berendsenPressure_count, BigReal, broadcast, checkpoint_stored, drudeBondTemp, drudeBondTempAvg, drudeComTemp, drudeComTempAvg, groupPressure_avg, groupPressure_tavg, Tensor::identity(), langevinPiston_strainRate, ReductionMgr::Object(), ppbonded, ppint, ppnonbonded, pressure_avg, pressure_tavg, SimParameters::pressureProfileAtomTypes, pressureProfileAverage, pressureProfileCount, SimParameters::pressureProfileEwaldOn, SimParameters::pressureProfileFreq, SimParameters::pressureProfileOn, SimParameters::pressureProfileSlabs, pressureProfileSlabs, random, SimParameters::randomSeed, reduction, REDUCTIONS_BASIC, REDUCTIONS_PPROF_BONDED, REDUCTIONS_PPROF_INTERNAL, REDUCTIONS_PPROF_NONBONDED, rescaleVelocities_numTemps, rescaleVelocities_sumTemps, simParams, simParams, smooth2_avg, Random::split(), SimParameters::strainRate, SimParameters::strainRate2, Tensor::symmetric(), tavg_count, temp_avg, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, and ReductionMgr::willRequire(). 00123 : 00124 computeChecksum(0), marginViolations(0), pairlistWarnings(0), 00125 simParams(Node::Object()->simParameters), 00126 state(s), 00127 #ifdef MEM_OPT_VERSION 00128 collection(CollectionMasterHandler::Object()), 00129 #else 00130 collection(CollectionMaster::Object()), 00131 #endif 00132 startCTime(0), 00133 firstCTime(CmiTimer()), 00134 startWTime(0), 00135 firstWTime(CmiWallTimer()), 00136 startBenchTime(0), 00137 ldbSteps(0), 00138 fflush_count(3) 00139 { 00140 broadcast = new ControllerBroadcasts; 00141 reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC); 00142 // pressure profile reductions 00143 pressureProfileSlabs = 0; 00144 pressureProfileCount = 0; 00145 ppbonded = ppnonbonded = ppint = NULL; 00146 if (simParams->pressureProfileOn || simParams->pressureProfileEwaldOn) { 00147 int ntypes = simParams->pressureProfileAtomTypes; 00148 int nslabs = pressureProfileSlabs = simParams->pressureProfileSlabs; 00149 // number of partitions for pairwise interactions 00150 int npairs = (ntypes * (ntypes+1))/2; 00151 pressureProfileAverage = new BigReal[3*nslabs]; 00152 memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal)); 00153 int freq = simParams->pressureProfileFreq; 00154 if (simParams->pressureProfileOn) { 00155 ppbonded = new PressureProfileReduction(REDUCTIONS_PPROF_BONDED, 00156 nslabs, npairs, "BONDED", freq); 00157 ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 00158 nslabs, npairs, "NONBONDED", freq); 00159 // internal partitions by atom type, but not in a pairwise fashion 00160 ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL, 00161 nslabs, ntypes, "INTERNAL", freq); 00162 } else { 00163 // just doing Ewald, so only need nonbonded 00164 ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 00165 nslabs, npairs, "NONBONDED", freq); 00166 } 00167 } 00168 random = new Random(simParams->randomSeed); 00169 random->split(0,PatchMap::Object()->numPatches()+1); 00170 00171 rescaleVelocities_sumTemps = 0; rescaleVelocities_numTemps = 0; 00172 berendsenPressure_avg = 0; berendsenPressure_count = 0; 00173 // strainRate tensor is symmetric to avoid rotation 00174 langevinPiston_strainRate = 00175 Tensor::symmetric(simParams->strainRate,simParams->strainRate2); 00176 if ( ! simParams->useFlexibleCell ) { 00177 BigReal avg = trace(langevinPiston_strainRate) / 3.; 00178 langevinPiston_strainRate = Tensor::identity(avg); 00179 } else if ( simParams->useConstantRatio ) { 00180 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\ 00181 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz ); 00182 AVGXY(langevinPiston_strainRate); 00183 #undef AVGXY 00184 } 00185 smooth2_avg = XXXBIGREAL; 00186 temp_avg = 0; 00187 pressure_avg = 0; 00188 groupPressure_avg = 0; 00189 avg_count = 0; 00190 pressure_tavg = 0; 00191 groupPressure_tavg = 0; 00192 tavg_count = 0; 00193 checkpoint_stored = 0; 00194 drudeComTemp = 0; 00195 drudeBondTemp = 0; 00196 drudeComTempAvg = 0; 00197 drudeBondTempAvg = 0; 00198 }
|
|
|
Definition at line 200 of file Controller.C. 00201 {
00202 delete broadcast;
00203 delete reduction;
00204 delete ppbonded;
00205 delete ppnonbonded;
00206 delete ppint;
00207 delete [] pressureProfileAverage;
00208 delete random;
00209 }
|
|
|
Definition at line 229 of file Controller.C. References BackEnd::awaken(), berendsenPressure_avg, berendsenPressure_count, broadcast, checkpoint_berendsenPressure_avg, checkpoint_berendsenPressure_count, checkpoint_langevinPiston_strainRate, checkpoint_lattice, checkpoint_stored, END_OF_RUN, enqueueCollections(), EVAL_MEASURE, FILE_OUTPUT, SimParameters::firstTimestep, SimpleBroadcastObject< T >::get(), SimParameters::initialTemp, integrate(), iout, langevinPiston_strainRate, NamdState::lattice, minimize(), NAMD_die(), outputExtendedSystem(), SCRIPT_CHECKPOINT, SCRIPT_MEASURE, SCRIPT_MINIMIZE, SCRIPT_OUTPUT, SCRIPT_REINITVELS, SCRIPT_RESCALEVELS, SCRIPT_REVERT, SCRIPT_RUN, SimParameters::scriptArg1, ControllerBroadcasts::scriptBarrier, simParams, state, and terminate(). 00230 {
00231 int scriptTask;
00232 int scriptSeq = 0;
00233 BackEnd::awaken();
00234 while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
00235 switch ( scriptTask ) {
00236 case SCRIPT_OUTPUT:
00237 enqueueCollections(FILE_OUTPUT);
00238 outputExtendedSystem(FILE_OUTPUT);
00239 break;
00240 case SCRIPT_MEASURE:
00241 enqueueCollections(EVAL_MEASURE);
00242 break;
00243 case SCRIPT_REINITVELS:
00244 iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
00245 << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
00246 break;
00247 case SCRIPT_RESCALEVELS:
00248 iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
00249 << " BY " << simParams->scriptArg1 << "\n" << endi;
00250 break;
00251 case SCRIPT_CHECKPOINT:
00252 iout << "CHECKPOINTING POSITIONS AT STEP " << simParams->firstTimestep
00253 << "\n" << endi;
00254 checkpoint_stored = 1;
00255 checkpoint_lattice = state->lattice;
00256 checkpoint_langevinPiston_strainRate = langevinPiston_strainRate;
00257 checkpoint_berendsenPressure_avg = berendsenPressure_avg;
00258 checkpoint_berendsenPressure_count = berendsenPressure_count;
00259 break;
00260 case SCRIPT_REVERT:
00261 iout << "REVERTING POSITIONS AT STEP " << simParams->firstTimestep
00262 << "\n" << endi;
00263 if ( ! checkpoint_stored )
00264 NAMD_die("Unable to revert, checkpoint was never called!");
00265 state->lattice = checkpoint_lattice;
00266 langevinPiston_strainRate = checkpoint_langevinPiston_strainRate;
00267 berendsenPressure_avg = checkpoint_berendsenPressure_avg;
00268 berendsenPressure_count = checkpoint_berendsenPressure_count;
00269 break;
00270 case SCRIPT_MINIMIZE:
00271 minimize();
00272 break;
00273 case SCRIPT_RUN:
00274 integrate();
00275 break;
00276 }
00277 BackEnd::awaken();
00278 }
00279 enqueueCollections(END_OF_RUN);
00280 outputExtendedSystem(END_OF_RUN);
00281 // note: this is a Converse interface call that gets translated to a
00282 // null method call if CMK_OPTIMIZE is set.
00283 if(traceAvailable())
00284 traceClose();
00285 terminate();
00286 }
|
|
|
Definition at line 35 of file Controller.h. Referenced by LdbCoordinator::awakenSequencers(), and run(). 00035 { CthAwaken(thread); };
|
|
|
Definition at line 598 of file Controller.C. References berendsenPressure_avg, berendsenPressure_count, SimParameters::berendsenPressureCompressibility, SimParameters::berendsenPressureFreq, SimParameters::berendsenPressureOn, SimParameters::berendsenPressureRelaxationTime, SimParameters::berendsenPressureTarget, broadcast, Tensor::diagonal(), SimParameters::dt, Tensor::identity(), iERROR(), iout, NamdState::lattice, LIMIT_SCALING, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), Lattice::rescale(), simParams, state, Tensor::xx, Tensor::yy, and Tensor::zz. Referenced by integrate(). 00599 {
00600 if ( simParams->berendsenPressureOn ) {
00601 berendsenPressure_count += 1;
00602 berendsenPressure_avg += controlPressure;
00603 const int freq = simParams->berendsenPressureFreq;
00604 if ( ! (berendsenPressure_count % freq) ) {
00605 Tensor factor = berendsenPressure_avg / berendsenPressure_count;
00606 berendsenPressure_avg = 0;
00607 berendsenPressure_count = 0;
00608 // We only use on-diagonal terms (for now)
00609 factor = Tensor::diagonal(diagonal(factor));
00610 factor -= Tensor::identity(simParams->berendsenPressureTarget);
00611 factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
00612 simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
00613 factor += Tensor::identity(1.0);
00614 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
00615 if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
00616 if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
00617 int limited = 0;
00618 LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
00619 LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
00620 LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
00621 #undef LIMIT_SCALING
00622 if ( limited ) {
00623 iout << iERROR << "Step " << step <<
00624 " cell rescaling factor limited.\n" << endi;
00625 }
00626 broadcast->positionRescaleFactor.publish(step,factor);
00627 state->lattice.rescale(factor);
00628 }
00629 } else {
00630 berendsenPressure_avg = 0;
00631 berendsenPressure_count = 0;
00632 }
00633 }
|
|
||||||||||||
|
Definition at line 1192 of file Controller.C. References BigReal, computeChecksum, iERROR(), iINFO(), iout, RequireReduction::item(), iWARN(), marginViolations, Node::molecule, SimParameters::mollyOn, NAMD_bug(), Molecule::numAtoms, Molecule::numCalcAngles, Molecule::numCalcBonds, Molecule::numCalcCrossterms, Molecule::numCalcDihedrals, Molecule::numCalcExclusions, Molecule::numCalcImpropers, Node::Object(), SimParameters::outputPairlists, pairlistWarnings, reduction, REDUCTION_ANGLE_CHECKSUM, REDUCTION_ATOM_CHECKSUM, REDUCTION_BOND_CHECKSUM, REDUCTION_COMPUTE_CHECKSUM, REDUCTION_CROSSTERM_CHECKSUM, REDUCTION_DIHEDRAL_CHECKSUM, REDUCTION_EXCLUSION_CHECKSUM, REDUCTION_IMPROPER_CHECKSUM, REDUCTION_MARGIN_VIOLATIONS, REDUCTION_PAIRLIST_WARNINGS, REDUCTION_STRAY_CHARGE_ERRORS, and simParams. Referenced by printDynamicsEnergies(), and printMinimizeEnergies(). 01192 {
01193 Node *node = Node::Object();
01194 Molecule *molecule = node->molecule;
01195
01196 // Some basic correctness checking
01197 BigReal checksum, checksum_b;
01198
01199 checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
01200 if ( ((int)checksum) != molecule->numAtoms )
01201 NAMD_bug("Bad global atom count!\n");
01202
01203 checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
01204 if ( ((int)checksum) != computeChecksum ) {
01205 if ( computeChecksum )
01206 NAMD_bug("Bad global compute count!\n");
01207 else
01208 computeChecksum = ((int)checksum);
01209 }
01210
01211 checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
01212 if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
01213 if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
01214 iout << iWARN << "Bad global bond count!\n" << endi;
01215 else NAMD_bug("Bad global bond count!\n");
01216 }
01217
01218 checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
01219 if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
01220 if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
01221 iout << iWARN << "Bad global angle count!\n" << endi;
01222 else NAMD_bug("Bad global angle count!\n");
01223 }
01224
01225 checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
01226 if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
01227 if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
01228 iout << iWARN << "Bad global dihedral count!\n" << endi;
01229 else NAMD_bug("Bad global dihedral count!\n");
01230 }
01231
01232 checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
01233 if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
01234 if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
01235 iout << iWARN << "Bad global improper count!\n" << endi;
01236 else NAMD_bug("Bad global improper count!\n");
01237 }
01238
01239 checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
01240 if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
01241 if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
01242 iout << iWARN << "Bad global crossterm count!\n" << endi;
01243 else NAMD_bug("Bad global crossterm count!\n");
01244 }
01245
01246 checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01247 if ( ((int)checksum) > molecule->numCalcExclusions &&
01248 ( ! simParams->mollyOn || step % slowFreq ) )
01249 iout << iWARN << "Not all atoms have unique coordinates.\n" << endi;
01250
01251 checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
01252 if ( ((int)checksum) && ! marginViolations ) {
01253 iout << iERROR << "Margin is too small for " << ((int)checksum) <<
01254 " atoms during timestep " << step << ".\n" << iERROR <<
01255 "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
01256 }
01257 marginViolations += (int)checksum;
01258
01259 checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
01260 if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
01261 iout << iINFO <<
01262 "Pairlistdist is too small for " << ((int)checksum) <<
01263 " computes during timestep " << step << ".\n" << endi;
01264 }
01265 if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
01266
01267 checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
01268 if ( checksum ) {
01269 if ( forgiving )
01270 iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
01271 else NAMD_bug("Stray PME grid charges detected!\n");
01272 }
01273 }
|
|
|
Definition at line 856 of file Controller.C. References BigReal, broadcast, iERROR(), iout, RequireReduction::item(), Vector::length2(), ControllerBroadcasts::momentumCorrection, NAMD_die(), SimpleBroadcastObject< T >::publish(), reduction, REDUCTION_MOMENTUM_MASS, slowFreq, Vector::x, Vector::y, and Vector::z. Referenced by integrate(). 00856 {
00857
00858 Vector momentum;
00859 momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
00860 momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
00861 momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
00862 const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
00863
00864 if ( momentum.length2() == 0. )
00865 iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
00866 if ( mass == 0. )
00867 NAMD_die("Total mass is zero in Controller::correctMomentum");
00868
00869 momentum *= (-1./mass);
00870
00871 broadcast->momentumCorrection.publish(step+slowFreq,momentum);
00872 }
|
|
||||||||||||
|
Definition at line 2075 of file Controller.C. References broadcast. Referenced by integrate(). 02075 {
02076 #if USE_BARRIER
02077 if (doBarrier) {
02078 broadcast->cycleBarrier.publish(step,1);
02079 CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
02080 CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02081 }
02082 #endif
02083 }
|
|
|
Definition at line 1783 of file Controller.C. References collection, Output::coordinateNeeded(), CollectionMaster::enqueuePositions(), CollectionMaster::enqueueVelocities(), EnqueueDataMsg::l, NamdState::lattice, state, EnqueueDataMsg::timestep, and Output::velocityNeeded(). Referenced by algorithm(), integrate(), and minimize(). 01784 {
01785 if ( Output::coordinateNeeded(timestep) ){
01786 #ifdef MEM_OPT_VERSION
01787 EnqueueDataMsg *msg = new EnqueueDataMsg;
01788 msg->timestep = timestep;
01789 msg->l = state->lattice;
01790 collection->enqueuePositions(msg);
01791 #else
01792 collection->enqueuePositions(timestep,state->lattice);
01793 #endif
01794 }
01795 if ( Output::velocityNeeded(timestep) )
01796 collection->enqueueVelocities(timestep);
01797 }
|
|
|
Definition at line 304 of file Controller.C. References berendsenPressure(), correctMomentum(), cycleBarrier(), enqueueCollections(), eventEndOfTimeStep, SimParameters::firstTimestep, SimParameters::FMAOn, SimParameters::fullDirectOn, SimParameters::fullElectFrequency, langevinPiston1(), langevinPiston2(), SimParameters::N, nbondFreq, SimParameters::nonbondedFrequency, outputExtendedSystem(), outputFepEnergy(), outputTiEnergy(), SimParameters::PMEOn, printDynamicsEnergies(), printFepMessage(), printTiMessage(), reassignVelocities(), rebalanceLoad(), receivePressure(), rescaleVelocities(), simParams, slowFreq, SimParameters::stepsPerCycle, tcoupleVelocities(), and SimParameters::zeroMomentum. Referenced by algorithm(). 00304 {
00305
00306 int step = simParams->firstTimestep;
00307
00308 const int numberOfSteps = simParams->N;
00309 const int stepsPerCycle = simParams->stepsPerCycle;
00310
00311 const int zeroMomentum = simParams->zeroMomentum;
00312
00313 nbondFreq = simParams->nonbondedFrequency;
00314 const int dofull = ( simParams->fullDirectOn ||
00315 simParams->FMAOn || simParams->PMEOn );
00316 if (dofull)
00317 slowFreq = simParams->fullElectFrequency;
00318 else
00319 slowFreq = simParams->nonbondedFrequency;
00320
00321 reassignVelocities(step); // only for full-step velecities
00322 receivePressure(step);
00323 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00324 correctMomentum(step);
00325 printFepMessage(step);
00326 printTiMessage(step);
00327 printDynamicsEnergies(step);
00328 outputFepEnergy(step);
00329 outputTiEnergy(step);
00330 traceUserEvent(eventEndOfTimeStep);
00331 outputExtendedSystem(step);
00332 rebalanceLoad(step);
00333
00334 // Handling SIGINT doesn't seem to be working on Lemieux, and it
00335 // sometimes causes the net-xxx versions of NAMD to segfault on exit,
00336 // so disable it for now.
00337 // namd_sighandler_t oldhandler = signal(SIGINT,
00338 // (namd_sighandler_t)my_sigint_handler);
00339 for ( ++step ; step <= numberOfSteps; ++step )
00340 {
00341
00342 rescaleVelocities(step);
00343 tcoupleVelocities(step);
00344 berendsenPressure(step);
00345 langevinPiston1(step);
00346 enqueueCollections(step); // after lattice scaling!
00347 receivePressure(step);
00348 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00349 correctMomentum(step);
00350 langevinPiston2(step);
00351 reassignVelocities(step);
00352 printDynamicsEnergies(step);
00353 outputFepEnergy(step);
00354 outputTiEnergy(step);
00355 traceUserEvent(eventEndOfTimeStep);
00356 // if (gotsigint) {
00357 // iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
00358 // NAMD_quit();
00359 // }
00360 outputExtendedSystem(step);
00361 #if CYCLE_BARRIER
00362 cycleBarrier(!((step+1) % stepsPerCycle),step);
00363 #elif PME_BARRIER
00364 cycleBarrier(dofull && !(step%slowFreq),step);
00365 #elif STEP_BARRIER
00366 cycleBarrier(1, step);
00367 #endif
00368
00369 rebalanceLoad(step);
00370
00371 #if PME_BARRIER
00372 cycleBarrier(dofull && !((step+1)%slowFreq),step); // step before PME
00373 #endif
00374 }
00375 // signal(SIGINT, oldhandler);
00376 }
|
|
|
Definition at line 635 of file Controller.C. References BigReal, BOLTZMAN, broadcast, Lattice::c(), controlNumDegFreedom, controlPressure, Tensor::diagonal(), SimParameters::dt, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, Random::gaussian(), Random::gaussian_vector(), Tensor::identity(), iINFO(), iout, SimParameters::langevinPistonDecay, SimParameters::langevinPistonOn, SimParameters::langevinPistonPeriod, SimParameters::langevinPistonTarget, SimParameters::langevinPistonTemp, NamdState::lattice, nbondFreq, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), random, Lattice::rescale(), simParams, slowFreq, state, SimParameters::surfaceTensionTarget, SimParameters::useConstantArea, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Tensor::xx, Tensor::yy, Vector::z, and Tensor::zz. Referenced by integrate(). 00636 {
00637 if ( simParams->langevinPistonOn )
00638 {
00639 Tensor &strainRate = langevinPiston_strainRate;
00640 int cellDims = simParams->useFlexibleCell ? 1 : 3;
00641 BigReal dt = simParams->dt;
00642 BigReal dt_long = slowFreq * dt;
00643 BigReal kT = BOLTZMAN * simParams->langevinPistonTemp;
00644 BigReal tau = simParams->langevinPistonPeriod;
00645 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00646
00647 #ifdef DEBUG_PRESSURE
00648 iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
00649 #endif
00650
00651 if ( ! ( (step-1) % slowFreq ) )
00652 {
00653 BigReal gamma = 1 / simParams->langevinPistonDecay;
00654 BigReal f1 = exp( -0.5 * dt_long * gamma );
00655 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00656 strainRate *= f1;
00657 if ( simParams->useFlexibleCell ) {
00658 // We only use on-diagonal terms (for now)
00659 if ( simParams->useConstantRatio ) {
00660 BigReal r = f2 * random->gaussian();
00661 strainRate.xx += r;
00662 strainRate.yy += r;
00663 strainRate.zz += f2 * random->gaussian();
00664 } else {
00665 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00666 }
00667 } else {
00668 strainRate += f2 * Tensor::identity(random->gaussian());
00669 }
00670
00671 #ifdef DEBUG_PRESSURE
00672 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00673 #endif
00674 }
00675
00676 // Apply surface tension. If surfaceTensionTarget is zero, we get
00677 // the default (isotropic pressure) case.
00678
00679 Tensor ptarget;
00680 ptarget.zz = simParams->langevinPistonTarget;
00681 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
00682 simParams->surfaceTensionTarget / state->lattice.c().z;
00683
00684 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00685 ( controlPressure - ptarget );
00686
00687 if (simParams->fixCellDims) {
00688 if (simParams->fixCellDimX) strainRate.xx = 0;
00689 if (simParams->fixCellDimY) strainRate.yy = 0;
00690 if (simParams->fixCellDimZ) strainRate.zz = 0;
00691 }
00692
00693
00694 #ifdef DEBUG_PRESSURE
00695 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00696 #endif
00697
00698 if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
00699 {
00700 // We only use on-diagonal terms (for now)
00701 Tensor factor;
00702 if ( !simParams->useConstantArea ) {
00703 factor.xx = exp( dt_long * strainRate.xx );
00704 factor.yy = exp( dt_long * strainRate.yy );
00705 } else {
00706 factor.xx = factor.yy = 1;
00707 }
00708 factor.zz = exp( dt_long * strainRate.zz );
00709 broadcast->positionRescaleFactor.publish(step,factor);
00710 state->lattice.rescale(factor);
00711 #ifdef DEBUG_PRESSURE
00712 iout << iINFO << "rescaling by: " << factor << "\n";
00713 #endif
00714 }
00715
00716 // corrections to integrator
00717 if ( ! ( step % nbondFreq ) )
00718 {
00719 #ifdef DEBUG_PRESSURE
00720 iout << iINFO << "correcting strain rate for nbond, ";
00721 #endif
00722 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00723 ( (nbondFreq - 1) * controlPressure_nbond );
00724 #ifdef DEBUG_PRESSURE
00725 iout << "strain rate: " << strainRate << "\n";
00726 #endif
00727 }
00728 if ( ! ( step % slowFreq ) )
00729 {
00730 #ifdef DEBUG_PRESSURE
00731 iout << iINFO << "correcting strain rate for slow, ";
00732 #endif
00733 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00734 ( (slowFreq - 1) * controlPressure_slow );
00735 #ifdef DEBUG_PRESSURE
00736 iout << "strain rate: " << strainRate << "\n";
00737 #endif
00738 }
00739 if (simParams->fixCellDims) {
00740 if (simParams->fixCellDimX) strainRate.xx = 0;
00741 if (simParams->fixCellDimY) strainRate.yy = 0;
00742 if (simParams->fixCellDimZ) strainRate.zz = 0;
00743 }
00744
00745 }
00746
00747 }
|
|
|
Definition at line 749 of file Controller.C. References BigReal, BOLTZMAN, Lattice::c(), controlNumDegFreedom, controlPressure, Tensor::diagonal(), SimParameters::dt, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, Random::gaussian(), Random::gaussian_vector(), Tensor::identity(), iINFO(), iout, SimParameters::langevinPistonDecay, SimParameters::langevinPistonOn, SimParameters::langevinPistonPeriod, SimParameters::langevinPistonTarget, SimParameters::langevinPistonTemp, NamdState::lattice, nbondFreq, random, simParams, slowFreq, state, SimParameters::surfaceTensionTarget, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Tensor::xx, Tensor::yy, Vector::z, and Tensor::zz. Referenced by integrate(). 00750 {
00751 if ( simParams->langevinPistonOn )
00752 {
00753 Tensor &strainRate = langevinPiston_strainRate;
00754 int cellDims = simParams->useFlexibleCell ? 1 : 3;
00755 BigReal dt = simParams->dt;
00756 BigReal dt_long = slowFreq * dt;
00757 BigReal kT = BOLTZMAN * simParams->langevinPistonTemp;
00758 BigReal tau = simParams->langevinPistonPeriod;
00759 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00760
00761 // corrections to integrator
00762 if ( ! ( step % nbondFreq ) )
00763 {
00764 #ifdef DEBUG_PRESSURE
00765 iout << iINFO << "correcting strain rate for nbond, ";
00766 #endif
00767 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00768 ( (nbondFreq - 1) * controlPressure_nbond );
00769 #ifdef DEBUG_PRESSURE
00770 iout << "strain rate: " << strainRate << "\n";
00771 #endif
00772 }
00773 if ( ! ( step % slowFreq ) )
00774 {
00775 #ifdef DEBUG_PRESSURE
00776 iout << iINFO << "correcting strain rate for slow, ";
00777 #endif
00778 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00779 ( (slowFreq - 1) * controlPressure_slow );
00780 #ifdef DEBUG_PRESSURE
00781 iout << "strain rate: " << strainRate << "\n";
00782 #endif
00783 }
00784
00785 // Apply surface tension. If surfaceTensionTarget is zero, we get
00786 // the default (isotropic pressure) case.
00787
00788 Tensor ptarget;
00789 ptarget.zz = simParams->langevinPistonTarget;
00790 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
00791 simParams->surfaceTensionTarget / state->lattice.c().z;
00792
00793 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00794 ( controlPressure - ptarget );
00795
00796
00797 #ifdef DEBUG_PRESSURE
00798 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00799 #endif
00800
00801 if ( ! ( step % slowFreq ) )
00802 {
00803 BigReal gamma = 1 / simParams->langevinPistonDecay;
00804 BigReal f1 = exp( -0.5 * dt_long * gamma );
00805 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00806 strainRate *= f1;
00807 if ( simParams->useFlexibleCell ) {
00808 // We only use on-diagonal terms (for now)
00809 if ( simParams->useConstantRatio ) {
00810 BigReal r = f2 * random->gaussian();
00811 strainRate.xx += r;
00812 strainRate.yy += r;
00813 strainRate.zz += f2 * random->gaussian();
00814 } else {
00815 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00816 }
00817 } else {
00818 strainRate += f2 * Tensor::identity(random->gaussian());
00819 }
00820 #ifdef DEBUG_PRESSURE
00821 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00822 #endif
00823 }
00824
00825 #ifdef DEBUG_PRESSURE
00826 iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
00827 #endif
00828 if (simParams->fixCellDims) {
00829 if (simParams->fixCellDimX) strainRate.xx = 0;
00830 if (simParams->fixCellDimY) strainRate.yy = 0;
00831 if (simParams->fixCellDimZ) strainRate.zz = 0;
00832 }
00833 }
00834
00835
00836 }
|
|
|
Definition at line 419 of file Controller.C. References BigReal, broadcast, CALCULATE, minpoint::dudx, enqueueCollections(), SimParameters::firstTimestep, iout, min_f_dot_f, min_f_dot_v, min_huge_count, min_v_dot_v, SimParameters::minBabyStep, ControllerBroadcasts::minimizeCoefficient, SimParameters::minLineGoal, SimParameters::minTinyStep, SimParameters::minVerbose, MOVETO, SimParameters::N, nbondFreq, minpoint::noGradient, PRINT_BRACKET, SimpleBroadcastObject< T >::publish(), simParams, slowFreq, minpoint::u, and minpoint::x. Referenced by algorithm(). 00419 {
00420 // iout << "Controller::minimize() called.\n" << endi;
00421
00422 const int minVerbose = simParams->minVerbose;
00423 const int numberOfSteps = simParams->N;
00424 int step = simParams->firstTimestep;
00425 slowFreq = nbondFreq = 1;
00426 BigReal tinystep = simParams->minTinyStep; // 1.0e-6
00427 BigReal babystep = simParams->minBabyStep; // 1.0e-2
00428 BigReal linegoal = simParams->minLineGoal; // 1.0e-3
00429 BigReal initstep = tinystep;
00430 const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
00431
00432 CALCULATE
00433
00434 int minSeq = 0;
00435
00436 // just move downhill until initial bad contacts go away or 100 steps
00437 int old_min_huge_count = -1;
00438 int steps_at_same_min_huge_count = 0;
00439 for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
00440 if ( min_huge_count == old_min_huge_count ) {
00441 if ( ++steps_at_same_min_huge_count > 10 ) break;
00442 } else {
00443 old_min_huge_count = min_huge_count;
00444 steps_at_same_min_huge_count = 0;
00445 }
00446 iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
00447 broadcast->minimizeCoefficient.publish(minSeq++,1.);
00448 enqueueCollections(step);
00449 CALCULATE
00450 }
00451 if ( min_huge_count ) {
00452 iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
00453 }
00454 iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00455
00456 int atStart = 2;
00457 int errorFactor = 10;
00458 BigReal old_f_dot_f = min_f_dot_f;
00459 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00460 broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
00461 int newDir = 1;
00462 min_f_dot_v = min_f_dot_f;
00463 min_v_dot_v = min_f_dot_f;
00464 while ( 1 ) {
00465 // line minimization
00466 // bracket minimum on line
00467 minpoint lo,hi,mid,last;
00468 BigReal x = 0;
00469 lo.x = x;
00470 lo.u = min_energy;
00471 lo.dudx = -1. * min_f_dot_v;
00472 lo.noGradient = min_huge_count;
00473 mid = lo;
00474 last = mid;
00475 if ( minVerbose ) {
00476 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
00477 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
00478 iout << "\n" << endi;
00479 }
00480 BigReal tol = fabs( linegoal * min_f_dot_v );
00481 if ( initstep > babystep ) initstep = babystep;
00482 if ( initstep < 1.0e-300 ) initstep = 1.0e-300;
00483 iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
00484 fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
00485 x = initstep;
00486 x *= sqrt( min_f_dot_f / min_v_dot_v ); MOVETO(x)
00487 while ( min_huge_count ) {
00488 x *= 0.25; MOVETO(x);
00489 initstep *= 0.25;
00490 }
00491 // bracket minimum on line
00492 initstep *= 0.25;
00493 BigReal maxinitstep = initstep * 16.0;
00494 while ( last.u < mid.u ) {
00495 initstep *= 2.0;
00496 lo = mid; mid = last;
00497 x *= 2.0; MOVETO(x)
00498 }
00499 if ( initstep > maxinitstep ) initstep = maxinitstep;
00500 hi = last;
00501 #define PRINT_BRACKET \
00502 iout << "LINE MINIMIZER BRACKET: DX " \
00503 << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
00504 " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
00505 lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
00506 PRINT_BRACKET
00507 // converge on minimum on line
00508 int itcnt;
00509 for ( itcnt = 10; itcnt > 0; --itcnt ) {
00510 int progress = 1;
00511 // select new position
00512 if ( mid.noGradient ) {
00513 if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) { // subdivide left side
00514 x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
00515 MOVETO(x)
00516 if ( last.u <= mid.u ) {
00517 hi = mid; mid = last;
00518 } else {
00519 lo = last;
00520 }
00521 } else { // subdivide right side
00522 x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
00523 MOVETO(x)
00524 if ( last.u <= mid.u ) {
00525 lo = mid; mid = last;
00526 } else {
00527 hi = last;
00528 }
00529 }
00530 } else {
00531 if ( mid.dudx > 0. ) { // subdivide left side
00532 BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
00533 BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
00534 x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
00535 BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
00536 if ( xdiv ) x /= xdiv; else x = last.x;
00537 if ( x > altxhi ) x = altxhi;
00538 if ( x < altxlo ) x = altxlo;
00539 if ( x-last.x == 0 ) break;
00540 MOVETO(x)
00541 if ( last.u <= mid.u ) {
00542 hi = mid; mid = last;
00543 } else {
00544 if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
00545 lo = last;
00546 }
00547 } else { // subdivide right side
00548 BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
00549 BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
00550 x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
00551 BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
00552 if ( xdiv ) x /= xdiv; else x = last.x;
00553 if ( x < altxlo ) x = altxlo;
00554 if ( x > altxhi ) x = altxhi;
00555 if ( x-last.x == 0 ) break;
00556 MOVETO(x)
00557 if ( last.u <= mid.u ) {
00558 lo = mid; mid = last;
00559 } else {
00560 if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
00561 hi = last;
00562 }
00563 }
00564 }
00565 PRINT_BRACKET
00566 if ( ! progress ) {
00567 MOVETO(mid.x);
00568 break;
00569 }
00570 if ( fabs(last.dudx) < tol ) break;
00571 if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
00572 if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
00573 }
00574 // new direction
00575 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00576 BigReal c = min_f_dot_f / old_f_dot_f;
00577 c = ( c > 1.5 ? 1.5 : c );
00578 if ( atStart ) { c = 0; --atStart; }
00579 if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
00580 c = 0;
00581 if ( errorFactor < 100 ) errorFactor += 10;
00582 }
00583 if ( c == 0 ) {
00584 iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00585 }
00586 broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
00587 newDir = 1;
00588 old_f_dot_f = min_f_dot_f;
00589 min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
00590 min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
00591 }
00592
00593 }
|
|
|
Definition at line 1971 of file Controller.C. References FILE_OUTPUT, SimParameters::firstTimestep, iout, SimParameters::N, NAMD_backup_file(), NAMD_err(), SimParameters::outputFilename, SimParameters::restartFilename, SimParameters::restartFrequency, SimParameters::restartSave, simParams, writeExtendedSystemData(), writeExtendedSystemLabels(), xstFile, SimParameters::xstFilename, and SimParameters::xstFrequency. Referenced by algorithm(), and integrate(). 01972 {
01973
01974 if ( step >= 0 ) {
01975
01976 // Write out eXtended System Trajectory (XST) file
01977 if ( simParams->xstFrequency &&
01978 ((step % simParams->xstFrequency) == 0) )
01979 {
01980 if ( ! xstFile.rdbuf()->is_open() )
01981 {
01982 NAMD_backup_file(simParams->xstFilename);
01983 xstFile.open(simParams->xstFilename);
01984 iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
01985 xstFile << "# NAMD extended system trajectory file" << std::endl;
01986 writeExtendedSystemLabels(xstFile);
01987 }
01988 writeExtendedSystemData(step,xstFile);
01989 xstFile.flush();
01990 }
01991
01992 // Write out eXtended System Configuration (XSC) files
01993 // Output a restart file
01994 if ( simParams->restartFrequency &&
01995 ((step % simParams->restartFrequency) == 0) &&
01996 (step != simParams->firstTimestep) )
01997 {
01998 iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
01999 << step << "\n" << endi;
02000 char fname[140];
02001 strcpy(fname, simParams->restartFilename);
02002 if ( simParams->restartSave ) {
02003 char timestepstr[20];
02004 sprintf(timestepstr,".%d",step);
02005 strcat(fname, timestepstr);
02006 }
02007 strcat(fname, ".xsc");
02008 NAMD_backup_file(fname,".old");
02009 std::ofstream xscFile(fname);
02010 if (!xscFile) {
02011 char err_msg[257];
02012 sprintf(err_msg, "Error opening XSC restart file %s",fname);
02013 NAMD_err(err_msg);
02014 }
02015 xscFile << "# NAMD extended system configuration restart file" << std::endl;
02016 writeExtendedSystemLabels(xscFile);
02017 writeExtendedSystemData(step,xscFile);
02018 if (!xscFile) {
02019 char err_msg[257];
02020 sprintf(err_msg, "Error writing XSC restart file %s",fname);
02021 NAMD_err(err_msg);
02022 }
02023 }
02024
02025 }
02026
02027 // Output final coordinates
02028 if (step == FILE_OUTPUT || step == END_OF_RUN)
02029 {
02030 iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
02031 << simParams->N << "\n" << endi;
02032 static char fname[140];
02033 strcpy(fname, simParams->outputFilename);
02034 strcat(fname, ".xsc");
02035 NAMD_backup_file(fname);
02036 std::ofstream xscFile(fname);
02037 if (!xscFile) {
02038 char err_msg[257];
02039 sprintf(err_msg, "Error opening XSC output file %s",fname);
02040 NAMD_err(err_msg);
02041 }
02042 xscFile << "# NAMD extended system configuration output file" << std::endl;
02043 writeExtendedSystemLabels(xscFile);
02044 writeExtendedSystemData(simParams->N,xscFile);
02045 if (!xscFile) {
02046 char err_msg[257];
02047 sprintf(err_msg, "Error writing XSC output file %s",fname);
02048 NAMD_err(err_msg);
02049 }
02050 }
02051
02052 // Close trajectory file
02053 if (step == END_OF_RUN) {
02054 if ( xstFile.rdbuf()->is_open() ) {
02055 xstFile.close();
02056 iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
02057 }
02058 }
02059
02060 }
|
|
|
Definition at line 1814 of file Controller.C. References SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchLambda, SimParameters::alchLambda2, SimParameters::alchOutFile, SimParameters::alchOutFreq, SimParameters::alchTemp, BigReal, BOLTZMAN, electEnergy, electEnergy_f, electEnergySlow, electEnergySlow_f, exp_dE_ByRT, fepFile, FepNo, fepSum, SimParameters::firstTimestep, iout, ljEnergy_f, SimParameters::N, NAMD_backup_file(), net_dE, simParams, and writeFepEnergyData(). Referenced by integrate(). 01814 {
01815 if (simParams->alchFepOn) {
01816 const int stepInRun = step - simParams->firstTimestep;
01817 const int alchEquilSteps = simParams->alchEquilSteps;
01818 const BigReal alchLambda = simParams->alchLambda;
01819 const BigReal alchLambda2 = simParams->alchLambda2;
01820 if (stepInRun == 0 || stepInRun == alchEquilSteps) {
01821 FepNo = 0;
01822 exp_dE_ByRT = 0.0;
01823 net_dE = 0.0;
01824 }
01825 BigReal dE = electEnergy_f + electEnergySlow_f + ljEnergy_f
01826 - (electEnergy + electEnergySlow + ljEnergy);
01827 BigReal RT = BOLTZMAN * simParams->alchTemp;
01828 FepNo++;
01829 exp_dE_ByRT += exp(-dE/RT);
01830 net_dE += dE;
01831
01832 if (stepInRun == 0) {
01833 if (!fepFile.rdbuf()->is_open()) {
01834 fepSum = 0.0;
01835 NAMD_backup_file(simParams->alchOutFile);
01836 fepFile.open(simParams->alchOutFile);
01837 iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
01838 fepFile << "# STEP Elec "
01839 << "vdW dE dE_avg Temp dG\n"
01840 << "# l l+dl "
01841 << " l l+dl E(l+dl)-E(l)" << std::endl;
01842 }
01843 fepFile << "#NEW FEP WINDOW: "
01844 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
01845 << alchLambda2 << std::endl;
01846 }
01847 if (stepInRun == alchEquilSteps) {
01848 fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
01849 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
01850 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
01851 }
01852 if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) {
01853 writeFepEnergyData(step, fepFile);
01854 fepFile.flush();
01855 }
01856 if (step == simParams->N) {
01857 fepSum = fepSum + dG;
01858 fepFile << "#Free energy change for lambda window [ " << alchLambda
01859 << " " << alchLambda2 << " ] is " << dG << " ; net change until now is " << fepSum << std::endl;
01860 }
01861 }
01862 }
|
|
|
Definition at line 1864 of file Controller.C. References SimParameters::alchElecLambdaStart, SimParameters::alchEquilSteps, SimParameters::alchLambda, SimParameters::alchOutFile, SimParameters::alchOutFreq, SimParameters::alchThermIntOn, SimParameters::alchVdwLambdaEnd, BigReal, electEnergy_ti_1, electEnergy_ti_2, electEnergySlow_ti_1, electEnergySlow_ti_2, SimParameters::firstTimestep, iout, NAMD_backup_file(), net_dEdl_elec_1, net_dEdl_elec_2, net_dEdl_lj_1, net_dEdl_lj_2, recent_dEdl_elec_1, recent_dEdl_elec_2, recent_dEdl_lj_1, recent_dEdl_lj_2, recent_TiNo, simParams, tiFile, TiNo, and writeTiEnergyData(). Referenced by integrate(). 01864 {
01865 if (simParams->alchThermIntOn) {
01866 const int stepInRun = step - simParams->firstTimestep;
01867 const int alchEquilSteps = simParams->alchEquilSteps;
01868 const BigReal alchLambda = simParams->alchLambda;
01869
01870 if (stepInRun == 0 || stepInRun == alchEquilSteps) {
01871 TiNo = 0;
01872 net_dEdl_elec_1 = 0;
01873 net_dEdl_elec_2 = 0;
01874 net_dEdl_lj_1 = 0;
01875 net_dEdl_lj_2 = 0;
01876 }
01877 if (stepInRun == 0 || (! ((step - 1) % simParams->alchOutFreq))) {
01878 // output of instantaneous dU/dl now replaced with running average
01879 // over last alchOutFreq steps (except for step 0)
01880 recent_TiNo = 0;
01881 recent_dEdl_elec_1 = 0;
01882 recent_dEdl_elec_2 = 0;
01883 recent_dEdl_lj_1 = 0;
01884 recent_dEdl_lj_2 = 0;
01885 }
01886 TiNo++;
01887 recent_TiNo++;
01888 // FB - PME is no longer scaled by global lambda, but by the respective
01889 // lambda as dictated by elecLambdaStart. All electrostatics now go together.
01890 net_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1;
01891 net_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2;
01892 net_dEdl_lj_1 += ljEnergy_ti_1;
01893 net_dEdl_lj_2 += ljEnergy_ti_2;
01894 recent_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1;
01895 recent_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2;
01896 recent_dEdl_lj_1 += ljEnergy_ti_1;
01897 recent_dEdl_lj_2 += ljEnergy_ti_2;
01898
01899 if (stepInRun == 0) {
01900 BigReal alchElecLambdaStart = simParams->alchElecLambdaStart;
01901 BigReal alchVdwLambdaEnd = simParams->alchVdwLambdaEnd;
01902 BigReal elec_lambda_1 = (alchLambda <= alchElecLambdaStart)? 0. : \
01903 (alchLambda - alchElecLambdaStart) / (1. - alchElecLambdaStart);
01904 BigReal elec_lambda_2 = ((1-alchLambda) <= alchElecLambdaStart)? 0. : \
01905 ((1-alchLambda) - alchElecLambdaStart) / (1. - alchElecLambdaStart);
01906 BigReal vdw_lambda_1 = (alchLambda >= alchVdwLambdaEnd)? 1. : \
01907 alchLambda / alchVdwLambdaEnd;
01908 BigReal vdw_lambda_2 = ((1-alchLambda) >= alchVdwLambdaEnd)? 1. : \
01909 (1-alchLambda) / alchVdwLambdaEnd;
01910 if (!tiFile.rdbuf()->is_open()) {
01911 //tiSum = 0.0;
01912 NAMD_backup_file(simParams->alchOutFile);
01913 tiFile.open(simParams->alchOutFile);
01914 iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
01915 tiFile << "# STEP Elec_dU/dl Elec_avg vdW_dU/dl vdw_avg Elec_dU/dl Elec_avg vdW_dU/dl vdw_avg PME_dU/dl PME_avg\n"
01916 << "# <---------------------PARTITION 1------------------------> <---------------------PARTITION 2--------------------->"
01917 << std::endl;
01918 }
01919 tiFile << "#NEW TI WINDOW: "
01920 << "LAMBDA " << alchLambda
01921 << "\n#PARTITION 1 VDW LAMBDA " << vdw_lambda_1
01922 << "\n#PARTITION 1 ELEC LAMBDA " << elec_lambda_1
01923 << "\n#PARTITION 2 VDW LAMBDA " << vdw_lambda_2
01924 << "\n#PARTITION 2 ELEC LAMBDA " << elec_lambda_2
01925 << "\n" << std::endl;
01926 }
01927 if (stepInRun == alchEquilSteps) {
01928 tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
01929 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
01930 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
01931 }
01932 if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) {
01933 writeTiEnergyData(step, tiFile);
01934 tiFile.flush();
01935 }
01936 }
01937 }
|
|
|
Definition at line 1330 of file Controller.C. References BigReal, compareChecksums(), RequireReduction::item(), NamdState::lattice, Node::molecule, NAMD_bug(), Molecule::numCalcExclusions, Node::Object(), printEnergies(), reduction, REDUCTION_EXCLUSION_CHECKSUM, Node::simParameters, and state. Referenced by integrate(). 01330 {
01331
01332 Node *node = Node::Object();
01333 Molecule *molecule = node->molecule;
01334 SimParameters *simParameters = node->simParameters;
01335 Lattice &lattice = state->lattice;
01336
01337 BigReal checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01338 if ( ((int)checksum) && ((int)checksum) < molecule->numCalcExclusions )
01339 NAMD_bug("Bad global exclusion count!\n");
01340 compareChecksums(step);
01341
01342 printEnergies(step,0);
01343 }
|
|
||||||||||||
|
Definition at line 1345 of file Controller.C. References Lattice::a(), Lattice::a_p(), avg_count, Lattice::b(), Lattice::b_p(), BigReal, Lattice::c(), Lattice::c_p(), CALLBACKDATA, CALLBACKLIST, ScriptTcl::doCallback(), drudeBondTemp, drudeBondTempAvg, drudeComTemp, drudeComTempAvg, SimParameters::drudeOn, SimParameters::dt, IMDEnergies::Eangle, IMDEnergies::Ebond, IMDEnergies::Edihe, IMDEnergies::Eelec, IMDEnergies::Eimpr, electEnergy, electEnergy_f, electEnergy_ti_1, electEnergy_ti_2, electEnergyPME_ti_1, electEnergyPME_ti_2, electEnergySlow, electEnergySlow_f, electEnergySlow_ti_1, electEnergySlow_ti_2, IMDEnergies::Epot, ETITLE(), IMDEnergies::Etot, IMDEnergies::Evdw, SimParameters::firstLdbStep, SimParameters::firstTimestep, FORMAT(), IMDOutput::gather_energies(), GET_VECTOR, PressureProfileReduction::getData(), Node::getScript(), groupPressure, groupPressure_avg, groupPressure_tavg, iERROR(), iINFO(), Node::imd, SimParameters::IMDfreq, SimParameters::IMDon, iout, RequireReduction::item(), iWARN(), kineticEnergyHalfstep, NamdState::lattice, SimParameters::LJcorrection, ljEnergy, ljEnergy_f, ljEnergy_ti_1, ljEnergy_ti_2, marginViolations, memusage_MB(), SimParameters::mergeCrossterms, Node::molecule, Node::Object(), Lattice::origin(), SimParameters::outputEnergies, SimParameters::outputMomenta, SimParameters::outputPairlists, SimParameters::outputPressure, SimParameters::pairInteractionOn, pairlistWarnings, ppbonded, ppint, ppnonbonded, pressure, pressure_avg, pressure_tavg, PRESSUREFACTOR, pressureProfileAverage, pressureProfileCount, SimParameters::pressureProfileFreq, printTiming(), reduction, REDUCTION_ANGLE_ENERGY, REDUCTION_BC_ENERGY, REDUCTION_BOND_ENERGY, REDUCTION_CROSSTERM_ENERGY, REDUCTION_DIHEDRAL_ENERGY, REDUCTION_ELECT_ENERGY, REDUCTION_ELECT_ENERGY_F, REDUCTION_ELECT_ENERGY_PME_TI_1, REDUCTION_ELECT_ENERGY_PME_TI_2, REDUCTION_ELECT_ENERGY_SLOW, REDUCTION_ELECT_ENERGY_SLOW_F, REDUCTION_ELECT_ENERGY_SLOW_TI_1, REDUCTION_ELECT_ENERGY_SLOW_TI_2, REDUCTION_ELECT_ENERGY_TI_1, REDUCTION_ELECT_ENERGY_TI_2, REDUCTION_IMPROPER_ENERGY, REDUCTION_LJ_ENERGY, REDUCTION_LJ_ENERGY_F, REDUCTION_LJ_ENERGY_TI_1, REDUCTION_LJ_ENERGY_TI_2, REDUCTION_MISC_ENERGY, REDUCTION_PAIR_ELECT_FORCE, REDUCTION_PAIR_VDW_FORCE, Node::simParameters, simParams, smooth2_avg, state, IMDEnergies::T, Molecule::tail_corr_ener, tavg_count, temp_avg, temperature, totalEnergy, IMDEnergies::tstep, Lattice::volume(), Vector::x, Vector::y, and Vector::z. Referenced by printDynamicsEnergies(), and printMinimizeEnergies(). 01346 {
01347 Node *node = Node::Object();
01348 Molecule *molecule = node->molecule;
01349 SimParameters *simParameters = node->simParameters;
01350 Lattice &lattice = state->lattice;
01351
01352 BigReal bondEnergy;
01353 BigReal angleEnergy;
01354 BigReal dihedralEnergy;
01355 BigReal improperEnergy;
01356 BigReal crosstermEnergy;
01357 BigReal boundaryEnergy;
01358 BigReal miscEnergy;
01359 BigReal potentialEnergy;
01360 BigReal flatEnergy;
01361 BigReal smoothEnergy;
01362 Vector momentum;
01363 Vector angularMomentum;
01364 BigReal volume = lattice.volume();
01365
01366 bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
01367 angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
01368 dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
01369 improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
01370 crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
01371 boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
01372 miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
01373
01374 if ( minimize || ! ( step % nbondFreq ) )
01375 {
01376 electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
01377 ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
01378 //fepb
01379 electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
01380 ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
01381
01382 electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
01383 ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
01384 electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
01385 ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
01386 //fepe
01387 }
01388
01389 if ( minimize || ! ( step % slowFreq ) )
01390 {
01391 electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
01392 //fepb
01393 electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
01394
01395 electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
01396 electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
01397 electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
01398 electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
01399 //fepe
01400 }
01401
01402 if (simParameters->LJcorrection && volume) {
01403 // Apply tail correction to energy
01404 //printf("Volume is %f\n", volume);
01405 //printf("Applying tail correction of %f to energy\n", molecule->tail_corr_ener / volume);
01406 ljEnergy += molecule->tail_corr_ener / volume;
01407 }
01408
01409
01410 momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
01411 momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
01412 momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
01413 angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
01414 angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
01415 angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
01416
01417 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
01418 improperEnergy + electEnergy + electEnergySlow + ljEnergy +
01419 crosstermEnergy + boundaryEnergy + miscEnergy;
01420 totalEnergy = potentialEnergy + kineticEnergy;
01421 flatEnergy = totalEnergy +
01422 (1.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
01423 if ( !(step%slowFreq) ) {
01424 // only adjust based on most accurate energies
01425 BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
01426 if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
01427 smooth2_avg *= 0.9375;
01428 smooth2_avg += 0.0625 * s;
01429 }
01430 smoothEnergy = flatEnergy + smooth2_avg -
01431 (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
01432
01433 if ( simParameters->outputMomenta && ! minimize &&
01434 ! ( step % simParameters->outputMomenta ) )
01435 {
01436 iout << "MOMENTUM: " << step
01437 << " P: " << momentum
01438 << " L: " << angularMomentum
01439 << "\n" << endi;
01440 }
01441
01442 if ( simParameters->outputPressure ) {
01443 pressure_tavg += pressure;
01444 groupPressure_tavg += groupPressure;
01445 tavg_count += 1;
01446 if ( minimize || ! ( step % simParameters->outputPressure ) ) {
01447 #ifdef NAMD_CUDA
01448 if (
01449 step < (simParams->firstTimestep + 10*simParams->outputPressure) ) {
01450 iout << iWARN << "Pressure tensor is incorrect. CUDA currently evaluates scalar pressure only.\n" << endi;
01451 }
01452 #endif
01453 iout << "PRESSURE: " << step << " "
01454 << PRESSUREFACTOR * pressure << "\n"
01455 << "GPRESSURE: " << step << " "
01456 << PRESSUREFACTOR * groupPressure << "\n";
01457 if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
01458 << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
01459 << "GPRESSAVG: " << step << " "
01460 << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n" << endi;
01461 pressure_tavg = 0;
01462 groupPressure_tavg = 0;
01463 tavg_count = 0;
01464 }
01465 }
01466
01467 // pressure profile reductions
01468 if (pressureProfileSlabs) {
01469 const int freq = simParams->pressureProfileFreq;
01470 const int arraysize = 3*pressureProfileSlabs;
01471
01472 BigReal *total = new BigReal[arraysize];
01473 memset(total, 0, arraysize*sizeof(BigReal));
01474 const int first = simParams->firstTimestep;
01475
01476 if (ppbonded) ppbonded->getData(first, step, lattice, total);
01477 if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
01478 if (ppint) ppint->getData(first, step, lattice, total);
01479 for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
01480 pressureProfileCount++;
01481
01482 if (!(step % freq)) {
01483 // convert NAMD internal virial to pressure in units of bar
01484 BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
01485
01486 iout << "PRESSUREPROFILE: " << step << " ";
01487 if (step == first) {
01488 // output pressure profile for this step
01489 for (int i=0; i<arraysize; i++) {
01490 iout << total[i] * scalefac << " ";
01491 }
01492 } else {
01493 // output pressure profile averaged over the last count steps.
01494 scalefac /= pressureProfileCount;
01495 for (int i=0; i<arraysize; i++)
01496 iout << pressureProfileAverage[i]*scalefac << " ";
01497 }
01498 iout << "\n" << endi;
01499
01500 // Clear the average for the next block
01501 memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
01502 pressureProfileCount = 0;
01503 }
01504 delete [] total;
01505 }
01506
01507 int stepInRun = step - simParams->firstTimestep;
01508 if ( stepInRun % simParams->firstLdbStep == 0 ) {
01509 int benchPhase = stepInRun / simParams->firstLdbStep;
01510 if ( benchPhase > 0 && benchPhase < 7 ) {
01511 #ifdef NAMD_CUDA
01512 if ( simParams->outputEnergies < 60 ) {
01513 iout << iWARN << "Energy evaluation is done on CPU, increase outputEnergies to improve performance.\n";
01514 }
01515 #endif
01516 iout << iINFO;
01517 if ( benchPhase < 4 ) iout << "Initial time: ";
01518 else iout << "Benchmark time: ";
01519 iout << CkNumPes() << " CPUs ";
01520 {
01521 BigReal wallTime = CmiWallTimer() - startBenchTime;
01522 BigReal wallPerStep =
01523 (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
01524 BigReal ns = simParams->dt / 1000000.0;
01525 BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
01526 BigReal daysPerNano = wallPerStep * days / ns;
01527 iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
01528 iout << memusage_MB() << " MB memory\n" << endi;
01529 }
01530 }
01531 startBenchTime = CmiWallTimer();
01532 }
01533
01534 printTiming(step);
01535
01536 Vector pairVDWForce, pairElectForce;
01537 if ( simParameters->pairInteractionOn ) {
01538 GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
01539 GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
01540 }
01541
01542 // callback to Tcl with whatever we can
01543 #ifdef NAMD_TCL
01544 #define CALLBACKDATA(LABEL,VALUE) \
01545 labels << (LABEL) << " "; values << (VALUE) << " ";
01546 #define CALLBACKLIST(LABEL,VALUE) \
01547 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
01548 if (node->getScript() && node->getScript()->doCallback()) {
01549 std::ostringstream labels, values;
01550 #if CMK_BLUEGENEL
01551 // the normal version below gives a compiler error
01552 values.precision(16);
01553 #else
01554 values << std::setprecision(16);
01555 #endif
01556 CALLBACKDATA("TS",step);
01557 CALLBACKDATA("BOND",bondEnergy);
01558 CALLBACKDATA("ANGLE",angleEnergy);
01559 CALLBACKDATA("DIHED",dihedralEnergy);
01560 CALLBACKDATA("CROSS",crosstermEnergy);
01561 CALLBACKDATA("IMPRP",improperEnergy);
01562 CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
01563 CALLBACKDATA("VDW",ljEnergy);
01564 CALLBACKDATA("BOUNDARY",boundaryEnergy);
01565 CALLBACKDATA("MISC",miscEnergy);
01566 CALLBACKDATA("POTENTIAL",potentialEnergy);
01567 CALLBACKDATA("KINETIC",kineticEnergy);
01568 CALLBACKDATA("TOTAL",totalEnergy);
01569 CALLBACKDATA("TEMP",temperature);
01570 CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
01571 CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
01572 CALLBACKDATA("VOLUME",lattice.volume());
01573 CALLBACKLIST("CELL_A",lattice.a());
01574 CALLBACKLIST("CELL_B",lattice.b());
01575 CALLBACKLIST("CELL_C",lattice.c());
01576 CALLBACKLIST("CELL_O",lattice.origin());
01577 labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
01578 << lattice.b_p() << " " << lattice.c_p() << "} ";
01579 if ( simParameters->pairInteractionOn ) {
01580 CALLBACKLIST("VDW_FORCE",pairVDWForce);
01581 CALLBACKLIST("ELECT_FORCE",pairElectForce);
01582 }
01583
01584 labels << '\0'; values << '\0'; // insane but makes Linux work
01585 std::string labelstring = labels.str();
01586 std::string valuestring = values.str();
01587 node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
01588 }
01589 #undef CALLBACKDATA
01590 #endif
01591
01592 drudeComTempAvg += drudeComTemp;
01593 drudeBondTempAvg += drudeBondTemp;
01594
01595 temp_avg += temperature;
01596 pressure_avg += trace(pressure)/3.;
01597 groupPressure_avg += trace(groupPressure)/3.;
01598 avg_count += 1;
01599
01600 if ( simParams->outputPairlists && pairlistWarnings &&
01601 ! (step % simParams->outputPairlists) ) {
01602 iout << iINFO << pairlistWarnings <<
01603 " pairlist warnings in past " << simParams->outputPairlists <<
01604 " steps.\n" << endi;
01605 pairlistWarnings = 0;
01606 }
01607
01608 // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
01609 if ( ! minimize && step % simParameters->outputEnergies ) return;
01610 // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
01611
01612 if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
01613 IMDEnergies energies;
01614 energies.tstep = step;
01615 energies.T = temp_avg/avg_count;
01616 energies.Etot = totalEnergy;
01617 energies.Epot = potentialEnergy;
01618 energies.Evdw = ljEnergy;
01619 energies.Eelec = electEnergy + electEnergySlow;
01620 energies.Ebond = bondEnergy;
01621 energies.Eangle = angleEnergy;
01622 energies.Edihe = dihedralEnergy + crosstermEnergy;
01623 energies.Eimpr = improperEnergy;
01624 Node::Object()->imd->gather_energies(&energies);
01625 }
01626
01627 if ( marginViolations ) {
01628 iout << iERROR << marginViolations <<
01629 " margin violations detected since previous energy output.\n" << endi;
01630 }
01631 marginViolations = 0;
01632
01633 if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
01634 {
01635 iout << "ETITLE: TS";
01636 iout << FORMAT("BOND");
01637 iout << FORMAT("ANGLE");
01638 iout << FORMAT("DIHED");
01639 if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
01640 iout << FORMAT("IMPRP");
01641 iout << " ";
01642 iout << FORMAT("ELECT");
01643 iout << FORMAT("VDW");
01644 iout << FORMAT("BOUNDARY");
01645 iout << FORMAT("MISC");
01646 iout << FORMAT("KINETIC");
01647 iout << " ";
01648 iout << FORMAT("TOTAL");
01649 iout << FORMAT("TEMP");
01650 iout << FORMAT("POTENTIAL");
01651 // iout << FORMAT("TOTAL2");
01652 iout << FORMAT("TOTAL3");
01653 iout << FORMAT("TEMPAVG");
01654 if ( volume != 0. ) {
01655 iout << " ";
01656 iout << FORMAT("PRESSURE");
01657 iout << FORMAT("GPRESSURE");
01658 iout << FORMAT("VOLUME");
01659 iout << FORMAT("PRESSAVG");
01660 iout << FORMAT("GPRESSAVG");
01661 }
01662 if (simParameters->drudeOn) {
01663 iout << " ";
01664 iout << FORMAT("DRUDECOM");
01665 iout << FORMAT("DRUDEBOND");
01666 iout << FORMAT("DRCOMAVG");
01667 iout << FORMAT("DRBONDAVG");
01668 }
01669 iout << "\n\n" << endi;
01670 }
01671
01672 // N.B. HP's aCC compiler merges FORMAT calls in the same expression.
01673 // Need separate statements because data returned in static array.
01674
01675 iout << ETITLE(step);
01676 iout << FORMAT(bondEnergy);
01677 iout << FORMAT(angleEnergy);
01678 if ( simParameters->mergeCrossterms ) {
01679 iout << FORMAT(dihedralEnergy+crosstermEnergy);
01680 } else {
01681 iout << FORMAT(dihedralEnergy);
01682 iout << FORMAT(crosstermEnergy);
01683 }
01684 iout << FORMAT(improperEnergy);
01685 iout << " ";
01686 iout << FORMAT(electEnergy+electEnergySlow);
01687 iout << FORMAT(ljEnergy);
01688 iout << FORMAT(boundaryEnergy);
01689 iout << FORMAT(miscEnergy);
01690 iout << FORMAT(kineticEnergy);
01691 iout << " ";
01692 iout << FORMAT(totalEnergy);
01693 iout << FORMAT(temperature);
01694 iout << FORMAT(potentialEnergy);
01695 // iout << FORMAT(flatEnergy);
01696 iout << FORMAT(smoothEnergy);
01697 iout << FORMAT(temp_avg/avg_count);
01698 if ( volume != 0. )
01699 {
01700 iout << " ";
01701 iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
01702 iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
01703 iout << FORMAT(volume);
01704 iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
01705 iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
01706 }
01707 if (simParameters->drudeOn) {
01708 iout << " ";
01709 iout << FORMAT(drudeComTemp);
01710 iout << FORMAT(drudeBondTemp);
01711 iout << FORMAT(drudeComTempAvg/avg_count);
01712 iout << FORMAT(drudeBondTempAvg/avg_count);
01713 }
01714 iout << "\n\n" << endi;
01715
01716 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
01717 char webout[80];
01718 sprintf(webout,"%d %d %d %d",(int)totalEnergy,
01719 (int)(potentialEnergy),
01720 (int)kineticEnergy,(int)temperature);
01721 CApplicationDepositNode0Data(webout);
01722 #endif
01723
01724 if (simParameters->pairInteractionOn) {
01725 iout << "PAIR INTERACTION:";
01726 iout << " STEP: " << step;
01727 iout << " VDW_FORCE: ";
01728 iout << FORMAT(pairVDWForce.x);
01729 iout << FORMAT(pairVDWForce.y);
01730 iout << FORMAT(pairVDWForce.z);
01731 iout << " ELECT_FORCE: ";
01732 iout << FORMAT(pairElectForce.x);
01733 iout << FORMAT(pairElectForce.y);
01734 iout << FORMAT(pairElectForce.z);
01735 iout << "\n" << endi;
01736 }
01737 drudeComTempAvg = 0;
01738 drudeBondTempAvg = 0;
01739 temp_avg = 0;
01740 pressure_avg = 0;
01741 groupPressure_avg = 0;
01742 avg_count = 0;
01743
01744 if ( fflush_count ) {
01745 --fflush_count;
01746 fflush(stdout);
01747 }
01748 }
|
|
|
Definition at line 875 of file Controller.C. References SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchLambda, SimParameters::alchLambda2, SimParameters::alchTemp, BigReal, iout, and simParams. Referenced by integrate(). 00876 {
00877 if (simParams->alchFepOn) {
00878 const BigReal alchLambda = simParams->alchLambda;
00879 const BigReal alchLambda2 = simParams->alchLambda2;
00880 const BigReal alchTemp = simParams->alchTemp;
00881 const int alchEquilSteps = simParams->alchEquilSteps;
00882 iout << "FEP: RESETTING FOR NEW FEP WINDOW "
00883 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2
00884 << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
00885 << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
00886 << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
00887 << " K FOR FEP CALCULATION\n" << endi;
00888 }
00889 }
|
|
|
Definition at line 1308 of file Controller.C. References BigReal, compareChecksums(), iout, RequireReduction::item(), iWARN(), min_energy, min_f_dot_f, min_f_dot_v, min_huge_count, min_v_dot_v, Node::molecule, Molecule::numCalcExclusions, Node::Object(), printEnergies(), receivePressure(), reduction, REDUCTION_EXCLUSION_CHECKSUM, REDUCTION_MIN_F_DOT_F, REDUCTION_MIN_F_DOT_V, REDUCTION_MIN_HUGE_COUNT, and REDUCTION_MIN_V_DOT_V. 01308 {
01309
01310 receivePressure(step,1);
01311
01312 Node *node = Node::Object();
01313 Molecule *molecule = node->molecule;
01314 BigReal checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01315 if ( ((int)checksum) && ((int)checksum) < molecule->numCalcExclusions )
01316 iout << iWARN << "Bad global exclusion count, possible error!\n" << iWARN
01317 << "Increasing cutoff during minimization may avoid this.\n" << endi;
01318 compareChecksums(step,1);
01319
01320 printEnergies(step,1);
01321
01322 min_energy = totalEnergy;
01323 min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
01324 min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
01325 min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
01326 min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
01327
01328 }
|
|
|
Definition at line 890 of file Controller.C. References SimParameters::alchEquilSteps, SimParameters::alchLambda, SimParameters::alchTemp, SimParameters::alchThermIntOn, BigReal, iout, and simParams. Referenced by integrate(). 00891 {
00892 if (simParams->alchThermIntOn) {
00893 const BigReal alchLambda = simParams->alchLambda;
00894 const BigReal alchTemp = simParams->alchTemp;
00895 const int alchEquilSteps = simParams->alchEquilSteps;
00896 iout << "TI: RESETTING FOR NEW WINDOW "
00897 << "LAMBDA SET TO " << alchLambda
00898 << "\nTI: WINDOW TO HAVE " << alchEquilSteps
00899 << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
00900 }
00901 }
|
|
|
Definition at line 1275 of file Controller.C. References SimParameters::firstTimestep, iout, iWARN(), memusage_MB(), SimParameters::N, SimParameters::outputEnergies, SimParameters::outputTiming, and simParams. Referenced by printEnergies(). 01275 {
01276
01277 if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
01278 {
01279 const double endWTime = CmiWallTimer() - firstWTime;
01280 const double endCTime = CmiTimer() - firstCTime;
01281
01282 const double elapsedW =
01283 (endWTime - startWTime) / simParams->outputTiming;
01284 const double elapsedC =
01285 (endCTime - startCTime) / simParams->outputTiming;
01286
01287 const double remainingW = elapsedW * (simParams->N - step);
01288 const double remainingW_hours = remainingW / 3600;
01289
01290 startWTime = endWTime;
01291 startCTime = endCTime;
01292
01293 #ifdef NAMD_CUDA
01294 if ( simParams->outputEnergies < 60 &&
01295 step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
01296 iout << iWARN << "Energy evaluation is done on CPU, increase outputEnergies to improve performance.\n" << endi;
01297 }
01298 #endif
01299 if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
01300 CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
01301 ", %g hours remaining, %f MB of memory in use.\n",
01302 step, endCTime, elapsedC, endWTime, elapsedW,
01303 remainingW_hours, memusage_MB());
01304 }
01305 }
01306 }
|
|
|
Definition at line 904 of file Controller.C. References BigReal, iout, SimParameters::reassignFreq, SimParameters::reassignHold, SimParameters::reassignIncr, SimParameters::reassignTemp, and simParams. Referenced by integrate(). 00905 {
00906 const int reassignFreq = simParams->reassignFreq;
00907 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
00908 BigReal newTemp = simParams->reassignTemp;
00909 newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
00910 if ( simParams->reassignIncr > 0.0 ) {
00911 if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
00912 newTemp = simParams->reassignHold;
00913 } else {
00914 if ( newTemp < simParams->reassignHold )
00915 newTemp = simParams->reassignHold;
00916 }
00917 iout << "REASSIGNING VELOCITIES AT STEP " << step
00918 << " TO " << newTemp << " KELVIN.\n" << endi;
00919 }
00920 }
|
|
|
Definition at line 2062 of file Controller.C. References fflush_count, LdbCoordinator::getNumStepsToRun(), ldbSteps, LdbCoordinator::Object(), and LdbCoordinator::rebalance(). Referenced by integrate(). 02063 {
02064 if ( ! ldbSteps ) {
02065 ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
02066 }
02067 if ( ! --ldbSteps ) {
02068 startBenchTime -= CmiWallTimer();
02069 LdbCoordinator::Object()->rebalance(this);
02070 startBenchTime += CmiWallTimer();
02071 fflush_count = 3;
02072 }
02073 }
|
|
||||||||||||
|
Definition at line 957 of file Controller.C. References AVGXY, BigReal, SimParameters::comMove, controlNumDegFreedom, controlPressure, controlPressure_nbond, controlPressure_normal, controlPressure_slow, Tensor::diagonal(), drudeBondTemp, drudeComTemp, SimParameters::drudeOn, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, SimParameters::fixedAtomsOn, GET_TENSOR, GET_VECTOR, groupPressure, groupPressure_nbond, groupPressure_normal, groupPressure_slow, Tensor::identity(), iINFO(), iout, RequireReduction::item(), kineticEnergy, kineticEnergyCentered, kineticEnergyHalfstep, SimParameters::langevinOn, NamdState::lattice, SimParameters::LJcorrection, Node::molecule, Molecule::num_deg_freedom(), Molecule::num_fixed_atoms(), Molecule::num_fixed_groups(), Molecule::num_group_deg_freedom(), Molecule::numAtoms, Molecule::numConstraints, numDegFreedom, Molecule::numDrudeAtoms, Molecule::numFepInitial, Molecule::numFixedAtoms, Molecule::numFixedGroups, Molecule::numFixedRigidBonds, Molecule::numHydrogenGroups, Molecule::numLonepairs, Molecule::numRigidBonds, Node::Object(), Lattice::origin(), outer(), SimParameters::pairInteractionOn, pressure, pressure_nbond, pressure_normal, pressure_slow, reduction, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY, REDUCTION_EXT_FORCE_NBOND, REDUCTION_EXT_FORCE_NORMAL, REDUCTION_EXT_FORCE_SLOW, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, RequireReduction::require(), Node::simParameters, state, Molecule::tail_corr_virial, temperature, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, SimParameters::useGroupPressure, and Lattice::volume(). Referenced by integrate(), and printMinimizeEnergies(). 00958 {
00959 Node *node = Node::Object();
00960 Molecule *molecule = node->molecule;
00961 SimParameters *simParameters = node->simParameters;
00962 Lattice &lattice = state->lattice;
00963
00964 reduction->require();
00965
00966 Tensor virial;
00967 Tensor virial_normal;
00968 Tensor virial_nbond;
00969 Tensor virial_slow;
00970 #ifdef ALTVIRIAL
00971 Tensor altVirial_normal;
00972 Tensor altVirial_nbond;
00973 Tensor altVirial_slow;
00974 #endif
00975 Tensor intVirial;
00976 Tensor intVirial_normal;
00977 Tensor intVirial_nbond;
00978 Tensor intVirial_slow;
00979 Vector extForce_normal;
00980 Vector extForce_nbond;
00981 Vector extForce_slow;
00982 BigReal volume;
00983
00984 #if 1
00985 numDegFreedom = molecule->num_deg_freedom();
00986 int numGroupDegFreedom = molecule->num_group_deg_freedom();
00987 int numFixedGroups = molecule->num_fixed_groups();
00988 int numFixedAtoms = molecule->num_fixed_atoms();
00989 #endif
00990 #if 0
00991 int numAtoms = molecule->numAtoms;
00992 numDegFreedom = 3 * numAtoms;
00993 int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
00994 int numFixedAtoms =
00995 ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
00996 int numLonepairs = molecule->numLonepairs;
00997 int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
00998 if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
00999 if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01000 if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01001 if ( ! ( numFixedAtoms || molecule->numConstraints
01002 || simParameters->comMove || simParameters->langevinOn ) ) {
01003 numDegFreedom -= 3;
01004 numGroupDegFreedom -= 3;
01005 }
01006 if (simParameters->pairInteractionOn) {
01007 // this doesn't attempt to deal with fixed atoms or constraints
01008 numDegFreedom = 3 * molecule->numFepInitial;
01009 }
01010 int numRigidBonds = molecule->numRigidBonds;
01011 int numFixedRigidBonds =
01012 ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01013
01014 // numLonepairs is subtracted here because all lonepairs have a rigid bond
01015 // to oxygen, but all of the LP degrees of freedom are dealt with above
01016 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01017 #endif
01018
01019 kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01020 kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01021
01022 BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01023 reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01024 BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01025 reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01026
01027 BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01028 / ( numDegFreedom * BOLTZMAN );
01029 BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01030 / ( numDegFreedom * BOLTZMAN );
01031 BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01032 / ( numGroupDegFreedom * BOLTZMAN );
01033 BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01034 / ( numGroupDegFreedom * BOLTZMAN );
01035
01036 /* test code for comparing different temperatures
01037 iout << "TEMPTEST: " << step << " " <<
01038 atomTempHalfstep << " " <<
01039 atomTempCentered << " " <<
01040 groupTempHalfstep << " " <<
01041 groupTempCentered << "\n" << endi;
01042 iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
01043 numGroupDegFreedom << "\n" << endi;
01044 */
01045
01046 GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01047 GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01048 GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01049
01050 #ifdef ALTVIRIAL
01051 GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01052 GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01053 GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01054 #endif
01055
01056 GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01057 GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01058 GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01059
01060 GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01061 GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01062 GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01063 Vector extPosition = lattice.origin();
01064 virial_normal -= outer(extForce_normal,extPosition);
01065 virial_nbond -= outer(extForce_nbond,extPosition);
01066 virial_slow -= outer(extForce_slow,extPosition);
01067
01068
01069 kineticEnergy = kineticEnergyCentered;
01070 temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMAN );
01071
01072 if (simParameters->drudeOn) {
01073 BigReal drudeComKE
01074 = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01075 BigReal drudeBondKE
01076 = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01077 int g_bond = 3 * molecule->numDrudeAtoms;
01078 int g_com = numDegFreedom - g_bond;
01079 drudeComTemp = 2.0 * drudeComKE / (g_com * BOLTZMAN);
01080 drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMAN)) : 0.);
01081 }
01082
01083 if ( (volume=lattice.volume()) != 0. )
01084 {
01085
01086 if (simParameters->LJcorrection && volume) {
01087 // Apply tail correction to pressure
01088 //printf("Volume is %f\n", volume);
01089 //printf("Applying tail correction of %f to virial\n", molecule->tail_corr_virial / volume);
01090 virial_normal += Tensor::identity(molecule->tail_corr_virial / volume);
01091 }
01092
01093 // kinetic energy component included in virials
01094 pressure_normal = virial_normal / volume;
01095 groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01096
01097 if ( minimize || ! ( step % nbondFreq ) )
01098 {
01099 pressure_nbond = virial_nbond / volume;
01100 groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01101 }
01102
01103 if ( minimize || ! ( step % slowFreq ) )
01104 {
01105 pressure_slow = virial_slow / volume;
01106 groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01107 }
01108
01109 /*
01110 iout << "VIRIALS: " << virial_normal << " " << virial_nbond << " " <<
01111 virial_slow << " " << ( virial_normal - intVirial_normal ) << " " <<
01112 ( virial_nbond - intVirial_nbond ) << " " <<
01113 ( virial_slow - intVirial_slow ) << "\n";
01114 */
01115
01116 pressure = pressure_normal + pressure_nbond + pressure_slow;
01117 groupPressure = groupPressure_normal + groupPressure_nbond +
01118 groupPressure_slow;
01119 }
01120 else
01121 {
01122 pressure = Tensor();
01123 groupPressure = Tensor();
01124 }
01125
01126 if ( simParameters->useGroupPressure )
01127 {
01128 controlPressure_normal = groupPressure_normal;
01129 controlPressure_nbond = groupPressure_nbond;
01130 controlPressure_slow = groupPressure_slow;
01131 controlPressure = groupPressure;
01132 controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01133 if ( ! ( numFixedAtoms || molecule->numConstraints
01134 || simParameters->comMove || simParameters->langevinOn ) ) {
01135 controlNumDegFreedom -= 1;
01136 }
01137 }
01138 else
01139 {
01140 controlPressure_normal = pressure_normal;
01141 controlPressure_nbond = pressure_nbond;
01142 controlPressure_slow = pressure_slow;
01143 controlPressure = pressure;
01144 controlNumDegFreedom = numDegFreedom / 3;
01145 }
01146
01147 if (simParameters->fixCellDims) {
01148 if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01149 if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01150 if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01151 }
01152
01153 if ( simParameters->useFlexibleCell ) {
01154 // use symmetric pressure to control rotation
01155 // controlPressure_normal = symmetric(controlPressure_normal);
01156 // controlPressure_nbond = symmetric(controlPressure_nbond);
01157 // controlPressure_slow = symmetric(controlPressure_slow);
01158 // controlPressure = symmetric(controlPressure);
01159 // only use on-diagonal components for now
01160 controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01161 controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01162 controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01163 controlPressure = Tensor::diagonal(diagonal(controlPressure));
01164 if ( simParameters->useConstantRatio ) {
01165 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01166 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01167 AVGXY(controlPressure_normal);
01168 AVGXY(controlPressure_nbond);
01169 AVGXY(controlPressure_slow);
01170 AVGXY(controlPressure);
01171 #undef AVGXY
01172 }
01173 } else {
01174 controlPressure_normal =
01175 Tensor::identity(trace(controlPressure_normal)/3.);
01176 controlPressure_nbond =
01177 Tensor::identity(trace(controlPressure_nbond)/3.);
01178 controlPressure_slow =
01179 Tensor::identity(trace(controlPressure_slow)/3.);
01180 controlPressure =
01181 Tensor::identity(trace(controlPressure)/3.);
01182 }
01183
01184 #ifdef DEBUG_PRESSURE
01185 iout << iINFO << "Control pressure = " << controlPressure <<
01186 " = " << controlPressure_normal << " + " <<
01187 controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01188 #endif
01189
01190 }
|
|
|
Definition at line 838 of file Controller.C. References BigReal, broadcast, iout, SimpleBroadcastObject< T >::publish(), SimParameters::rescaleFreq, SimParameters::rescaleTemp, rescaleVelocities_numTemps, rescaleVelocities_sumTemps, simParams, and ControllerBroadcasts::velocityRescaleFactor. Referenced by integrate(). 00839 {
00840 const int rescaleFreq = simParams->rescaleFreq;
00841 if ( rescaleFreq > 0 ) {
00842 rescaleVelocities_sumTemps += temperature; ++rescaleVelocities_numTemps;
00843 if ( rescaleVelocities_numTemps == rescaleFreq ) {
00844 BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
00845 const BigReal rescaleTemp = simParams->rescaleTemp;
00846 BigReal factor = sqrt(rescaleTemp/avgTemp);
00847 broadcast->velocityRescaleFactor.publish(step,factor);
00848 iout << "RESCALING VELOCITIES AT STEP " << step
00849 << " FROM AVERAGE TEMPERATURE OF " << avgTemp
00850 << " TO " << rescaleTemp << " KELVIN.\n" << endi;
00851 rescaleVelocities_sumTemps = 0; rescaleVelocities_numTemps = 0;
00852 }
00853 }
00854 }
|
|
|
Definition at line 216 of file Controller.C. References awaken(), CTRL_STK_SZ, and DebugM. Referenced by NamdState::runController(). 00217 {
00218 // create a Thread and invoke it
00219 DebugM(4, "Starting thread in controller on this=" << this << "\n");
00220 thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
00221 CthSetStrategyDefault(thread);
00222 #if CMK_BLUEGENE_CHARM
00223 BgAttach(thread);
00224 #endif
00225 awaken();
00226 }
|
|
|
Definition at line 922 of file Controller.C. References BigReal, broadcast, SimpleBroadcastObject< T >::publish(), simParams, ControllerBroadcasts::tcoupleCoefficient, SimParameters::tCoupleOn, SimParameters::tCoupleTemp, and temperature. Referenced by integrate(). 00923 {
00924 if ( simParams->tCoupleOn )
00925 {
00926 const BigReal tCoupleTemp = simParams->tCoupleTemp;
00927 BigReal coefficient = 1.;
00928 if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
00929 broadcast->tcoupleCoefficient.publish(step,coefficient);
00930 }
00931 }
|
|
|
Definition at line 2085 of file Controller.C. References BackEnd::awaken(). Referenced by algorithm(). 02085 {
02086 BackEnd::awaken();
02087 CthFree(thread);
02088 CthSuspend();
02089 }
|
|
||||||||||||
|
Definition at line 1763 of file Controller.C. References Lattice::a(), Lattice::a_p(), Lattice::b(), Lattice::b_p(), Lattice::c(), Lattice::c_p(), langevinPiston_strainRate, SimParameters::langevinPistonOn, NamdState::lattice, Lattice::origin(), simParams, state, Vector::x, Vector::y, and Vector::z. Referenced by outputExtendedSystem(). 01763 {
01764 Lattice &lattice = state->lattice;
01765 file << step;
01766 if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
01767 if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
01768 if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
01769 file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
01770 if ( simParams->langevinPistonOn ) {
01771 Vector strainRate = diagonal(langevinPiston_strainRate);
01772 Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
01773 file << " " << strainRate.x;
01774 file << " " << strainRate.y;
01775 file << " " << strainRate.z;
01776 file << " " << strainRate2.x;
01777 file << " " << strainRate2.y;
01778 file << " " << strainRate2.z;
01779 }
01780 file << std::endl;
01781 }
|
|
|
Definition at line 1750 of file Controller.C. References Lattice::a_p(), Lattice::b_p(), Lattice::c_p(), SimParameters::langevinPistonOn, NamdState::lattice, simParams, and state. Referenced by outputExtendedSystem(). 01750 {
01751 Lattice &lattice = state->lattice;
01752 file << "#$LABELS step";
01753 if ( lattice.a_p() ) file << " a_x a_y a_z";
01754 if ( lattice.b_p() ) file << " b_x b_y b_z";
01755 if ( lattice.c_p() ) file << " c_x c_y c_z";
01756 file << " o_x o_y o_z";
01757 if ( simParams->langevinPistonOn ) {
01758 file << " s_x s_y s_z s_u s_v s_w";
01759 }
01760 file << std::endl;
01761 }
|
|
||||||||||||
|
Definition at line 1939 of file Controller.C. References SimParameters::alchTemp, BigReal, BOLTZMAN, dG, electEnergy, electEnergy_f, exp_dE_ByRT, fepFile, FepNo, FEPTITLE(), FORMAT(), ljEnergy_f, net_dE, simParams, and temperature. Referenced by outputFepEnergy(). 01939 {
01940 BigReal eeng = electEnergy+electEnergySlow;
01941 BigReal eeng_f = electEnergy_f + electEnergySlow_f;
01942 BigReal dE = eeng_f + ljEnergy_f - eeng - ljEnergy;
01943 BigReal RT = BOLTZMAN * simParams->alchTemp;
01944 dG = -(RT * log(exp_dE_ByRT/FepNo));
01945 BigReal dE_avg = net_dE/FepNo;
01946 fepFile << FEPTITLE(step);
01947 fepFile << FORMAT(eeng);
01948 fepFile << FORMAT(eeng_f);
01949 fepFile << FORMAT(ljEnergy);
01950 fepFile << FORMAT(ljEnergy_f);
01951 fepFile << FORMAT(dE);
01952 fepFile << FORMAT(dE_avg);
01953 fepFile << FORMAT(temperature);
01954 fepFile << FORMAT(dG);
01955 fepFile << std::endl;
01956 }
|
|
||||||||||||
|
Definition at line 1958 of file Controller.C. References FORMAT(), net_dEdl_elec_1, net_dEdl_elec_2, net_dEdl_lj_1, net_dEdl_lj_2, recent_dEdl_elec_1, recent_dEdl_elec_2, recent_dEdl_lj_1, recent_dEdl_lj_2, recent_TiNo, tiFile, TiNo, and TITITLE(). Referenced by outputTiEnergy(). 01958 {
01959 tiFile << TITITLE(step);
01960 tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
01961 tiFile << FORMAT(net_dEdl_elec_1/TiNo);
01962 tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
01963 tiFile << FORMAT(net_dEdl_lj_1/TiNo);
01964 tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
01965 tiFile << FORMAT(net_dEdl_elec_2/TiNo);
01966 tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
01967 tiFile << FORMAT(net_dEdl_lj_2/TiNo);
01968 tiFile << std::endl;
01969 }
|
|
|
Definition at line 38 of file Controller.h. |
|
|
Definition at line 59 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 136 of file Controller.h. Referenced by algorithm(), berendsenPressure(), and Controller(). |
|
|
Definition at line 137 of file Controller.h. Referenced by algorithm(), berendsenPressure(), and Controller(). |
|
|
Definition at line 168 of file Controller.h. Referenced by algorithm(), berendsenPressure(), Controller(), correctMomentum(), cycleBarrier(), langevinPiston1(), minimize(), rescaleVelocities(), and tcoupleVelocities(). |
|
|
Definition at line 187 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 188 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 186 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 185 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 184 of file Controller.h. Referenced by algorithm(), and Controller(). |
|
|
Definition at line 166 of file Controller.h. Referenced by enqueueCollections(). |
|
|
Definition at line 64 of file Controller.h. Referenced by compareChecksums(). |
|
|
Definition at line 126 of file Controller.h. Referenced by langevinPiston1(), langevinPiston2(), and receivePressure(). |
|
|
Definition at line 127 of file Controller.h. Referenced by langevinPiston1(), langevinPiston2(), and receivePressure(). |
|
|
Definition at line 52 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 51 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 53 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 87 of file Controller.h. Referenced by writeFepEnergyData(). |
|
|
Definition at line 114 of file Controller.h. Referenced by Controller(), printEnergies(), and receivePressure(). |
|
|
Definition at line 116 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 113 of file Controller.h. Referenced by Controller(), printEnergies(), and receivePressure(). |
|
|
Definition at line 115 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 78 of file Controller.h. Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData(). |
|
|
Definition at line 82 of file Controller.h. Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData(). |
|
|
Definition at line 93 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 96 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 103 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 104 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 79 of file Controller.h. Referenced by outputFepEnergy(), and printEnergies(). |
|
|
Definition at line 83 of file Controller.h. Referenced by outputFepEnergy(), and printEnergies(). |
|
|
Definition at line 94 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 97 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 85 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 175 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 88 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 90 of file Controller.h. Referenced by outputFepEnergy(). |
|
|
Definition at line 144 of file Controller.h. Referenced by rebalanceLoad(). |
|
|
Definition at line 125 of file Controller.h. Referenced by printEnergies(), and receivePressure(). |
|
|
Definition at line 58 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 49 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 48 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 50 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 61 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 118 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 120 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 119 of file Controller.h. Referenced by printEnergies(), and receivePressure(). |
|
|
Definition at line 140 of file Controller.h. Referenced by algorithm(), Controller(), and writeExtendedSystemData(). |
|
|
Definition at line 142 of file Controller.h. Referenced by rebalanceLoad(). |
|
|
Definition at line 80 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 84 of file Controller.h. Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData(). |
|
|
Definition at line 95 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 98 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 65 of file Controller.h. Referenced by compareChecksums(), and printEnergies(). |
|
|
Definition at line 69 of file Controller.h. Referenced by printMinimizeEnergies(). |
|
|
Definition at line 70 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 71 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 73 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 72 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 54 of file Controller.h. Referenced by integrate(), langevinPiston1(), langevinPiston2(), and minimize(). |
|
|
Definition at line 86 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 99 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 100 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 101 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 102 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 76 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 66 of file Controller.h. Referenced by compareChecksums(), and printEnergies(). |
|
|
Definition at line 156 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 158 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 157 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 124 of file Controller.h. Referenced by printEnergies(), and receivePressure(). |
|
|
Definition at line 57 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 46 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 45 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 47 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 60 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 161 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 160 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 159 of file Controller.h. Referenced by Controller(). |
|
|
Definition at line 150 of file Controller.h. Referenced by Controller(), langevinPiston1(), and langevinPiston2(). |
|
|
Definition at line 106 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 107 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 108 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 109 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 110 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 153 of file Controller.h. Referenced by compareChecksums(), Controller(), correctMomentum(), printDynamicsEnergies(), printEnergies(), printMinimizeEnergies(), and receivePressure(). |
|
|
Definition at line 132 of file Controller.h. Referenced by Controller(), and rescaleVelocities(). |
|
|
Definition at line 131 of file Controller.h. Referenced by Controller(), and rescaleVelocities(). |
|
|
|
Definition at line 55 of file Controller.h. Referenced by correctMomentum(), integrate(), langevinPiston1(), langevinPiston2(), and minimize(). |
|
|
Definition at line 122 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 123 of file Controller.h. |
|
|
Definition at line 152 of file Controller.h. Referenced by algorithm(), berendsenPressure(), enqueueCollections(), langevinPiston1(), langevinPiston2(), printDynamicsEnergies(), printEnergies(), receivePressure(), writeExtendedSystemData(), and writeExtendedSystemLabels(). |
|
|
Definition at line 62 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 56 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 121 of file Controller.h. Referenced by printEnergies(), receivePressure(), tcoupleVelocities(), and writeFepEnergyData(). |
|
|
Definition at line 179 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 105 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 77 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 169 of file Controller.h. Referenced by outputExtendedSystem(). |
1.3.9.1