00001
00007
00008
00009
00010
00011
00012
00013
00014 #include "InfoStream.h"
00015 #include "memusage.h"
00016 #include "Node.h"
00017 #include "Molecule.h"
00018 #include "SimParameters.h"
00019 #include "Controller.h"
00020 #include "ReductionMgr.h"
00021 #include "CollectionMaster.h"
00022 #include "Output.h"
00023 #include "strlib.h"
00024 #include "BroadcastObject.h"
00025 #include "NamdState.h"
00026 #include "ScriptTcl.h"
00027 #include "Broadcasts.h"
00028 #include "LdbCoordinator.h"
00029 #include "Thread.h"
00030 #include <math.h>
00031 #include <signal.h>
00032 #include "NamdOneTools.h"
00033 #include "PatchMap.h"
00034 #include "PatchMap.inl"
00035 #include "Random.h"
00036 #include "imd.h"
00037 #include "IMDOutput.h"
00038 #include "BackEnd.h"
00039 #include <fstream>
00040 #include <iomanip>
00041 #include <errno.h>
00042 #include "qd.h"
00043
00044 #include "ComputeNonbondedMICKernel.h"
00045
00046 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
00047 extern "C" void CApplicationDepositNode0Data(char *);
00048 #endif
00049
00050 #ifndef cbrt
00051
00052 #define cbrt(x) pow(x,(double)(1.0/3.0))
00053 #endif
00054
00055
00056 #define MIN_DEBUG_LEVEL 3
00057
00058 #include "Debug.h"
00059
00060 #define XXXBIGREAL 1.0e32
00061
00062 class PressureProfileReduction {
00063 private:
00064 RequireReduction *reduction;
00065 const int nslabs;
00066 const int freq;
00067 int nelements;
00068 int count;
00069 char *name;
00070
00071 BigReal *average;
00072
00073 public:
00074 PressureProfileReduction(int rtag, int numslabs, int numpartitions,
00075 const char *myname, int outputfreq)
00076 : nslabs(numslabs), freq(outputfreq) {
00077 name = strdup(myname);
00078 nelements = 3*nslabs * numpartitions;
00079 reduction = ReductionMgr::Object()->willRequire(rtag,nelements);
00080
00081 average = new BigReal[nelements];
00082 count = 0;
00083 }
00084 ~PressureProfileReduction() {
00085 delete [] average;
00086 delete reduction;
00087 free(name);
00088 }
00089
00090 void getData(int firsttimestep, int step, const Lattice &lattice,
00091 BigReal *total) {
00092 reduction->require();
00093
00094 int i;
00095 double inv_volume = 1.0 / lattice.volume();
00096
00097 int arraysize = 3*nslabs;
00098 for (i=0; i<nelements; i++) {
00099 BigReal val = reduction->item(i) * inv_volume;
00100 average[i] += val;
00101 total[i % arraysize] += val;
00102 }
00103 count++;
00104 if (!(step % freq)) {
00105
00106 BigReal scalefac = PRESSUREFACTOR * nslabs;
00107
00108 iout << "PPROFILE" << name << ": " << step << " ";
00109 if (step == firsttimestep) {
00110
00111 for (i=0; i<nelements; i++)
00112 iout << reduction->item(i)*scalefac*inv_volume << " ";
00113 } else {
00114
00115 scalefac /= count;
00116 for (i=0; i<nelements; i++)
00117 iout << average[i]*scalefac << " ";
00118 }
00119 iout << "\n" << endi;
00120
00121 memset(average, 0, nelements*sizeof(BigReal));
00122 count = 0;
00123 }
00124 }
00125 };
00126
00127
00128 Controller::Controller(NamdState *s) :
00129 computeChecksum(0), marginViolations(0), pairlistWarnings(0),
00130 simParams(Node::Object()->simParameters),
00131 state(s),
00132 collection(CollectionMaster::Object()),
00133 startCTime(0),
00134 firstCTime(CmiTimer()),
00135 startWTime(0),
00136 firstWTime(CmiWallTimer()),
00137 startBenchTime(0),
00138 stepInFullRun(0),
00139 ldbSteps(0),
00140 fflush_count(3)
00141 {
00142 broadcast = new ControllerBroadcasts;
00143 min_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MINIMIZER,1);
00144
00145 if (simParams->accelMDOn) {
00146
00147 amd_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
00148
00149 submit_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_AMD);
00150
00151 reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_AMD);
00152 } else {
00153 amd_reduction = NULL;
00154 submit_reduction = NULL;
00155 reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
00156 }
00157
00158 pressureProfileSlabs = 0;
00159 pressureProfileCount = 0;
00160 ppbonded = ppnonbonded = ppint = NULL;
00161 if (simParams->pressureProfileOn || simParams->pressureProfileEwaldOn) {
00162 int ntypes = simParams->pressureProfileAtomTypes;
00163 int nslabs = pressureProfileSlabs = simParams->pressureProfileSlabs;
00164
00165 int npairs = (ntypes * (ntypes+1))/2;
00166 pressureProfileAverage = new BigReal[3*nslabs];
00167 memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
00168 int freq = simParams->pressureProfileFreq;
00169 if (simParams->pressureProfileOn) {
00170 ppbonded = new PressureProfileReduction(REDUCTIONS_PPROF_BONDED,
00171 nslabs, npairs, "BONDED", freq);
00172 ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED,
00173 nslabs, npairs, "NONBONDED", freq);
00174
00175 ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL,
00176 nslabs, ntypes, "INTERNAL", freq);
00177 } else {
00178
00179 ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED,
00180 nslabs, npairs, "NONBONDED", freq);
00181 }
00182 }
00183 random = new Random(simParams->randomSeed);
00184 random->split(0,PatchMap::Object()->numPatches()+1);
00185
00186 heat = totalEnergy0 = 0.0;
00187 rescaleVelocities_sumTemps = 0;
00188 rescaleVelocities_numTemps = 0;
00189 stochRescale_count = 0;
00190 if (simParams->stochRescaleOn) {
00191 stochRescaleTimefactor = exp(-simParams->stochRescaleFreq *
00192 simParams->dt * 0.001 / simParams->stochRescalePeriod);
00193 }
00194 berendsenPressure_avg = 0; berendsenPressure_count = 0;
00195
00196 langevinPiston_strainRate =
00197 Tensor::symmetric(simParams->strainRate,simParams->strainRate2);
00198 if ( ! simParams->useFlexibleCell ) {
00199 BigReal avg = trace(langevinPiston_strainRate) / 3.;
00200 langevinPiston_strainRate = Tensor::identity(avg);
00201 } else if ( simParams->useConstantRatio ) {
00202 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
00203 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
00204 AVGXY(langevinPiston_strainRate);
00205 #undef AVGXY
00206 }
00207 langevinPiston_origStrainRate = langevinPiston_strainRate;
00208 if (simParams->multigratorOn) {
00209 multigratorXi = 0.0;
00210 int n = simParams->multigratorNoseHooverChainLength;
00211 BigReal tau = simParams->multigratorTemperatureRelaxationTime;
00212 Node *node = Node::Object();
00213 Molecule *molecule = node->molecule;
00214 BigReal Nf = molecule->num_deg_freedom();
00215 BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00216 multigratorNu.resize(n);
00217 multigratorNuT.resize(n);
00218 multigratorZeta.resize(n);
00219 multigratorOmega.resize(n);
00220 for (int i=0;i < n;i++) {
00221 multigratorNu[i] = 0.0;
00222 multigratorZeta[i] = 0.0;
00223 if (i == 0) {
00224 multigratorOmega[i] = Nf*kT0*tau*tau;
00225 } else {
00226 multigratorOmega[i] = kT0*tau*tau;
00227 }
00228 }
00229 multigratorReduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
00230 } else {
00231 multigratorReduction = NULL;
00232 }
00233 origLattice = state->lattice;
00234 smooth2_avg = XXXBIGREAL;
00235 temp_avg = 0;
00236 pressure_avg = 0;
00237 groupPressure_avg = 0;
00238 avg_count = 0;
00239 pressure_tavg = 0;
00240 groupPressure_tavg = 0;
00241 tavg_count = 0;
00242 checkpoint_stored = 0;
00243 drudeBondTemp = 0;
00244 drudeBondTempAvg = 0;
00245 cumAlchWork = 0;
00246 }
00247
00248 Controller::~Controller(void)
00249 {
00250 delete broadcast;
00251 delete reduction;
00252 delete min_reduction;
00253 delete amd_reduction;
00254 delete submit_reduction;
00255 delete ppbonded;
00256 delete ppnonbonded;
00257 delete ppint;
00258 delete [] pressureProfileAverage;
00259 delete random;
00260 if (multigratorReduction) delete multigratorReduction;
00261 }
00262
00263 void Controller::threadRun(Controller* arg)
00264 {
00265 arg->algorithm();
00266 }
00267
00268 void Controller::run(void)
00269 {
00270
00271 DebugM(4, "Starting thread in controller on this=" << this << "\n");
00272 thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
00273 CthSetStrategyDefault(thread);
00274 #if CMK_BLUEGENE_CHARM
00275 BgAttach(thread);
00276 #endif
00277 awaken();
00278 }
00279
00280
00281 void Controller::algorithm(void)
00282 {
00283 int scriptTask;
00284 int scriptSeq = 0;
00285 BackEnd::awaken();
00286 while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
00287 switch ( scriptTask ) {
00288 case SCRIPT_OUTPUT:
00289 enqueueCollections(FILE_OUTPUT);
00290 outputExtendedSystem(FILE_OUTPUT);
00291 break;
00292 case SCRIPT_FORCEOUTPUT:
00293 enqueueCollections(FORCE_OUTPUT);
00294 break;
00295 case SCRIPT_MEASURE:
00296 enqueueCollections(EVAL_MEASURE);
00297 break;
00298 case SCRIPT_REINITVELS:
00299 iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
00300 << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
00301 break;
00302 case SCRIPT_RESCALEVELS:
00303 iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
00304 << " BY " << simParams->scriptArg1 << "\n" << endi;
00305 break;
00306 case SCRIPT_RESCALESOLUTECHARGES:
00307
00308
00309 break;
00310 case SCRIPT_CHECKPOINT:
00311 iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
00312 << "\n" << endi;
00313 checkpoint_stored = 1;
00314 checkpoint_lattice = state->lattice;
00315 checkpoint_state = *(ControllerState*)this;
00316 break;
00317 case SCRIPT_REVERT:
00318 iout << "REVERTING AT STEP " << simParams->firstTimestep
00319 << "\n" << endi;
00320 if ( ! checkpoint_stored )
00321 NAMD_die("Unable to revert, checkpoint was never called!");
00322 state->lattice = checkpoint_lattice;
00323 *(ControllerState*)this = checkpoint_state;
00324 break;
00325 case SCRIPT_CHECKPOINT_STORE:
00326 iout << "STORING CHECKPOINT AT STEP " << simParams->firstTimestep
00327 << " TO KEY " << simParams->scriptStringArg1;
00328 if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00329 iout << "\n" << endi;
00330 if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00331 if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00332 checkpoints[simParams->scriptStringArg1] = new checkpoint;
00333 }
00334 checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
00335 cp.lattice = state->lattice;
00336 cp.state = *(ControllerState*)this;
00337 } else {
00338 Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00339 scriptTask, state->lattice, *(ControllerState*)this);
00340 }
00341 break;
00342 case SCRIPT_CHECKPOINT_LOAD:
00343 iout << "LOADING CHECKPOINT AT STEP " << simParams->firstTimestep
00344 << " FROM KEY " << simParams->scriptStringArg1;
00345 if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00346 iout << "\n" << endi;
00347 if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00348 if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00349 NAMD_die("Unable to load checkpoint, requested key was never stored.");
00350 }
00351 checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
00352 state->lattice = cp.lattice;
00353 *(ControllerState*)this = cp.state;
00354 } else {
00355 Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00356 scriptTask, state->lattice, *(ControllerState*)this);
00357 }
00358 break;
00359 case SCRIPT_CHECKPOINT_SWAP:
00360 iout << "SWAPPING CHECKPOINT AT STEP " << simParams->firstTimestep
00361 << " FROM KEY " << simParams->scriptStringArg1;
00362 if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00363 iout << "\n" << endi;
00364 if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00365 if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00366 NAMD_die("Unable to swap checkpoint, requested key was never stored.");
00367 }
00368 checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
00369 std::swap(state->lattice,cp.lattice);
00370 std::swap(*(ControllerState*)this,cp.state);
00371 } else {
00372 Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00373 scriptTask, state->lattice, *(ControllerState*)this);
00374 }
00375 break;
00376 case SCRIPT_CHECKPOINT_FREE:
00377 iout << "FREEING CHECKPOINT AT STEP " << simParams->firstTimestep
00378 << " FROM KEY " << simParams->scriptStringArg1;
00379 if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00380 iout << "\n" << endi;
00381 if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00382 if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00383 NAMD_die("Unable to free checkpoint, requested key was never stored.");
00384 }
00385 delete checkpoints[simParams->scriptStringArg1];
00386 checkpoints.erase(simParams->scriptStringArg1);
00387 } else {
00388 Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00389 scriptTask, state->lattice, *(ControllerState*)this);
00390 }
00391 break;
00392 case SCRIPT_ATOMSENDRECV:
00393 case SCRIPT_ATOMSEND:
00394 case SCRIPT_ATOMRECV:
00395 break;
00396 case SCRIPT_MINIMIZE:
00397 minimize();
00398 break;
00399 case SCRIPT_RUN:
00400 case SCRIPT_CONTINUE:
00401 integrate(scriptTask);
00402 break;
00403 default:
00404 NAMD_bug("Unknown task in Controller::algorithm");
00405 }
00406 BackEnd::awaken();
00407 }
00408 enqueueCollections(END_OF_RUN);
00409 outputExtendedSystem(END_OF_RUN);
00410 terminate();
00411 }
00412
00413
00414
00415
00416
00417
00418 extern int eventEndOfTimeStep;
00419
00420
00421 static int gotsigint = 0;
00422 static void my_sigint_handler(int sig) {
00423 if (sig == SIGINT) gotsigint = 1;
00424 }
00425 extern "C" {
00426 typedef void (*namd_sighandler_t)(int);
00427 }
00428
00429 void Controller::integrate(int scriptTask) {
00430 char traceNote[24];
00431
00432 int step = simParams->firstTimestep;
00433
00434 const int numberOfSteps = simParams->N;
00435 const int stepsPerCycle = simParams->stepsPerCycle;
00436
00437 const int zeroMomentum = simParams->zeroMomentum;
00438
00439 nbondFreq = simParams->nonbondedFrequency;
00440 const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00441 if (dofull)
00442 slowFreq = simParams->fullElectFrequency;
00443 else
00444 slowFreq = simParams->nonbondedFrequency;
00445 if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
00446
00447 if ( scriptTask == SCRIPT_RUN ) {
00448
00449 reassignVelocities(step);
00450 rescaleaccelMD(step);
00451 adaptTempUpdate(step);
00452
00453 receivePressure(step);
00454 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00455 correctMomentum(step);
00456 printFepMessage(step);
00457 printTiMessage(step);
00458 printDynamicsEnergies(step);
00459 outputFepEnergy(step);
00460 outputTiEnergy(step);
00461 if(traceIsOn()){
00462 traceUserEvent(eventEndOfTimeStep);
00463 sprintf(traceNote, "s:%d", step);
00464 traceUserSuppliedNote(traceNote);
00465 }
00466 outputExtendedSystem(step);
00467 rebalanceLoad(step);
00468
00469 }
00470
00471
00472
00473
00474
00475
00476 for ( ++step ; step <= numberOfSteps; ++step )
00477 {
00478
00479 adaptTempUpdate(step);
00480 rescaleVelocities(step);
00481 tcoupleVelocities(step);
00482 stochRescaleVelocities(step);
00483 berendsenPressure(step);
00484 langevinPiston1(step);
00485 rescaleaccelMD(step);
00486 enqueueCollections(step);
00487 receivePressure(step);
00488 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00489 correctMomentum(step);
00490 langevinPiston2(step);
00491 reassignVelocities(step);
00492
00493 multigratorTemperature(step, 1);
00494 multigratorPressure(step, 1);
00495 multigratorPressure(step, 2);
00496 multigratorTemperature(step, 2);
00497
00498 printDynamicsEnergies(step);
00499 outputFepEnergy(step);
00500 outputTiEnergy(step);
00501 if(traceIsOn()){
00502 traceUserEvent(eventEndOfTimeStep);
00503 sprintf(traceNote, "s:%d", step);
00504 traceUserSuppliedNote(traceNote);
00505 }
00506
00507
00508
00509
00510 outputExtendedSystem(step);
00511 #if CYCLE_BARRIER
00512 cycleBarrier(!((step+1) % stepsPerCycle),step);
00513 #elif PME_BARRIER
00514 cycleBarrier(dofull && !(step%slowFreq),step);
00515 #elif STEP_BARRIER
00516 cycleBarrier(1, step);
00517 #endif
00518
00519 if(Node::Object()->specialTracing || simParams->statsOn){
00520 int bstep = simParams->traceStartStep;
00521 int estep = bstep + simParams->numTraceSteps;
00522 if(step == bstep){
00523 traceBarrier(1, step);
00524 }
00525 if(step == estep){
00526 traceBarrier(0, step);
00527 }
00528 }
00529
00530 #ifdef MEASURE_NAMD_WITH_PAPI
00531 if(simParams->papiMeasure) {
00532 int bstep = simParams->papiMeasureStartStep;
00533 int estep = bstep + simParams->numPapiMeasureSteps;
00534 if(step == bstep) {
00535 papiMeasureBarrier(1, step);
00536 }
00537 if(step == estep) {
00538 papiMeasureBarrier(0, step);
00539 }
00540 }
00541 #endif
00542
00543 rebalanceLoad(step);
00544
00545 #if PME_BARRIER
00546 cycleBarrier(dofull && !((step+1)%slowFreq),step);
00547 #endif
00548 }
00549
00550 }
00551
00552
00553 #define CALCULATE \
00554 printMinimizeEnergies(step); \
00555 outputExtendedSystem(step); \
00556 rebalanceLoad(step); \
00557 if ( step == numberOfSteps ) return; \
00558 else ++step;
00559
00560 #define MOVETO(X) \
00561 if ( step == numberOfSteps ) { \
00562 if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \
00563 if ( newDir || (mid.x-last.x) ) { \
00564 broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \
00565 } else { \
00566 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00567 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00568 min_reduction->require(); \
00569 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00570 } \
00571 enqueueCollections(step); \
00572 CALCULATE \
00573 } else if ( (X)-last.x ) { \
00574 broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \
00575 newDir = 0; \
00576 last.x = (X); \
00577 enqueueCollections(step); \
00578 CALCULATE \
00579 last.u = min_energy; \
00580 last.dudx = -1. * min_f_dot_v; \
00581 last.noGradient = min_huge_count; \
00582 if ( minVerbose ) { \
00583 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \
00584 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \
00585 iout << "\n" << endi; \
00586 } \
00587 }
00588
00589 struct minpoint {
00590 BigReal x, u, dudx, noGradient;
00591 minpoint() : x(0), u(0), dudx(0), noGradient(1) { ; }
00592 };
00593
00594 void Controller::minimize() {
00595
00596
00597 const int minVerbose = simParams->minVerbose;
00598 const int numberOfSteps = simParams->N;
00599 int step = simParams->firstTimestep;
00600 slowFreq = nbondFreq = 1;
00601 BigReal linegoal = simParams->minLineGoal;
00602 const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
00603
00604 CALCULATE
00605
00606 int minSeq = 0;
00607
00608
00609 int old_min_huge_count = 2000000000;
00610 int steps_at_same_min_huge_count = 0;
00611 for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
00612 if ( min_huge_count >= old_min_huge_count ) {
00613 if ( ++steps_at_same_min_huge_count > 10 ) break;
00614 } else {
00615 old_min_huge_count = min_huge_count;
00616 steps_at_same_min_huge_count = 0;
00617 }
00618 iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
00619 broadcast->minimizeCoefficient.publish(minSeq++,1.);
00620 enqueueCollections(step);
00621 CALCULATE
00622 }
00623 if ( min_huge_count ) {
00624 iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
00625 }
00626 iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00627
00628 int atStart = 2;
00629 int errorFactor = 10;
00630 BigReal old_f_dot_f = min_f_dot_f;
00631 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00632 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00633 int newDir = 1;
00634 min_f_dot_v = min_f_dot_f;
00635 min_v_dot_v = min_f_dot_f;
00636 while ( 1 ) {
00637
00638
00639 minpoint lo,hi,mid,last;
00640 BigReal x = 0;
00641 lo.x = x;
00642 lo.u = min_energy;
00643 lo.dudx = -1. * min_f_dot_v;
00644 lo.noGradient = min_huge_count;
00645 mid = lo;
00646 last = mid;
00647 if ( minVerbose ) {
00648 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
00649 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
00650 iout << "\n" << endi;
00651 }
00652 BigReal tol = fabs( linegoal * min_f_dot_v );
00653 iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
00654 fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
00655 int start_with_huge = last.noGradient;
00656 min_reduction->require();
00657 BigReal maxstep = 0.1 / sqrt(min_reduction->item(0));
00658 x = maxstep; MOVETO(x);
00659
00660 while ( last.u < mid.u ) {
00661 if ( last.dudx < mid.dudx * (5.5 - x/maxstep)/5. ) {
00662 x = 6 * maxstep; break;
00663 }
00664 lo = mid; mid = last;
00665 x += maxstep;
00666 if ( x > 5.5 * maxstep ) break;
00667 MOVETO(x)
00668 }
00669 if ( x > 5.5 * maxstep ) {
00670 iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" << endi;
00671 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00672 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00673 newDir = 1;
00674 old_f_dot_f = min_f_dot_f;
00675 min_f_dot_v = min_f_dot_f;
00676 min_v_dot_v = min_f_dot_f;
00677 continue;
00678 }
00679 hi = last;
00680 #define PRINT_BRACKET \
00681 iout << "LINE MINIMIZER BRACKET: DX " \
00682 << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
00683 " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
00684 lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
00685 PRINT_BRACKET
00686
00687 int itcnt;
00688 for ( itcnt = 10; itcnt > 0; --itcnt ) {
00689 int progress = 1;
00690
00691 if ( mid.noGradient ) {
00692 if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) {
00693 x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
00694 MOVETO(x)
00695 if ( last.u <= mid.u ) {
00696 hi = mid; mid = last;
00697 } else {
00698 lo = last;
00699 }
00700 } else {
00701 x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
00702 MOVETO(x)
00703 if ( last.u <= mid.u ) {
00704 lo = mid; mid = last;
00705 } else {
00706 hi = last;
00707 }
00708 }
00709 } else {
00710 if ( mid.dudx > 0. ) {
00711 BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
00712 BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
00713 x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
00714 BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
00715 if ( xdiv ) x /= xdiv; else x = last.x;
00716 if ( x > altxhi ) x = altxhi;
00717 if ( x < altxlo ) x = altxlo;
00718 if ( x-last.x == 0 ) break;
00719 MOVETO(x)
00720 if ( last.u <= mid.u ) {
00721 hi = mid; mid = last;
00722 } else {
00723 if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
00724 lo = last;
00725 }
00726 } else {
00727 BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
00728 BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
00729 x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
00730 BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
00731 if ( xdiv ) x /= xdiv; else x = last.x;
00732 if ( x < altxlo ) x = altxlo;
00733 if ( x > altxhi ) x = altxhi;
00734 if ( x-last.x == 0 ) break;
00735 MOVETO(x)
00736 if ( last.u <= mid.u ) {
00737 lo = mid; mid = last;
00738 } else {
00739 if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
00740 hi = last;
00741 }
00742 }
00743 }
00744 PRINT_BRACKET
00745 if ( ! progress ) {
00746 MOVETO(mid.x);
00747 break;
00748 }
00749 if ( fabs(last.dudx) < tol ) break;
00750 if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
00751 if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
00752 }
00753
00754 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00755 BigReal c = min_f_dot_f / old_f_dot_f;
00756 c = ( c > 1.5 ? 1.5 : c );
00757 if ( atStart ) { c = 0; --atStart; }
00758 if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
00759 c = 0;
00760 if ( errorFactor < 100 ) errorFactor += 10;
00761 }
00762 if ( c == 0 ) {
00763 iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00764 }
00765 broadcast->minimizeCoefficient.publish(minSeq++,c);
00766 newDir = 1;
00767 old_f_dot_f = min_f_dot_f;
00768 min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
00769 min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
00770 }
00771
00772 }
00773
00774 #undef MOVETO
00775 #undef CALCULATE
00776
00777
00778 void Controller::multigratorPressure(int step, int callNumber) {
00779 if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
00780 BigReal P0 = simParams->multigratorPressureTarget;
00781 BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00782 BigReal NG = controlNumDegFreedom;
00783 BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
00784 BigReal tau = simParams->multigratorPressureRelaxationTime;
00785 BigReal s = 0.5*simParams->multigratorPressureFreq*simParams->dt/tau;
00786 {
00787
00788 BigReal V = state->lattice.volume();
00789 BigReal Pinst = trace(controlPressure)/3.0;
00790 BigReal PGsum = trace(momentumSqrSum);
00791
00792 multigratorXiT = multigratorXi + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
00793 BigReal scale = exp(s*NGfac*multigratorXiT);
00794 BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
00795
00796 Tensor scaleTensor = Tensor::identity(scale);
00797 Tensor volScaleTensor = Tensor::identity(scale);
00798 Tensor velScaleTensor = Tensor::identity(velScale);
00799 state->lattice.rescale(volScaleTensor);
00800 if (callNumber == 1) {
00801 broadcast->positionRescaleFactor.publish(step,scaleTensor);
00802 broadcast->velocityRescaleTensor.publish(step,velScaleTensor);
00803 } else {
00804 broadcast->positionRescaleFactor2.publish(step,scaleTensor);
00805 broadcast->velocityRescaleTensor2.publish(step,velScaleTensor);
00806 }
00807 }
00808
00809 {
00810
00811 reduction->require();
00812 Tensor virial_normal;
00813 Tensor virial_nbond;
00814 Tensor virial_slow;
00815 Tensor intVirial_normal;
00816 Tensor intVirial_nbond;
00817 Tensor intVirial_slow;
00818 Vector extForce_normal;
00819 Vector extForce_nbond;
00820 Vector extForce_slow;
00821 GET_TENSOR(momentumSqrSum, reduction, REDUCTION_MOMENTUM_SQUARED);
00822 GET_TENSOR(virial_normal, reduction, REDUCTION_VIRIAL_NORMAL);
00823 GET_TENSOR(virial_nbond, reduction, REDUCTION_VIRIAL_NBOND);
00824 GET_TENSOR(virial_slow, reduction, REDUCTION_VIRIAL_SLOW);
00825 GET_TENSOR(intVirial_normal, reduction, REDUCTION_INT_VIRIAL_NORMAL);
00826 GET_TENSOR(intVirial_nbond, reduction, REDUCTION_INT_VIRIAL_NBOND);
00827 GET_TENSOR(intVirial_slow, reduction, REDUCTION_INT_VIRIAL_SLOW);
00828 GET_VECTOR(extForce_normal, reduction, REDUCTION_EXT_FORCE_NORMAL);
00829 GET_VECTOR(extForce_nbond, reduction, REDUCTION_EXT_FORCE_NBOND);
00830 GET_VECTOR(extForce_slow, reduction, REDUCTION_EXT_FORCE_SLOW);
00831 calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
00832 intVirial_normal, intVirial_nbond, intVirial_slow,
00833 extForce_normal, extForce_nbond, extForce_slow);
00834 if (callNumber == 2) {
00835
00836 BigReal kineticEnergy = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
00837 temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
00838 }
00839 }
00840
00841 {
00842
00843 BigReal V = state->lattice.volume();
00844 BigReal Pinst = trace(controlPressure)/3.0;
00845 BigReal PGsum = trace(momentumSqrSum);
00846
00847 multigratorXi = multigratorXiT + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
00848 }
00849
00850 }
00851 }
00852
00853 void Controller::multigratorTemperature(int step, int callNumber) {
00854 if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
00855 BigReal tau = simParams->multigratorTemperatureRelaxationTime;
00856 BigReal t = 0.5*simParams->multigratorTemperatureFreq*simParams->dt;
00857 BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00858 BigReal Nf = numDegFreedom;
00859 int n = simParams->multigratorNoseHooverChainLength;
00860 BigReal T1, T2, v;
00861 {
00862 T1 = temperature;
00863 BigReal kTinst = BOLTZMANN * temperature;
00864 for (int i=n-1;i >= 0;i--) {
00865 if (i == 0) {
00866 BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
00867 multigratorNuT[0] = exp(-0.5*t*NuOmega)*multigratorNu[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
00868 } else if (i == n-1) {
00869 multigratorNuT[n-1] = multigratorNu[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
00870 } else {
00871 BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
00872 multigratorNuT[i] = exp(-0.5*t*NuOmega)*multigratorNu[i] +
00873 0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
00874 }
00875 }
00876 BigReal velScale = exp(-t*multigratorNuT[0]/multigratorOmega[0]);
00877 v = velScale;
00878 if (callNumber == 1)
00879 broadcast->velocityRescaleFactor.publish(step,velScale);
00880 else
00881 broadcast->velocityRescaleFactor2.publish(step,velScale);
00882 }
00883
00884 {
00885
00886 multigratorReduction->require();
00887 BigReal kineticEnergy = multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY);
00888 temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
00889 T2 = temperature;
00890 if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
00891
00892 GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
00893 }
00894 }
00895
00896
00897
00898 {
00899 BigReal kTinst = BOLTZMANN * temperature;
00900 for (int i=0;i < n;i++) {
00901 if (i == 0) {
00902 BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
00903 multigratorNu[0] = exp(-0.5*t*NuOmega)*multigratorNuT[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
00904 } else if (i == n-1) {
00905 multigratorNu[n-1] = multigratorNuT[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
00906 } else {
00907 BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
00908 multigratorNu[i] = exp(-0.5*t*NuOmega)*multigratorNuT[i] +
00909 0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
00910 }
00911 multigratorZeta[i] += t*multigratorNuT[i]/multigratorOmega[i];
00912 }
00913 }
00914
00915 }
00916 }
00917
00918
00919 BigReal Controller::multigatorCalcEnthalpy(BigReal potentialEnergy, int step, int minimize) {
00920 Node *node = Node::Object();
00921 Molecule *molecule = node->molecule;
00922 SimParameters *simParameters = node->simParameters;
00923
00924 BigReal V = state->lattice.volume();
00925 BigReal P0 = simParams->multigratorPressureTarget;
00926 BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00927 BigReal NG = controlNumDegFreedom;
00928 BigReal Nf = numDegFreedom;
00929 BigReal tauV = simParams->multigratorPressureRelaxationTime;
00930 BigReal sumZeta = 0.0;
00931 for (int i=1;i < simParams->multigratorNoseHooverChainLength;i++) {
00932 sumZeta += multigratorZeta[i];
00933 }
00934 BigReal nuOmegaSum = 0.0;
00935 for (int i=0;i < simParams->multigratorNoseHooverChainLength;i++) {
00936 nuOmegaSum += multigratorNu[i]*multigratorNu[i]/(2.0*multigratorOmega[i]);
00937 }
00938 BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
00939 BigReal eta = sqrt(kT0*W)*multigratorXi;
00940
00941 BigReal enthalpy = kineticEnergy + potentialEnergy + eta*eta/(2.0*W) + P0*V + nuOmegaSum + kT0*(Nf*multigratorZeta[0] + sumZeta);
00942
00943
00944
00945
00946
00947 return enthalpy;
00948 }
00949
00950 void Controller::berendsenPressure(int step)
00951 {
00952 if ( simParams->berendsenPressureOn ) {
00953 berendsenPressure_count += 1;
00954 berendsenPressure_avg += controlPressure;
00955 const int freq = simParams->berendsenPressureFreq;
00956 if ( ! (berendsenPressure_count % freq) ) {
00957 Tensor factor = berendsenPressure_avg / berendsenPressure_count;
00958 berendsenPressure_avg = 0;
00959 berendsenPressure_count = 0;
00960
00961 factor = Tensor::diagonal(diagonal(factor));
00962 factor -= Tensor::identity(simParams->berendsenPressureTarget);
00963 factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
00964 simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
00965 factor += Tensor::identity(1.0);
00966 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
00967 if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
00968 if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
00969 int limited = 0;
00970 LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
00971 LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
00972 LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
00973 #undef LIMIT_SCALING
00974 if ( limited ) {
00975 iout << iERROR << "Step " << step <<
00976 " cell rescaling factor limited.\n" << endi;
00977 }
00978 broadcast->positionRescaleFactor.publish(step,factor);
00979 state->lattice.rescale(factor);
00980 }
00981 } else {
00982 berendsenPressure_avg = 0;
00983 berendsenPressure_count = 0;
00984 }
00985 }
00986
00987 void Controller::langevinPiston1(int step)
00988 {
00989 if ( simParams->langevinPistonOn )
00990 {
00991 Tensor &strainRate = langevinPiston_strainRate;
00992 int cellDims = simParams->useFlexibleCell ? 1 : 3;
00993 BigReal dt = simParams->dt;
00994 BigReal dt_long = slowFreq * dt;
00995 BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
00996 BigReal tau = simParams->langevinPistonPeriod;
00997 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00998
00999 #ifdef DEBUG_PRESSURE
01000 iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
01001 #endif
01002
01003 if ( ! ( (step-1) % slowFreq ) )
01004 {
01005 BigReal gamma = 1 / simParams->langevinPistonDecay;
01006 BigReal f1 = exp( -0.5 * dt_long * gamma );
01007 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
01008 strainRate *= f1;
01009 if ( simParams->useFlexibleCell ) {
01010
01011 if ( simParams->useConstantRatio ) {
01012 BigReal r = f2 * random->gaussian();
01013 strainRate.xx += r;
01014 strainRate.yy += r;
01015 strainRate.zz += f2 * random->gaussian();
01016 } else {
01017 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
01018 }
01019 } else {
01020 strainRate += f2 * Tensor::identity(random->gaussian());
01021 }
01022
01023 #ifdef DEBUG_PRESSURE
01024 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
01025 #endif
01026 }
01027
01028
01029
01030
01031 Tensor ptarget;
01032 ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
01033 if ( simParams->surfaceTensionTarget != 0. ) {
01034 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
01035 simParams->surfaceTensionTarget / state->lattice.c().z;
01036 }
01037
01038 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01039 ( controlPressure - ptarget );
01040
01041 if (simParams->fixCellDims) {
01042 if (simParams->fixCellDimX) strainRate.xx = 0;
01043 if (simParams->fixCellDimY) strainRate.yy = 0;
01044 if (simParams->fixCellDimZ) strainRate.zz = 0;
01045 }
01046
01047
01048 #ifdef DEBUG_PRESSURE
01049 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
01050 #endif
01051
01052 if ( simParams->langevinPistonBarrier ) {
01053 if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
01054 {
01055
01056 Tensor factor;
01057 if ( !simParams->useConstantArea ) {
01058 factor.xx = exp( dt_long * strainRate.xx );
01059 factor.yy = exp( dt_long * strainRate.yy );
01060 } else {
01061 factor.xx = factor.yy = 1;
01062 }
01063 factor.zz = exp( dt_long * strainRate.zz );
01064 broadcast->positionRescaleFactor.publish(step,factor);
01065 state->lattice.rescale(factor);
01066 #ifdef DEBUG_PRESSURE
01067 iout << iINFO << "rescaling by: " << factor << "\n";
01068 #endif
01069 }
01070 } else {
01071 if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
01072 {
01073 if ( ! positionRescaleFactor.xx ) {
01074
01075 Tensor factor;
01076 if ( !simParams->useConstantArea ) {
01077 factor.xx = exp( dt_long * strainRate.xx );
01078 factor.yy = exp( dt_long * strainRate.yy );
01079 } else {
01080 factor.xx = factor.yy = 1;
01081 }
01082 factor.zz = exp( dt_long * strainRate.zz );
01083 broadcast->positionRescaleFactor.publish(step,factor);
01084 positionRescaleFactor = factor;
01085 strainRate_old = strainRate;
01086 }
01087 state->lattice.rescale(positionRescaleFactor);
01088 #ifdef DEBUG_PRESSURE
01089 iout << iINFO << "rescaling by: " << positionRescaleFactor << "\n";
01090 #endif
01091 }
01092 if ( ! ( (step-slowFreq/2) % slowFreq ) )
01093 {
01094
01095 Tensor factor;
01096 if ( !simParams->useConstantArea ) {
01097 factor.xx = exp( (dt+dt_long) * strainRate.xx - dt * strainRate_old.xx );
01098 factor.yy = exp( (dt+dt_long) * strainRate.yy - dt * strainRate_old.yy );
01099 } else {
01100 factor.xx = factor.yy = 1;
01101 }
01102 factor.zz = exp( (dt+dt_long) * strainRate.zz - dt * strainRate_old.zz );
01103 broadcast->positionRescaleFactor.publish(step+1,factor);
01104 positionRescaleFactor = factor;
01105 strainRate_old = strainRate;
01106 }
01107 }
01108
01109
01110 if ( ! ( step % nbondFreq ) )
01111 {
01112 #ifdef DEBUG_PRESSURE
01113 iout << iINFO << "correcting strain rate for nbond, ";
01114 #endif
01115 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01116 ( (nbondFreq - 1) * controlPressure_nbond );
01117 #ifdef DEBUG_PRESSURE
01118 iout << "strain rate: " << strainRate << "\n";
01119 #endif
01120 }
01121 if ( ! ( step % slowFreq ) )
01122 {
01123 #ifdef DEBUG_PRESSURE
01124 iout << iINFO << "correcting strain rate for slow, ";
01125 #endif
01126 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01127 ( (slowFreq - 1) * controlPressure_slow );
01128 #ifdef DEBUG_PRESSURE
01129 iout << "strain rate: " << strainRate << "\n";
01130 #endif
01131 }
01132 if (simParams->fixCellDims) {
01133 if (simParams->fixCellDimX) strainRate.xx = 0;
01134 if (simParams->fixCellDimY) strainRate.yy = 0;
01135 if (simParams->fixCellDimZ) strainRate.zz = 0;
01136 }
01137
01138 }
01139
01140 }
01141
01142 void Controller::langevinPiston2(int step)
01143 {
01144 if ( simParams->langevinPistonOn )
01145 {
01146 Tensor &strainRate = langevinPiston_strainRate;
01147 int cellDims = simParams->useFlexibleCell ? 1 : 3;
01148 BigReal dt = simParams->dt;
01149 BigReal dt_long = slowFreq * dt;
01150 BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
01151 BigReal tau = simParams->langevinPistonPeriod;
01152 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
01153
01154
01155 if ( ! ( step % nbondFreq ) )
01156 {
01157 #ifdef DEBUG_PRESSURE
01158 iout << iINFO << "correcting strain rate for nbond, ";
01159 #endif
01160 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01161 ( (nbondFreq - 1) * controlPressure_nbond );
01162 #ifdef DEBUG_PRESSURE
01163 iout << "strain rate: " << strainRate << "\n";
01164 #endif
01165 }
01166 if ( ! ( step % slowFreq ) )
01167 {
01168 #ifdef DEBUG_PRESSURE
01169 iout << iINFO << "correcting strain rate for slow, ";
01170 #endif
01171 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01172 ( (slowFreq - 1) * controlPressure_slow );
01173 #ifdef DEBUG_PRESSURE
01174 iout << "strain rate: " << strainRate << "\n";
01175 #endif
01176 }
01177
01178
01179
01180
01181 Tensor ptarget;
01182 ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
01183 if ( simParams->surfaceTensionTarget != 0. ) {
01184 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
01185 simParams->surfaceTensionTarget / state->lattice.c().z;
01186 }
01187
01188 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01189 ( controlPressure - ptarget );
01190
01191
01192 #ifdef DEBUG_PRESSURE
01193 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
01194 #endif
01195
01196 if ( ! ( step % slowFreq ) )
01197 {
01198 BigReal gamma = 1 / simParams->langevinPistonDecay;
01199 BigReal f1 = exp( -0.5 * dt_long * gamma );
01200 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
01201 strainRate *= f1;
01202 if ( simParams->useFlexibleCell ) {
01203
01204 if ( simParams->useConstantRatio ) {
01205 BigReal r = f2 * random->gaussian();
01206 strainRate.xx += r;
01207 strainRate.yy += r;
01208 strainRate.zz += f2 * random->gaussian();
01209 } else {
01210 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
01211 }
01212 } else {
01213 strainRate += f2 * Tensor::identity(random->gaussian());
01214 }
01215 #ifdef DEBUG_PRESSURE
01216 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
01217 #endif
01218 }
01219
01220 #ifdef DEBUG_PRESSURE
01221 iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
01222 #endif
01223 if (simParams->fixCellDims) {
01224 if (simParams->fixCellDimX) strainRate.xx = 0;
01225 if (simParams->fixCellDimY) strainRate.yy = 0;
01226 if (simParams->fixCellDimZ) strainRate.zz = 0;
01227 }
01228 }
01229
01230
01231 }
01232
01233 void Controller::rescaleVelocities(int step)
01234 {
01235 const int rescaleFreq = simParams->rescaleFreq;
01236 if ( rescaleFreq > 0 ) {
01237 rescaleVelocities_sumTemps += temperature; ++rescaleVelocities_numTemps;
01238 if ( rescaleVelocities_numTemps == rescaleFreq ) {
01239 BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
01240 BigReal rescaleTemp = simParams->rescaleTemp;
01241 if ( simParams->adaptTempOn && simParams->adaptTempRescale && (step > simParams->adaptTempFirstStep ) &&
01242 (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
01243 rescaleTemp = adaptTempT;
01244 }
01245 BigReal factor = sqrt(rescaleTemp/avgTemp);
01246 broadcast->velocityRescaleFactor.publish(step,factor);
01247
01248
01249
01250 rescaleVelocities_sumTemps = 0; rescaleVelocities_numTemps = 0;
01251 }
01252 }
01253 }
01254
01255 void Controller::correctMomentum(int step) {
01256
01257 Vector momentum;
01258 momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
01259 momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
01260 momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
01261 const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
01262
01263 if ( momentum.length2() == 0. )
01264 iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
01265 if ( mass == 0. )
01266 NAMD_die("Total mass is zero in Controller::correctMomentum");
01267
01268 momentum *= (-1./mass);
01269
01270 broadcast->momentumCorrection.publish(step+slowFreq,momentum);
01271 }
01272
01273
01274 void Controller::printFepMessage(int step)
01275 {
01276 if (simParams->alchOn && simParams->alchFepOn
01277 && !simParams->alchLambdaFreq) {
01278 const BigReal alchLambda = simParams->alchLambda;
01279 const BigReal alchLambda2 = simParams->alchLambda2;
01280 const BigReal alchLambdaIDWS = simParams->alchLambdaIDWS;
01281 const BigReal alchTemp = simParams->alchTemp;
01282 const int alchEquilSteps = simParams->alchEquilSteps;
01283 iout << "FEP: RESETTING FOR NEW FEP WINDOW "
01284 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2;
01285 if ( alchLambdaIDWS >= 0. ) {
01286 iout << " LAMBDA_IDWS " << alchLambdaIDWS;
01287 }
01288 iout << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
01289 << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
01290 << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
01291 << " K FOR FEP CALCULATION\n" << endi;
01292 }
01293 }
01294 void Controller::printTiMessage(int step)
01295 {
01296 if (simParams->alchOn && simParams->alchThermIntOn
01297 && !simParams->alchLambdaFreq) {
01298 const BigReal alchLambda = simParams->alchLambda;
01299 const int alchEquilSteps = simParams->alchEquilSteps;
01300 iout << "TI: RESETTING FOR NEW WINDOW "
01301 << "LAMBDA SET TO " << alchLambda
01302 << "\nTI: WINDOW TO HAVE " << alchEquilSteps
01303 << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
01304 }
01305 }
01306
01307
01308 void Controller::reassignVelocities(int step)
01309 {
01310 const int reassignFreq = simParams->reassignFreq;
01311 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
01312 BigReal newTemp = simParams->reassignTemp;
01313 newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
01314 if ( simParams->reassignIncr > 0.0 ) {
01315 if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
01316 newTemp = simParams->reassignHold;
01317 } else {
01318 if ( newTemp < simParams->reassignHold )
01319 newTemp = simParams->reassignHold;
01320 }
01321 iout << "REASSIGNING VELOCITIES AT STEP " << step
01322 << " TO " << newTemp << " KELVIN.\n" << endi;
01323 }
01324 }
01325
01326 void Controller::tcoupleVelocities(int step)
01327 {
01328 if ( simParams->tCoupleOn )
01329 {
01330 const BigReal tCoupleTemp = simParams->tCoupleTemp;
01331 BigReal coefficient = 1.;
01332 if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
01333 broadcast->tcoupleCoefficient.publish(step,coefficient);
01334 }
01335 }
01336
01348 void Controller::stochRescaleVelocities(int step)
01349 {
01350 if ( simParams->stochRescaleOn ) {
01351 ++stochRescale_count;
01352 if ( stochRescale_count == simParams->stochRescaleFreq ) {
01353 double coefficient = stochRescaleCoefficient();
01354 broadcast->stochRescaleCoefficient.publish(step,coefficient);
01355 stochRescale_count = 0;
01356 }
01357 }
01358 }
01359
01364 double Controller::stochRescaleCoefficient() {
01365 const double stochRescaleTemp = simParams->stochRescaleTemp;
01366 double coefficient = 1;
01367 if ( temperature > 0 ) {
01368 double R1 = random->gaussian();
01369
01370
01371
01372
01373
01374
01375
01376 double R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
01377 double tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
01378
01379 coefficient = sqrt(stochRescaleTimefactor +
01380 (1 - stochRescaleTimefactor)*tempfactor*R2sum +
01381 2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*
01382 stochRescaleTimefactor));
01383 }
01384 heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient-1);
01385 return coefficient;
01386 }
01387
01388 static char *FORMAT(BigReal X)
01389 {
01390 static char tmp_string[25];
01391 const double maxnum = 99999999999.9999;
01392 if ( X > maxnum ) X = maxnum;
01393 if ( X < -maxnum ) X = -maxnum;
01394 sprintf(tmp_string," %14.4f",X);
01395 return tmp_string;
01396 }
01397
01398 static char *FORMAT(const char *X)
01399 {
01400 static char tmp_string[25];
01401 sprintf(tmp_string," %14s",X);
01402 return tmp_string;
01403 }
01404
01405 static char *ETITLE(int X)
01406 {
01407 static char tmp_string[21];
01408 sprintf(tmp_string,"ENERGY: %7d",X);
01409 return tmp_string;
01410 }
01411
01412 void Controller::receivePressure(int step, int minimize)
01413 {
01414 Node *node = Node::Object();
01415 Molecule *molecule = node->molecule;
01416 SimParameters *simParameters = node->simParameters;
01417 Lattice &lattice = state->lattice;
01418
01419 reduction->require();
01420
01421 Tensor virial_normal;
01422 Tensor virial_nbond;
01423 Tensor virial_slow;
01424 #ifdef ALTVIRIAL
01425 Tensor altVirial_normal;
01426 Tensor altVirial_nbond;
01427 Tensor altVirial_slow;
01428 #endif
01429 Tensor intVirial_normal;
01430 Tensor intVirial_nbond;
01431 Tensor intVirial_slow;
01432 Vector extForce_normal;
01433 Vector extForce_nbond;
01434 Vector extForce_slow;
01435
01436 #if 1
01437 numDegFreedom = molecule->num_deg_freedom();
01438 int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
01439 int numFixedGroups = molecule->num_fixed_groups();
01440 int numFixedAtoms = molecule->num_fixed_atoms();
01441 #endif
01442 #if 0
01443 int numAtoms = molecule->numAtoms;
01444 numDegFreedom = 3 * numAtoms;
01445 int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
01446 int numFixedAtoms =
01447 ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
01448 int numLonepairs = molecule->numLonepairs;
01449 int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
01450 if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
01451 if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01452 if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01453 if ( ! ( numFixedAtoms || molecule->numConstraints
01454 || simParameters->comMove || simParameters->langevinOn ) ) {
01455 numDegFreedom -= 3;
01456 numGroupDegFreedom -= 3;
01457 }
01458 if (simParameters->pairInteractionOn) {
01459
01460 numDegFreedom = 3 * molecule->numFepInitial;
01461 }
01462 int numRigidBonds = molecule->numRigidBonds;
01463 int numFixedRigidBonds =
01464 ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01465
01466
01467
01468 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01469 #endif
01470
01471 kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01472 kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01473
01474 BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01475 reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01476 BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01477 reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01478
01479 BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01480 / ( numDegFreedom * BOLTZMANN );
01481 BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01482 / ( numDegFreedom * BOLTZMANN );
01483 BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01484 / ( numGroupDegFreedom * BOLTZMANN );
01485 BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01486 / ( numGroupDegFreedom * BOLTZMANN );
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498 GET_TENSOR(momentumSqrSum,reduction,REDUCTION_MOMENTUM_SQUARED);
01499
01500 GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01501 GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01502 GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01503
01504 #ifdef ALTVIRIAL
01505 GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01506 GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01507 GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01508 #endif
01509
01510 GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01511 GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01512 GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01513
01514 GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01515 GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01516 GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01517
01518
01519
01520
01521
01522
01523 kineticEnergy = kineticEnergyCentered;
01524 temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );
01525
01526 if (simParameters->drudeOn) {
01527 BigReal drudeComKE
01528 = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01529 BigReal drudeBondKE
01530 = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01531 int g_bond = 3 * molecule->numDrudeAtoms;
01532 int g_com = numDegFreedom - g_bond;
01533 temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
01534 drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
01535 }
01536
01537
01538 if ( simParameters->useGroupPressure )
01539 {
01540 controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01541 if ( ! ( numFixedAtoms || molecule->numConstraints
01542 || simParameters->comMove || simParameters->langevinOn ) ) {
01543 controlNumDegFreedom -= 1;
01544 }
01545 }
01546 else
01547 {
01548 controlNumDegFreedom = numDegFreedom / 3;
01549 }
01550 if (simParameters->fixCellDims) {
01551 if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01552 if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01553 if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01554 }
01555
01556
01557 calcPressure(step, minimize,
01558 virial_normal, virial_nbond, virial_slow,
01559 intVirial_normal, intVirial_nbond, intVirial_slow,
01560 extForce_normal, extForce_nbond, extForce_slow);
01561
01562 #ifdef DEBUG_PRESSURE
01563 iout << iINFO << "Control pressure = " << controlPressure <<
01564 " = " << controlPressure_normal << " + " <<
01565 controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01566 #endif
01567 if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
01568 iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
01569 << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
01570 << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
01571 << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
01572 << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
01573 << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
01574 << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
01575 << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
01576 << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
01577 << endi;
01578 }
01579 }
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590 void Controller::calcPressure(int step, int minimize,
01591 const Tensor& virial_normal_in, const Tensor& virial_nbond_in, const Tensor& virial_slow_in,
01592 const Tensor& intVirial_normal, const Tensor& intVirial_nbond, const Tensor& intVirial_slow,
01593 const Vector& extForce_normal, const Vector& extForce_nbond, const Vector& extForce_slow) {
01594
01595 Tensor virial_normal = virial_normal_in;
01596 Tensor virial_nbond = virial_nbond_in;
01597 Tensor virial_slow = virial_slow_in;
01598
01599
01600
01601
01602
01603 Node *node = Node::Object();
01604 Molecule *molecule = node->molecule;
01605 SimParameters *simParameters = node->simParameters;
01606 Lattice &lattice = state->lattice;
01607
01608 BigReal volume;
01609
01610 Vector extPosition = lattice.origin();
01611 virial_normal -= outer(extForce_normal,extPosition);
01612 virial_nbond -= outer(extForce_nbond,extPosition);
01613 virial_slow -= outer(extForce_slow,extPosition);
01614
01615 if ( (volume=lattice.volume()) != 0. )
01616 {
01617
01618 if (simParameters->LJcorrection && volume) {
01619 #ifdef MEM_OPT_VERSION
01620 NAMD_bug("LJcorrection not supported in memory optimized build.");
01621 #else
01622
01623 BigReal alchLambda = simParameters->getCurrentLambda(step);
01624 virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
01625 #endif
01626 }
01627
01628
01629 pressure_normal = virial_normal / volume;
01630 groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01631
01632 if (simParameters->accelMDOn) {
01633 pressure_amd = virial_amd / volume;
01634 pressure_normal += pressure_amd;
01635 groupPressure_normal += pressure_amd;
01636 }
01637
01638 if ( minimize || ! ( step % nbondFreq ) )
01639 {
01640 pressure_nbond = virial_nbond / volume;
01641 groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01642 }
01643
01644 if ( minimize || ! ( step % slowFreq ) )
01645 {
01646 pressure_slow = virial_slow / volume;
01647 groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01648 }
01649
01650 pressure = pressure_normal + pressure_nbond + pressure_slow;
01651 groupPressure = groupPressure_normal + groupPressure_nbond +
01652 groupPressure_slow;
01653 }
01654 else
01655 {
01656 pressure = Tensor();
01657 groupPressure = Tensor();
01658 }
01659
01660 if ( simParameters->useGroupPressure )
01661 {
01662 controlPressure_normal = groupPressure_normal;
01663 controlPressure_nbond = groupPressure_nbond;
01664 controlPressure_slow = groupPressure_slow;
01665 controlPressure = groupPressure;
01666 }
01667 else
01668 {
01669 controlPressure_normal = pressure_normal;
01670 controlPressure_nbond = pressure_nbond;
01671 controlPressure_slow = pressure_slow;
01672 controlPressure = pressure;
01673 }
01674
01675 if ( simParameters->useFlexibleCell ) {
01676
01677
01678
01679
01680
01681
01682 controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01683 controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01684 controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01685 controlPressure = Tensor::diagonal(diagonal(controlPressure));
01686 if ( simParameters->useConstantRatio ) {
01687 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01688 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01689 AVGXY(controlPressure_normal);
01690 AVGXY(controlPressure_nbond);
01691 AVGXY(controlPressure_slow);
01692 AVGXY(controlPressure);
01693 #undef AVGXY
01694 }
01695 } else {
01696 controlPressure_normal =
01697 Tensor::identity(trace(controlPressure_normal)/3.);
01698 controlPressure_nbond =
01699 Tensor::identity(trace(controlPressure_nbond)/3.);
01700 controlPressure_slow =
01701 Tensor::identity(trace(controlPressure_slow)/3.);
01702 controlPressure =
01703 Tensor::identity(trace(controlPressure)/3.);
01704 }
01705 }
01706
01707 inline void Controller :: calc_accelMDG_mean_std
01708 (BigReal testV, int n,
01709 BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV){
01710 BigReal delta;
01711
01712 if(testV > *Vmax){
01713 *Vmax = testV;
01714 }
01715 else if(testV < *Vmin){
01716 *Vmin = testV;
01717 }
01718
01719
01720 delta = testV - *Vavg;
01721 *Vavg += delta / (BigReal)n;
01722 *M2 += delta * (testV - (*Vavg));
01723
01724 *sigmaV = sqrt(*M2/n);
01725 }
01726
01727 inline void Controller :: calc_accelMDG_E_k
01728 (int iE, int step, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV,
01729 BigReal* k0, BigReal* k, BigReal* E, int *iEused, char *warn){
01730 switch(iE){
01731 case 2:
01732 *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
01733
01734 if(*k0 > 1.)
01735 *warn |= 1;
01736 else if(*k0 <= 0.)
01737 *warn |= 2;
01738
01739 else{
01740 *E = Vmin + (Vmax-Vmin)/(*k0);
01741 *iEused = 2;
01742 break;
01743 }
01744 case 1:
01745 *E = Vmax;
01746 *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
01747 if(*k0 > 1.)
01748 *k0 = 1.;
01749 *iEused = 1;
01750 break;
01751 }
01752
01753 *k = *k0 / (Vmax - Vmin);
01754 }
01755
01756 inline void Controller :: calc_accelMDG_force_factor
01757 (BigReal k, BigReal E, BigReal testV, Tensor vir_orig,
01758 BigReal *dV, BigReal *factor, Tensor *vir){
01759 BigReal VE = testV - E;
01760
01761 if(VE < 0){
01762 *factor = k * VE;
01763 *vir += vir_orig * (*factor);
01764 *dV += (*factor) * VE/2.;
01765 *factor += 1.;
01766 }
01767 else{
01768 *factor = 1.;
01769 }
01770 }
01771
01772 void Controller :: write_accelMDG_rest_file
01773 (int step_n, char type, int V_n, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal M2, BigReal E, BigReal k, bool write_topic, bool lasttime){
01774 FILE *rest;
01775 char timestepstr[20];
01776 char *filename;
01777 int baselen;
01778 char *restart_name;
01779 const char *bsuffix;
01780
01781 if(lasttime || simParams->restartFrequency == 0){
01782 filename = simParams->outputFilename;
01783 bsuffix = ".BAK";
01784 }
01785 else{
01786 filename = simParams->restartFilename;
01787 bsuffix = ".old";
01788 }
01789
01790 baselen = strlen(filename);
01791 restart_name = new char[baselen+26];
01792
01793 strcpy(restart_name, filename);
01794 if ( simParams->restartSave && !lasttime) {
01795 sprintf(timestepstr,".%d",step_n);
01796 strcat(restart_name, timestepstr);
01797 bsuffix = ".BAK";
01798 }
01799 strcat(restart_name, ".gamd");
01800
01801 if(write_topic){
01802 NAMD_backup_file(restart_name,bsuffix);
01803
01804 rest = fopen(restart_name, "w");
01805 if(rest){
01806 fprintf(rest, "# NAMD accelMDG restart file\n"
01807 "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
01808 "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01809 type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01810 fclose(rest);
01811 iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
01812 }
01813 else
01814 iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01815 }
01816 else{
01817 rest = fopen(restart_name, "a");
01818 if(rest){
01819 fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01820 type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01821 fclose(rest);
01822 }
01823 else
01824 iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01825 }
01826
01827 delete[] restart_name;
01828 }
01829
01830
01831 void Controller::rescaleaccelMD(int step, int minimize)
01832 {
01833 if ( !simParams->accelMDOn ) return;
01834
01835 amd_reduction->require();
01836
01837
01838 for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
01839 submit_reduction->item(i) += amd_reduction->item(i);
01840 }
01841 submit_reduction->submit();
01842
01843 if (step == simParams->firstTimestep) {
01844 accelMDdVAverage = 0;
01845 accelMDdV = 0;
01846 }
01847
01848 if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
01849
01850 Node *node = Node::Object();
01851 Molecule *molecule = node->molecule;
01852 Lattice &lattice = state->lattice;
01853
01854 const BigReal accelMDE = simParams->accelMDE;
01855 const BigReal accelMDalpha = simParams->accelMDalpha;
01856 const BigReal accelMDTE = simParams->accelMDTE;
01857 const BigReal accelMDTalpha = simParams->accelMDTalpha;
01858 const int accelMDOutFreq = simParams->accelMDOutFreq;
01859
01860
01861 static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
01862 static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
01863 static BigReal k0P, kP, EP;
01864 static BigReal k0D, kD, ED;
01865 static int V_n = 1, iEusedD, iEusedP;
01866 static char warnD, warnP;
01867 static int restfreq;
01868
01869 const int firststep = (simParams->accelMDFirstStep > simParams->firstTimestep) ? simParams->accelMDFirstStep : simParams->firstTimestep;
01870 const int ntcmdprep = (simParams->accelMDGcMDPrepSteps > 0) ? simParams->accelMDFirstStep + simParams->accelMDGcMDPrepSteps : 0;
01871 const int ntcmd = simParams->accelMDFirstStep + simParams->accelMDGcMDSteps;
01872 const int ntebprep = ntcmd + simParams->accelMDGEquiPrepSteps;
01873 const int nteb = ntcmd + simParams->accelMDGEquiSteps;
01874 const int ntave = simParams->accelMDGStatWindow;
01875
01876 BigReal volume;
01877 BigReal bondEnergy;
01878 BigReal angleEnergy;
01879 BigReal dihedralEnergy;
01880 BigReal improperEnergy;
01881 BigReal crosstermEnergy;
01882 BigReal boundaryEnergy;
01883 BigReal miscEnergy;
01884 BigReal amd_electEnergy;
01885 BigReal amd_ljEnergy;
01886 BigReal amd_ljEnergy_Corr = 0.;
01887 BigReal amd_electEnergySlow;
01888 BigReal amd_groLJEnergy;
01889 BigReal amd_groGaussEnergy;
01890 BigReal amd_goTotalEnergy;
01891 BigReal potentialEnergy;
01892 BigReal factor_dihe = 1;
01893 BigReal factor_tot = 1;
01894 BigReal testV;
01895 BigReal dV = 0.;
01896 Vector accelMDfactor;
01897 Tensor vir;
01898 Tensor vir_dihe;
01899 Tensor vir_normal;
01900 Tensor vir_nbond;
01901 Tensor vir_slow;
01902
01903 volume = lattice.volume();
01904
01905 bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
01906 angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
01907 dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
01908 improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
01909 crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
01910 boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
01911 miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);
01912
01913 GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
01914 GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
01915 GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
01916 GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
01917
01918 if ( !( step % nbondFreq ) ) {
01919 amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
01920 amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
01921 amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
01922 amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
01923 BigReal goNativeEnergy = amd_reduction->item(REDUCTION_GO_NATIVE_ENERGY);
01924 BigReal goNonnativeEnergy = amd_reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
01925 amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
01926 } else {
01927 amd_electEnergy = electEnergy;
01928 amd_ljEnergy = ljEnergy;
01929 amd_groLJEnergy = groLJEnergy;
01930 amd_groGaussEnergy = groGaussEnergy;
01931 amd_goTotalEnergy = goTotalEnergy;
01932 }
01933
01934 if( simParams->LJcorrection && volume ) {
01935
01936
01937
01938 amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
01939 }
01940
01941 if ( !( step % slowFreq ) ) {
01942 amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
01943 } else {
01944 amd_electEnergySlow = electEnergySlow;
01945 }
01946
01947 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
01948 improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
01949 crosstermEnergy + boundaryEnergy + miscEnergy +
01950 amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
01951
01952
01953 if(simParams->accelMDG){
01954
01955 if(step == firststep) {
01956 iEusedD = iEusedP = simParams->accelMDGiE;
01957 warnD = warnP = '\0';
01958
01959
01960 if(simParams->accelMDGRestart){
01961 FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
01962 char line[256];
01963 int dihe_n=0, tot_n=0;
01964 if(!rest){
01965 sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
01966 NAMD_die(line);
01967 }
01968
01969 while(fgets(line, 256, rest)){
01970 if(line[0] == 'D'){
01971 dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
01972 &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
01973 }
01974 else if(line[0] == 'T'){
01975 tot_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
01976 &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
01977 }
01978 }
01979
01980 fclose(rest);
01981
01982 V_n++;
01983
01984
01985 if(simParams->accelMDdihe || simParams->accelMDdual){
01986 if(dihe_n !=8)
01987 NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
01988 k0D = kD * (VmaxD - VminD);
01989 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
01990 << " Vmax " << VmaxD << " Vmin " << VminD
01991 << " Vavg " << VavgD << " sigmaV " << sigmaVD
01992 << " E " << ED << " k " << kD
01993 << "\n" << endi;
01994 }
01995
01996
01997 if(!simParams->accelMDdihe || simParams->accelMDdual){
01998 if(tot_n !=8)
01999 NAMD_die("Invalid total potential energy format in accelMDG restart file");
02000 k0P = kP * (VmaxP - VminP);
02001 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
02002 << " Vmax " << VmaxP << " Vmin " << VminP
02003 << " Vavg " << VavgP << " sigmaV " << sigmaVP
02004 << " E " << EP << " k " << kP
02005 << "\n" << endi;
02006 }
02007
02008 iEusedD = (ED == VmaxD) ? 1 : 2;
02009 iEusedP = (EP == VmaxP) ? 1 : 2;
02010 }
02011
02012 restfreq = simParams->restartFrequency ? simParams->restartFrequency : 500;
02013 }
02014
02015
02016 if(simParams->accelMDdihe || simParams->accelMDdual){
02017 testV = dihedralEnergy + crosstermEnergy;
02018
02019
02020 if(step > firststep && step % restfreq == 0)
02021 write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
02022 true, false);
02023
02024 if( simParams->accelMDLastStep > 0 ){
02025 if( step == simParams->accelMDLastStep )
02026 write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
02027 true, true);
02028 }
02029 else if(step == simParams->N)
02030 write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
02031 true, true);
02032
02033
02034 if(step < ntcmd){
02035
02036 if(step == firststep && !simParams->accelMDGRestart){
02037
02038 VmaxD = VminD = VavgD = testV;
02039 M2D = sigmaVD = 0.;
02040 }
02041
02042 else if(step == ntcmdprep && ntcmdprep != 0){
02043
02044 VmaxD = VminD = VavgD = testV;
02045 M2D = sigmaVD = 0.;
02046 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
02047 }
02048
02049 else if(ntave > 0 && step % ntave == 0){
02050
02051 if(testV > VmaxD) VmaxD = testV;
02052 if(testV < VminD) VminD = testV;
02053
02054 VavgD = testV;
02055 M2D = sigmaVD = 0.;
02056 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
02057 }
02058
02059 else
02060 calc_accelMDG_mean_std(testV, V_n,
02061 &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
02062
02063
02064 if(step == ntcmd - 1){
02065 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D,
02066 VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
02067 }
02068 }
02069
02070 else if(step < nteb){
02071 calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
02072 &dV, &factor_dihe, &vir);
02073
02074
02075 if(step == ntcmd && simParams->accelMDGresetVaftercmd){
02076
02077 VmaxD = VminD = VavgD = testV;
02078 M2D = sigmaVD = 0.;
02079 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
02080 }
02081
02082 else if(ntave > 0 && step % ntave == 0){
02083
02084 if(testV > VmaxD) VmaxD = testV;
02085 if(testV < VminD) VminD = testV;
02086
02087 VavgD = testV;
02088 M2D = sigmaVD = 0.;
02089 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
02090 }
02091 else
02092 calc_accelMDG_mean_std(testV, V_n,
02093 &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
02094
02095
02096 if(step >= ntebprep)
02097 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
02098 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D,
02099 VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
02100 }
02101
02102 else{
02103 calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
02104 &dV, &factor_dihe, &vir);
02105 }
02106
02107 }
02108
02109 if(!simParams->accelMDdihe || simParams->accelMDdual){
02110 Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
02111 testV = potentialEnergy;
02112 if(simParams->accelMDdual){
02113 testV -= dihedralEnergy + crosstermEnergy;
02114 vir_tot -= vir_dihe;
02115 }
02116
02117
02118 if(step > firststep && step % restfreq == 0)
02119 write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
02120 !simParams->accelMDdual, false);
02121
02122 if( simParams->accelMDLastStep > 0 ){
02123 if( step == simParams->accelMDLastStep )
02124 write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
02125 !simParams->accelMDdual, true);
02126 }
02127 else if(step == simParams->N)
02128 write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
02129 !simParams->accelMDdual, true);
02130
02131
02132 if(step < ntcmd){
02133
02134 if(step == firststep && !simParams->accelMDGRestart){
02135
02136 VmaxP = VminP = VavgP = testV;
02137 M2P = sigmaVP = 0.;
02138 }
02139
02140 else if(step == ntcmdprep && ntcmdprep != 0){
02141
02142 VmaxP = VminP = VavgP = testV;
02143 M2P = sigmaVP = 0.;
02144 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
02145 }
02146
02147 else if(ntave > 0 && step % ntave == 0){
02148
02149 if(testV > VmaxP) VmaxP = testV;
02150 if(testV < VminP) VminP = testV;
02151
02152 VavgP = testV;
02153 M2P = sigmaVP = 0.;
02154 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
02155 }
02156
02157 else
02158 calc_accelMDG_mean_std(testV, V_n,
02159 &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
02160
02161 if(step == ntcmd - 1){
02162 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P,
02163 VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
02164 }
02165 }
02166
02167 else if(step < nteb){
02168 calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
02169 &dV, &factor_tot, &vir);
02170
02171
02172 if(step == ntcmd && simParams->accelMDGresetVaftercmd){
02173
02174 VmaxP = VminP = VavgP = testV;
02175 M2P = sigmaVP = 0.;
02176 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
02177 }
02178
02179 else if(ntave > 0 && step % ntave == 0){
02180
02181 if(testV > VmaxP) VmaxP = testV;
02182 if(testV < VminP) VminP = testV;
02183
02184 VavgP = testV;
02185 M2P = sigmaVP = 0.;
02186 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
02187 }
02188 else
02189 calc_accelMDG_mean_std(testV, V_n,
02190 &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
02191
02192
02193 if(step >= ntebprep)
02194 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
02195 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P,
02196 VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
02197 }
02198
02199 else{
02200 calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
02201 &dV, &factor_tot, &vir);
02202 }
02203
02204 }
02205 accelMDdVAverage += dV;
02206
02207
02208 if((ntcmdprep > 0 && step == ntcmdprep) ||
02209 (step < nteb && ntave > 0 && step % ntave == 0) ||
02210 (simParams->accelMDGresetVaftercmd && step == ntcmd)){
02211 V_n = 1;
02212 }
02213
02214 if(step < nteb)
02215 V_n++;
02216
02217 }
02218
02219 else{
02220 if (simParams->accelMDdihe) {
02221
02222 testV = dihedralEnergy + crosstermEnergy;
02223 if ( testV < accelMDE ) {
02224 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
02225 factor_dihe *= factor_dihe;
02226 vir = vir_dihe * (factor_dihe - 1.0);
02227 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02228 accelMDdVAverage += dV;
02229 }
02230
02231 } else if (simParams->accelMDdual) {
02232
02233 testV = dihedralEnergy + crosstermEnergy;
02234 if ( testV < accelMDE ) {
02235 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
02236 factor_dihe *= factor_dihe;
02237 vir = vir_dihe * (factor_dihe - 1.0) ;
02238 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02239 }
02240
02241 testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
02242 if ( testV < accelMDTE ) {
02243 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
02244 factor_tot *= factor_tot;
02245 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
02246 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
02247 }
02248 accelMDdVAverage += dV;
02249
02250 } else {
02251
02252 testV = potentialEnergy;
02253 if ( testV < accelMDE ) {
02254 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
02255 factor_tot *= factor_tot;
02256 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
02257 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02258 accelMDdVAverage += dV;
02259 }
02260 }
02261 }
02262
02263 accelMDfactor[0]=factor_dihe;
02264 accelMDfactor[1]=factor_tot;
02265 accelMDfactor[2]=1;
02266 broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
02267 virial_amd = vir;
02268 accelMDdV = dV;
02269
02270 if ( factor_tot < 0.001 ) {
02271 iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
02272 << "\n" << endi;
02273 }
02274
02275 if ( ! (step % accelMDOutFreq) ) {
02276 if ( !(step == simParams->firstTimestep) ) {
02277 accelMDdVAverage = accelMDdVAverage/accelMDOutFreq;
02278 }
02279 iout << "ACCELERATED MD: STEP " << step
02280 << " dV " << dV
02281 << " dVAVG " << accelMDdVAverage
02282 << " BOND " << bondEnergy
02283 << " ANGLE " << angleEnergy
02284 << " DIHED " << dihedralEnergy+crosstermEnergy
02285 << " IMPRP " << improperEnergy
02286 << " ELECT " << amd_electEnergy+amd_electEnergySlow
02287 << " VDW " << amd_ljEnergy
02288 << " POTENTIAL " << potentialEnergy;
02289 if(amd_ljEnergy_Corr)
02290 iout << " LJcorr " << amd_ljEnergy_Corr;
02291 iout << "\n" << endi;
02292 if(simParams->accelMDG){
02293 if(simParams->accelMDdihe || simParams->accelMDdual)
02294 iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD
02295 << " Vmax " << VmaxD << " Vmin " << VminD
02296 << " Vavg " << VavgD << " sigmaV " << sigmaVD
02297 << " E " << ED << " k0 " << k0D << " k " << kD
02298 << "\n" << endi;
02299 if(warnD & 1)
02300 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
02301 << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02302 if(warnD & 2)
02303 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
02304 << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02305 if(!simParams->accelMDdihe || simParams->accelMDdual)
02306 iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP
02307 << " Vmax " << VmaxP << " Vmin " << VminP
02308 << " Vavg " << VavgP << " sigmaV " << sigmaVP
02309 << " E " << EP << " k0 " << k0P << " k " << kP
02310 << "\n" << endi;
02311 if(warnP & 1)
02312 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
02313 << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02314 if(warnP & 2)
02315 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
02316 << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02317 warnD = warnP = '\0';
02318 }
02319
02320 accelMDdVAverage = 0;
02321
02322 if (simParams->accelMDDebugOn) {
02323 Tensor p_normal;
02324 Tensor p_nbond;
02325 Tensor p_slow;
02326 Tensor p;
02327 if ( volume != 0. ) {
02328 p_normal = vir_normal/volume;
02329 p_nbond = vir_nbond/volume;
02330 p_slow = vir_slow/volume;
02331 p = vir/volume;
02332 }
02333 iout << " accelMD Scaling Factor: " << accelMDfactor << "\n"
02334 << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
02335 << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
02336 << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
02337 << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n"
02338 << endi;
02339 }
02340 }
02341 }
02342
02343 void Controller::adaptTempInit(int step) {
02344 if (!simParams->adaptTempOn) return;
02345 iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
02346 adaptTempDtMin = 0;
02347 adaptTempDtMax = 0;
02348 adaptTempAutoDt = false;
02349 if (simParams->adaptTempBins == 0) {
02350 iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
02351 std::ifstream adaptTempRead(simParams->adaptTempInFile);
02352 if (adaptTempRead) {
02353 int readInt;
02354 BigReal readReal;
02355 bool readBool;
02356
02357 adaptTempRead >> readInt;
02358
02359 adaptTempRead >> adaptTempT;
02360 adaptTempRead >> adaptTempBetaMin;
02361 adaptTempRead >> adaptTempBetaMax;
02362 adaptTempBetaMin = 1./adaptTempBetaMin;
02363 adaptTempBetaMax = 1./adaptTempBetaMax;
02364
02365 if (adaptTempBetaMin > adaptTempBetaMax){
02366 readReal = adaptTempBetaMax;
02367 adaptTempBetaMax = adaptTempBetaMin;
02368 adaptTempBetaMin = adaptTempBetaMax;
02369 }
02370 adaptTempRead >> adaptTempBins;
02371 adaptTempRead >> adaptTempCg;
02372 adaptTempRead >> adaptTempDt;
02373 adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02374 adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02375 adaptTempPotEnergySamples = new int[adaptTempBins];
02376 adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02377 adaptTempPotEnergyVar = new BigReal[adaptTempBins];
02378 adaptTempPotEnergyAve = new BigReal[adaptTempBins];
02379 adaptTempBetaN = new BigReal[adaptTempBins];
02380 adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02381 for(int j = 0; j < adaptTempBins; ++j) {
02382 adaptTempRead >> adaptTempPotEnergyAve[j];
02383 adaptTempRead >> adaptTempPotEnergyVar[j];
02384 adaptTempRead >> adaptTempPotEnergySamples[j];
02385 adaptTempRead >> adaptTempPotEnergyAveNum[j];
02386 adaptTempRead >> adaptTempPotEnergyVarNum[j];
02387 adaptTempRead >> adaptTempPotEnergyAveDen[j];
02388 adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02389 }
02390 adaptTempRead.close();
02391 }
02392 else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
02393 }
02394 else {
02395 adaptTempBins = simParams->adaptTempBins;
02396 adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02397 adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02398 adaptTempPotEnergySamples = new int[adaptTempBins];
02399 adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02400 adaptTempPotEnergyVar = new BigReal[adaptTempBins];
02401 adaptTempPotEnergyAve = new BigReal[adaptTempBins];
02402 adaptTempBetaN = new BigReal[adaptTempBins];
02403 adaptTempBetaMax = 1./simParams->adaptTempTmin;
02404 adaptTempBetaMin = 1./simParams->adaptTempTmax;
02405 adaptTempCg = simParams->adaptTempCgamma;
02406 adaptTempDt = simParams->adaptTempDt;
02407 adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02408 adaptTempT = simParams->initialTemp;
02409 if (simParams->langevinOn)
02410 adaptTempT = simParams->langevinTemp;
02411 else if (simParams->rescaleFreq > 0)
02412 adaptTempT = simParams->rescaleTemp;
02413 for(int j = 0; j < adaptTempBins; ++j){
02414 adaptTempPotEnergyAveNum[j] = 0.;
02415 adaptTempPotEnergyAveDen[j] = 0.;
02416 adaptTempPotEnergySamples[j] = 0;
02417 adaptTempPotEnergyVarNum[j] = 0.;
02418 adaptTempPotEnergyVar[j] = 0.;
02419 adaptTempPotEnergyAve[j] = 0.;
02420 adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02421 }
02422 }
02423 if (simParams->adaptTempAutoDt > 0.0) {
02424 adaptTempAutoDt = true;
02425 adaptTempDtMin = simParams->adaptTempAutoDt - 0.01;
02426 adaptTempDtMax = simParams->adaptTempAutoDt + 0.01;
02427 }
02428 adaptTempDTave = 0;
02429 adaptTempDTavenum = 0;
02430 iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
02431 iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
02432 iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
02433 iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
02434 iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
02435 if (simParams->adaptTempRestartFreq > 0) {
02436 NAMD_backup_file(simParams->adaptTempRestartFile);
02437 adaptTempRestartFile.open(simParams->adaptTempRestartFile);
02438 if(!adaptTempRestartFile)
02439 NAMD_die("Error opening restart file");
02440 }
02441 }
02442
02443 void Controller::adaptTempWriteRestart(int step) {
02444 if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
02445 adaptTempRestartFile.seekp(std::ios::beg);
02446 iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
02447 adaptTempRestartFile << step << " ";
02448
02449 adaptTempRestartFile << adaptTempT<< " ";
02450 adaptTempRestartFile << 1./adaptTempBetaMin << " ";
02451 adaptTempRestartFile << 1./adaptTempBetaMax << " ";
02452 adaptTempRestartFile << adaptTempBins << " ";
02453 adaptTempRestartFile << adaptTempCg << " ";
02454 adaptTempRestartFile << adaptTempDt ;
02455 adaptTempRestartFile << "\n" ;
02456 for(int j = 0; j < adaptTempBins; ++j) {
02457 adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
02458 adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
02459 adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
02460 adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
02461 adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
02462 adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
02463 adaptTempRestartFile << "\n";
02464 }
02465 adaptTempRestartFile.flush();
02466 }
02467 }
02468
02469 void Controller::adaptTempUpdate(int step, int minimize)
02470 {
02471
02472 if ( !simParams->adaptTempOn ) return;
02473 int j = 0;
02474 if (step == simParams->firstTimestep) {
02475 adaptTempInit(step);
02476 return;
02477 }
02478 if ( minimize || (step < simParams->adaptTempFirstStep ) ||
02479 ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
02480 const int adaptTempOutFreq = simParams->adaptTempOutFreq;
02481 const bool adaptTempDebug = simParams->adaptTempDebug;
02482
02483 BigReal adaptTempBeta = 1./adaptTempT;
02484 adaptTempBin = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
02485
02486 if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
02487 iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin
02488 << " adaptTempBeta: " << adaptTempBeta
02489 << " adaptTempDBeta: " << adaptTempDBeta
02490 << " betaMin:" << adaptTempBetaMin
02491 << " betaMax: " << adaptTempBetaMax << "\n";
02492 adaptTempPotEnergySamples[adaptTempBin] += 1;
02493 BigReal gammaAve = 1.-adaptTempCg/adaptTempPotEnergySamples[adaptTempBin];
02494
02495 BigReal potentialEnergy;
02496 BigReal potEnergyAverage;
02497 BigReal potEnergyVariance;
02498 potentialEnergy = totalEnergy-kineticEnergy;
02499
02500
02501 adaptTempPotEnergyAveNum[adaptTempBin] = adaptTempPotEnergyAveNum[adaptTempBin]*gammaAve + potentialEnergy;
02502 adaptTempPotEnergyAveDen[adaptTempBin] = adaptTempPotEnergyAveDen[adaptTempBin]*gammaAve + 1;
02503 adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
02504
02505 potEnergyAverage = adaptTempPotEnergyAveNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
02506 potEnergyVariance = adaptTempPotEnergyVarNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
02507 potEnergyVariance -= potEnergyAverage*potEnergyAverage;
02508
02509 adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
02510 adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
02511
02512
02513
02514
02515
02516 if ( ! ( step % simParams->adaptTempFreq ) ) {
02517
02518 if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
02519
02520 BigReal deltaBeta = 0.04*adaptTempBeta;
02521 BigReal betaPlus;
02522 BigReal betaMinus;
02523 int nMinus =0;
02524 int nPlus = 0;
02525 if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
02526 if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
02527 betaMinus = adaptTempBeta - deltaBeta;
02528 betaPlus = adaptTempBeta + deltaBeta;
02529 BigReal betaMinus2 = betaMinus*betaMinus;
02530 BigReal betaPlus2 = betaPlus*betaPlus;
02531 nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
02532 nPlus = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
02533
02534 BigReal potEnergyAve0 = 0.0;
02535 BigReal potEnergyAve1 = 0.0;
02536
02537 BigReal A0 = 0;
02538 BigReal A1 = 0;
02539 BigReal A2 = 0;
02540
02541 BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1];
02542 BigReal DeltaE2Ave;
02543 BigReal DeltaE2AveSum = 0;
02544 BigReal betaj;
02545 for (j = nMinus+1; j <= adaptTempBin; ++j) {
02546 potEnergyAve0 += adaptTempPotEnergyAve[j];
02547 DeltaE2Ave = adaptTempPotEnergyVar[j];
02548 DeltaE2AveSum += DeltaE2Ave;
02549 A0 += j*DeltaE2Ave;
02550 }
02551 A0 *= adaptTempDBeta;
02552 A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
02553 A0 *= adaptTempDBeta;
02554 betaj = adaptTempBetaN[nMinus+1]-betaMinus;
02555 betaj *= betaj;
02556 A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
02557 A0 /= (betaNp1-betaMinus);
02558
02559
02560 DeltaE2AveSum = 0;
02561 for (j = adaptTempBin+1; j < nPlus; j++) {
02562 potEnergyAve1 += adaptTempPotEnergyAve[j];
02563 DeltaE2Ave = adaptTempPotEnergyVar[j];
02564 DeltaE2AveSum += DeltaE2Ave;
02565 A1 += j*DeltaE2Ave;
02566 }
02567 A1 *= adaptTempDBeta;
02568 A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
02569 A1 *= adaptTempDBeta;
02570 if ( nPlus < adaptTempBins ) {
02571 betaj = betaPlus-adaptTempBetaN[nPlus];
02572 betaj *= betaj;
02573 A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
02574 }
02575 A1 /= (betaPlus-betaNp1);
02576
02577
02578 A2 = 0.5*adaptTempDBeta*potEnergyVariance;
02579
02580
02581 DeltaE2Ave = A0-A1;
02582 BigReal aplus = 0;
02583 BigReal aminus = 0;
02584 if (DeltaE2Ave != 0) {
02585 aplus = (A0-A2)/(A0-A1);
02586 if (aplus < 0) {
02587 aplus = 0;
02588 }
02589 if (aplus > 1) {
02590 aplus = 1;
02591 }
02592 aminus = 1-aplus;
02593 potEnergyAve0 *= adaptTempDBeta;
02594 potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
02595 potEnergyAve0 /= (betaNp1-betaMinus);
02596 potEnergyAve1 *= adaptTempDBeta;
02597 if ( nPlus < adaptTempBins ) {
02598 potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
02599 }
02600 potEnergyAve1 /= (betaPlus-betaNp1);
02601 potEnergyAverage = aminus*potEnergyAve0;
02602 potEnergyAverage += aplus*potEnergyAve1;
02603 }
02604 if (simParams->adaptTempDebug) {
02605 iout << "ADAPTEMP DEBUG:" << "\n"
02606 << " adaptTempBin: " << adaptTempBin << "\n"
02607 << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
02608 << " adaptTempBeta: " << adaptTempBeta << "\n"
02609 << " adaptTemp: " << adaptTempT<< "\n"
02610 << " betaMin: " << adaptTempBetaMin << "\n"
02611 << " betaMax: " << adaptTempBetaMax << "\n"
02612 << " gammaAve: " << gammaAve << "\n"
02613 << " deltaBeta: " << deltaBeta << "\n"
02614 << " betaMinus: " << betaMinus << "\n"
02615 << " betaPlus: " << betaPlus << "\n"
02616 << " nMinus: " << nMinus << "\n"
02617 << " nPlus: " << nPlus << "\n"
02618 << " A0: " << A0 << "\n"
02619 << " A1: " << A1 << "\n"
02620 << " A2: " << A2 << "\n"
02621 << " a+: " << aplus << "\n"
02622 << " a-: " << aminus << "\n"
02623 << endi;
02624 }
02625 }
02626 else {
02627 if (simParams->adaptTempDebug) {
02628 iout << "ADAPTEMP DEBUG:" << "\n"
02629 << " adaptTempBin: " << adaptTempBin << "\n"
02630 << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
02631 << " adaptTempBeta: " << adaptTempBeta << "\n"
02632 << " adaptTemp: " << adaptTempT<< "\n"
02633 << " betaMin: " << adaptTempBetaMin << "\n"
02634 << " betaMax: " << adaptTempBetaMax << "\n"
02635 << " gammaAve: " << gammaAve << "\n"
02636 << " deltaBeta: N/A\n"
02637 << " betaMinus: N/A\n"
02638 << " betaPlus: N/A\n"
02639 << " nMinus: N/A\n"
02640 << " nPlus: N/A\n"
02641 << " A0: N/A\n"
02642 << " A1: N/A\n"
02643 << " A2: N/A\n"
02644 << " a+: N/A\n"
02645 << " a-: N/A\n"
02646 << endi;
02647 }
02648 }
02649
02650
02651 BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
02652 dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
02653 dT += adaptTempT;
02654
02655
02656
02657 if ( dT > 1./adaptTempBetaMin || dT < 1./adaptTempBetaMax ) {
02658 dT = ((potentialEnergy-adaptTempPotEnergyAve[adaptTempBin])/BOLTZMANN+adaptTempT)*adaptTempDt;
02659 dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
02660 dT += adaptTempT;
02661
02662 if ( dT > 1./adaptTempBetaMin ) {
02663 if (!simParams->adaptTempRandom) {
02664
02665
02666
02667
02668 dT = adaptTempT;
02669 }
02670 else {
02671
02672
02673
02674 dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
02675 dT = 1./dT;
02676 }
02677 }
02678 else if ( dT < 1./adaptTempBetaMax ) {
02679 if (!simParams->adaptTempRandom) {
02680
02681
02682
02683 dT = adaptTempT;
02684 }
02685 else {
02686
02687
02688
02689 dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
02690 dT = 1./dT;
02691 }
02692 }
02693 else if (adaptTempAutoDt) {
02694
02695
02696 BigReal adaptTempTdiff = fabs(dT-adaptTempT);
02697 if (adaptTempTdiff > 0) {
02698 adaptTempDTave += adaptTempTdiff;
02699 adaptTempDTavenum++;
02700
02701 }
02702 if(adaptTempDTavenum == 100){
02703 BigReal Frac;
02704 adaptTempDTave /= adaptTempDTavenum;
02705 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
02706 Frac = adaptTempDTave/Frac;
02707
02708
02709 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
02710 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
02711 Frac = adaptTempDtMax/Frac;
02712 iout << "ADAPTEMP: Updating adaptTempDt to ";
02713 adaptTempDt *= Frac;
02714 iout << adaptTempDt << "\n" << endi;
02715 }
02716 adaptTempDTave = 0;
02717 adaptTempDTavenum = 0;
02718 }
02719 }
02720 }
02721 else if (adaptTempAutoDt) {
02722
02723
02724 BigReal adaptTempTdiff = fabs(dT-adaptTempT);
02725 if (adaptTempTdiff > 0) {
02726 adaptTempDTave += adaptTempTdiff;
02727 adaptTempDTavenum++;
02728
02729 }
02730 if(adaptTempDTavenum == 100){
02731 BigReal Frac;
02732 adaptTempDTave /= adaptTempDTavenum;
02733 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
02734 Frac = adaptTempDTave/Frac;
02735
02736
02737 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
02738 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
02739 Frac = adaptTempDtMax/Frac;
02740 iout << "ADAPTEMP: Updating adaptTempDt to ";
02741 adaptTempDt *= Frac;
02742 iout << adaptTempDt << "\n" << endi;
02743 }
02744 adaptTempDTave = 0;
02745 adaptTempDTavenum = 0;
02746
02747 }
02748
02749 }
02750
02751 adaptTempT = dT;
02752 broadcast->adaptTemperature.publish(step,adaptTempT);
02753 }
02754 adaptTempWriteRestart(step);
02755 if ( ! (step % adaptTempOutFreq) ) {
02756 iout << "ADAPTEMP: STEP " << step
02757 << " TEMP " << adaptTempT
02758 << " ENERGY " << std::setprecision(10) << potentialEnergy
02759 << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
02760 << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
02761 iout << "\n" << endi;
02762 }
02763
02764 }
02765
02766
02767 void Controller::compareChecksums(int step, int forgiving) {
02768 Node *node = Node::Object();
02769 Molecule *molecule = node->molecule;
02770
02771
02772 BigReal checksum, checksum_b;
02773 char errmsg[256];
02774
02775 checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
02776 if ( ((int)checksum) != molecule->numAtoms ) {
02777 sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
02778 (int)checksum, molecule->numAtoms);
02779 NAMD_bug(errmsg);
02780 }
02781
02782 checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
02783 if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
02784 if ( ! computeChecksum ) {
02785 computesPartitioned = 0;
02786 } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
02787 computesPartitioned = 1;
02788 } else {
02789 NAMD_bug("Bad global compute count!\n");
02790 }
02791 computeChecksum = ((int)checksum);
02792 }
02793
02794 checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
02795 if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
02796 sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
02797 (int)checksum, molecule->numCalcBonds);
02798 if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
02799 iout << iWARN << errmsg << endi;
02800 else NAMD_bug(errmsg);
02801 }
02802
02803 checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
02804 if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
02805 sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
02806 (int)checksum, molecule->numCalcAngles);
02807 if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
02808 iout << iWARN << errmsg << endi;
02809 else NAMD_bug(errmsg);
02810 }
02811
02812 checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
02813 if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
02814 sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
02815 (int)checksum, molecule->numCalcDihedrals);
02816 if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
02817 iout << iWARN << errmsg << endi;
02818 else NAMD_bug(errmsg);
02819 }
02820
02821 checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
02822 if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
02823 sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
02824 (int)checksum, molecule->numCalcImpropers);
02825 if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
02826 iout << iWARN << errmsg << endi;
02827 else NAMD_bug(errmsg);
02828 }
02829
02830 checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
02831 if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
02832 sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
02833 (int)checksum, molecule->numCalcTholes);
02834 if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
02835 iout << iWARN << errmsg << endi;
02836 else NAMD_bug(errmsg);
02837 }
02838
02839 checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
02840 if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
02841 sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
02842 (int)checksum, molecule->numCalcAnisos);
02843 if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
02844 iout << iWARN << errmsg << endi;
02845 else NAMD_bug(errmsg);
02846 }
02847
02848 checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
02849 if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
02850 sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
02851 (int)checksum, molecule->numCalcCrossterms);
02852 if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
02853 iout << iWARN << errmsg << endi;
02854 else NAMD_bug(errmsg);
02855 }
02856
02857 checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
02858 if ( ((int64)checksum) > molecule->numCalcExclusions &&
02859 ( ! simParams->mollyOn || step % slowFreq ) ) {
02860 if ( forgiving )
02861 iout << iWARN << "High global exclusion count ("
02862 << ((int64)checksum) << " vs "
02863 << molecule->numCalcExclusions << "), possible error!\n"
02864 << iWARN << "This warning is not unusual during minimization.\n"
02865 << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02866 else {
02867 char errmsg[256];
02868 sprintf(errmsg, "High global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02869 (long long)checksum, (long long)molecule->numCalcExclusions);
02870 NAMD_bug(errmsg);
02871 }
02872 }
02873
02874 if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
02875 ! simParams->goGroPair && ! simParams->qmForcesOn) {
02876 if ( forgiving || ! simParams->fullElectFrequency ) {
02877 iout << iWARN << "Low global exclusion count! ("
02878 << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
02879 if ( forgiving ) {
02880 iout << "\n"
02881 << iWARN << "This warning is not unusual during minimization.\n"
02882 << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02883 } else {
02884 iout << " System unstable or pairlistdist or cutoff too small.\n";
02885 }
02886 iout << endi;
02887 } else {
02888 char errmsg[256];
02889 sprintf(errmsg, "Low global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
02890 (long long)checksum, (long long)molecule->numCalcExclusions);
02891 NAMD_bug(errmsg);
02892 }
02893 }
02894
02895 #ifdef NAMD_CUDA
02896 checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
02897 if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
02898 ( ! simParams->mollyOn || step % slowFreq ) ) {
02899 if ( forgiving )
02900 iout << iWARN << "High global CUDA exclusion count ("
02901 << ((int64)checksum) << " vs "
02902 << molecule->numCalcFullExclusions << "), possible error!\n"
02903 << iWARN << "This warning is not unusual during minimization.\n"
02904 << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02905 else {
02906 char errmsg[256];
02907 sprintf(errmsg, "High global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02908 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02909 NAMD_bug(errmsg);
02910 }
02911 }
02912
02913 if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
02914 ! simParams->goGroPair && ! simParams->qmForcesOn) {
02915 if ( forgiving || ! simParams->fullElectFrequency ) {
02916 iout << iWARN << "Low global CUDA exclusion count! ("
02917 << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
02918 if ( forgiving ) {
02919 iout << "\n"
02920 << iWARN << "This warning is not unusual during minimization.\n"
02921 << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02922 } else {
02923 iout << " System unstable or pairlistdist or cutoff too small.\n";
02924 }
02925 iout << endi;
02926 } else {
02927 char errmsg[256];
02928 sprintf(errmsg, "Low global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
02929 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02930 NAMD_bug(errmsg);
02931 }
02932 }
02933 #endif
02934
02935 checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
02936 if ( ((int)checksum) && ! marginViolations ) {
02937 iout << iERROR << "Margin is too small for " << ((int)checksum) <<
02938 " atoms during timestep " << step << ".\n" << iERROR <<
02939 "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
02940 }
02941 marginViolations += (int)checksum;
02942
02943 checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
02944 if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
02945 iout << iINFO <<
02946 "Pairlistdist is too small for " << ((int)checksum) <<
02947 " computes during timestep " << step << ".\n" << endi;
02948 }
02949 if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
02950
02951 checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
02952 if ( checksum ) {
02953 if ( forgiving )
02954 iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
02955 else NAMD_bug("Stray PME grid charges detected!\n");
02956 }
02957 }
02958
02959 void Controller::printTiming(int step) {
02960
02961 if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
02962 {
02963 const double endWTime = CmiWallTimer() - firstWTime;
02964 const double endCTime = CmiTimer() - firstCTime;
02965
02966
02967 if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
02968
02969 const double elapsedW =
02970 (endWTime - startWTime) / simParams->outputTiming;
02971 const double elapsedC =
02972 (endCTime - startCTime) / simParams->outputTiming;
02973
02974 const double remainingW = elapsedW * (simParams->N - step);
02975 const double remainingW_hours = remainingW / 3600;
02976
02977 startWTime = endWTime;
02978 startCTime = endCTime;
02979
02980 #ifdef NAMD_CUDA
02981 if ( simParams->outputEnergies < 60 &&
02982 step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
02983 iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
02984 }
02985 #endif
02986 if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
02987 CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
02988 ", %g hours remaining, %f MB of memory in use.\n",
02989 step, endCTime, elapsedC, endWTime, elapsedW,
02990 remainingW_hours, memusage_MB());
02991 if ( fflush_count ) { --fflush_count; fflush(stdout); }
02992 }
02993 }
02994 }
02995
02996 void Controller::printMinimizeEnergies(int step) {
02997
02998 rescaleaccelMD(step,1);
02999 receivePressure(step,1);
03000
03001 Node *node = Node::Object();
03002 Molecule *molecule = node->molecule;
03003 compareChecksums(step,1);
03004
03005 printEnergies(step,1);
03006
03007 min_energy = totalEnergy;
03008 min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
03009 min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
03010 min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
03011 min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
03012
03013 }
03014
03015 void Controller::printDynamicsEnergies(int step) {
03016
03017 Node *node = Node::Object();
03018 Molecule *molecule = node->molecule;
03019 SimParameters *simParameters = node->simParameters;
03020 Lattice &lattice = state->lattice;
03021
03022 compareChecksums(step);
03023
03024 printEnergies(step,0);
03025 }
03026
03027 void Controller::printEnergies(int step, int minimize)
03028 {
03029 Node *node = Node::Object();
03030 Molecule *molecule = node->molecule;
03031 SimParameters *simParameters = node->simParameters;
03032 Lattice &lattice = state->lattice;
03033
03034
03035
03036
03037 BigReal bondEnergy;
03038 BigReal angleEnergy;
03039 BigReal dihedralEnergy;
03040 BigReal improperEnergy;
03041 BigReal crosstermEnergy;
03042 BigReal boundaryEnergy;
03043 BigReal miscEnergy;
03044 BigReal potentialEnergy;
03045 BigReal flatEnergy;
03046 BigReal smoothEnergy;
03047 BigReal work;
03048
03049 Vector momentum;
03050 Vector angularMomentum;
03051 BigReal volume = lattice.volume();
03052
03053 bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
03054 angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
03055 dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
03056 improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
03057 crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
03058 boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
03059 miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
03060
03061 if ( minimize || ! ( step % nbondFreq ) )
03062 {
03063 electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
03064 ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
03065
03066
03067 groLJEnergy = reduction->item(REDUCTION_GRO_LJ_ENERGY);
03068 groGaussEnergy = reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
03069
03070 goNativeEnergy = reduction->item(REDUCTION_GO_NATIVE_ENERGY);
03071 goNonnativeEnergy = reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
03072 goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
03073
03074
03075 bondedEnergyDiff_f = reduction->item(REDUCTION_BONDED_ENERGY_F);
03076 electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
03077 ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
03078
03079 bondedEnergy_ti_1 = reduction->item(REDUCTION_BONDED_ENERGY_TI_1);
03080 electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
03081 ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
03082 bondedEnergy_ti_2 = reduction->item(REDUCTION_BONDED_ENERGY_TI_2);
03083 electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
03084 ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
03085
03086 }
03087
03088 if ( minimize || ! ( step % slowFreq ) )
03089 {
03090 electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
03091
03092 electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
03093
03094 electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
03095 electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
03096 electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
03097 electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
03098
03099 }
03100
03101 if (simParameters->LJcorrection && volume) {
03102 #ifdef MEM_OPT_VERSION
03103 NAMD_bug("LJcorrection not supported in memory optimized build.");
03104 #else
03105
03106 BigReal alchLambda = simParameters->getCurrentLambda(step);
03107 ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
03108
03109 if (simParameters->alchOn) {
03110 if (simParameters->alchFepOn) {
03111 BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
03112 ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
03113 } else if (simParameters->alchThermIntOn) {
03114 ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
03115 ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
03116 }
03117 }
03118 #endif
03119 }
03120
03121
03122
03123 if (simParameters->alchLambdaFreq > 0) {
03124 if (step <=
03125 simParameters->firstTimestep + simParameters->alchEquilSteps) {
03126 cumAlchWork = 0.0;
03127 }
03128 alchWork = computeAlchWork(step);
03129 cumAlchWork += alchWork;
03130 }
03131
03132
03133 momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
03134 momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
03135 momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
03136 angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
03137 angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
03138 angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
03139
03140
03141 potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
03142 + improperEnergy + electEnergy + electEnergySlow + ljEnergy
03143 + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy
03144 + groLJEnergy + groGaussEnergy);
03145
03146 totalEnergy = potentialEnergy + kineticEnergy;
03147 flatEnergy = (totalEnergy
03148 + (1.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
03149 if ( !(step%slowFreq) ) {
03150
03151 BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
03152 if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
03153 if ( step != simParams->firstTimestep ) {
03154 smooth2_avg *= 0.9375;
03155 smooth2_avg += 0.0625 * s;
03156 }
03157 }
03158 smoothEnergy = (flatEnergy + smooth2_avg
03159 - (4.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
03160
03161
03162 if (step <= simParameters->firstTimestep &&
03163 (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
03164 heat = 0.0;
03165 totalEnergy0 = totalEnergy;
03166 }
03167 if ( simParameters->outputMomenta && ! minimize &&
03168 ! ( step % simParameters->outputMomenta ) )
03169 {
03170 iout << "MOMENTUM: " << step
03171 << " P: " << momentum
03172 << " L: " << angularMomentum
03173 << "\n" << endi;
03174 }
03175
03176 if ( simParameters->outputPressure ) {
03177 pressure_tavg += pressure;
03178 groupPressure_tavg += groupPressure;
03179 tavg_count += 1;
03180 if ( minimize || ! ( step % simParameters->outputPressure ) ) {
03181 iout << "PRESSURE: " << step << " "
03182 << PRESSUREFACTOR * pressure << "\n"
03183 << "GPRESSURE: " << step << " "
03184 << PRESSUREFACTOR * groupPressure << "\n";
03185 if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
03186 << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
03187 << "GPRESSAVG: " << step << " "
03188 << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
03189 iout << endi;
03190 pressure_tavg = 0;
03191 groupPressure_tavg = 0;
03192 tavg_count = 0;
03193 }
03194 }
03195
03196
03197 if (pressureProfileSlabs) {
03198 const int freq = simParams->pressureProfileFreq;
03199 const int arraysize = 3*pressureProfileSlabs;
03200
03201 BigReal *total = new BigReal[arraysize];
03202 memset(total, 0, arraysize*sizeof(BigReal));
03203 const int first = simParams->firstTimestep;
03204
03205 if (ppbonded) ppbonded->getData(first, step, lattice, total);
03206 if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
03207 if (ppint) ppint->getData(first, step, lattice, total);
03208 for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
03209 pressureProfileCount++;
03210
03211 if (!(step % freq)) {
03212
03213 BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
03214
03215 iout << "PRESSUREPROFILE: " << step << " ";
03216 if (step == first) {
03217
03218 for (int i=0; i<arraysize; i++) {
03219 iout << total[i] * scalefac << " ";
03220 }
03221 } else {
03222
03223 scalefac /= pressureProfileCount;
03224 for (int i=0; i<arraysize; i++)
03225 iout << pressureProfileAverage[i]*scalefac << " ";
03226 }
03227 iout << "\n" << endi;
03228
03229
03230 memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
03231 pressureProfileCount = 0;
03232 }
03233 delete [] total;
03234 }
03235
03236 if ( step != simParams->firstTimestep || stepInFullRun == 0 ) {
03237 if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
03238 int benchPhase = stepInFullRun / simParams->firstLdbStep;
03239 if ( benchPhase > 0 && benchPhase < 10 ) {
03240 #ifdef NAMD_CUDA
03241 if ( simParams->outputEnergies < 60 ) {
03242 iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
03243 }
03244 #endif
03245 iout << iINFO;
03246 if ( benchPhase < 4 ) iout << "Initial time: ";
03247 else iout << "Benchmark time: ";
03248 iout << CkNumPes() << " CPUs ";
03249 {
03250 BigReal wallPerStep =
03251 (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
03252 BigReal ns = simParams->dt / 1000000.0;
03253 BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
03254 BigReal daysPerNano = wallPerStep * days / ns;
03255 iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
03256 iout << memusage_MB() << " MB memory\n" << endi;
03257 }
03258 }
03259 startBenchTime = CmiWallTimer();
03260 }
03261 ++stepInFullRun;
03262 }
03263
03264 printTiming(step);
03265
03266 Vector pairVDWForce, pairElectForce;
03267 if ( simParameters->pairInteractionOn ) {
03268 GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
03269 GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
03270 }
03271
03272
03273 if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03274 work = totalEnergy - totalEnergy0 - heat;
03275 }
03276
03277
03278 #ifdef NAMD_TCL
03279 #define CALLBACKDATA(LABEL,VALUE) \
03280 labels << (LABEL) << " "; values << (VALUE) << " ";
03281 #define CALLBACKLIST(LABEL,VALUE) \
03282 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
03283 if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
03284 std::ostringstream labels, values;
03285 #if CMK_BLUEGENEL
03286
03287 values.precision(16);
03288 #else
03289 values << std::setprecision(16);
03290 #endif
03291 CALLBACKDATA("TS",step);
03292 CALLBACKDATA("BOND",bondEnergy);
03293 CALLBACKDATA("ANGLE",angleEnergy);
03294 CALLBACKDATA("DIHED",dihedralEnergy);
03295 CALLBACKDATA("CROSS",crosstermEnergy);
03296 CALLBACKDATA("IMPRP",improperEnergy);
03297 CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
03298 CALLBACKDATA("VDW",ljEnergy);
03299 CALLBACKDATA("BOUNDARY",boundaryEnergy);
03300 CALLBACKDATA("MISC",miscEnergy);
03301 CALLBACKDATA("POTENTIAL",potentialEnergy);
03302 CALLBACKDATA("KINETIC",kineticEnergy);
03303 CALLBACKDATA("TOTAL",totalEnergy);
03304 CALLBACKDATA("TEMP",temperature);
03305 CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
03306 CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
03307 CALLBACKDATA("VOLUME",lattice.volume());
03308 CALLBACKLIST("CELL_A",lattice.a());
03309 CALLBACKLIST("CELL_B",lattice.b());
03310 CALLBACKLIST("CELL_C",lattice.c());
03311 CALLBACKLIST("CELL_O",lattice.origin());
03312 labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
03313 << lattice.b_p() << " " << lattice.c_p() << "} ";
03314 if (simParameters->drudeOn) {
03315 CALLBACKDATA("DRUDEBOND",drudeBondTemp);
03316 }
03317 if ( simParameters->pairInteractionOn ) {
03318 CALLBACKLIST("VDW_FORCE",pairVDWForce);
03319 CALLBACKLIST("ELECT_FORCE",pairElectForce);
03320 }
03321 if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03322 CALLBACKLIST("HEAT",heat);
03323 CALLBACKLIST("WORK",work);
03324 }
03325 if (simParameters->alchOn) {
03326 if (simParameters->alchThermIntOn) {
03327 CALLBACKLIST("BOND1", bondedEnergy_ti_1);
03328 CALLBACKLIST("ELEC1", electEnergy_ti_1 + electEnergySlow_ti_1 +
03329 electEnergyPME_ti_1);
03330 CALLBACKLIST("VDW1", ljEnergy_ti_1);
03331 CALLBACKLIST("BOND2", bondedEnergy_ti_2);
03332 CALLBACKLIST("ELEC2", electEnergy_ti_2 + electEnergySlow_ti_2 +
03333 electEnergyPME_ti_2);
03334 CALLBACKLIST("VDW2", ljEnergy_ti_2);
03335 if (simParameters->alchLambdaFreq > 0) {
03336 CALLBACKLIST("CUMALCHWORK", cumAlchWork);
03337 }
03338 } else if (simParameters->alchFepOn) {
03339 CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
03340 improperEnergy + bondedEnergyDiff_f);
03341 CALLBACKLIST("ELEC2", electEnergy_f + electEnergySlow_f);
03342 CALLBACKLIST("VDW2", ljEnergy_f);
03343 }
03344 }
03345
03346 labels << '\0'; values << '\0';
03347 state->callback_labelstring = labels.str();
03348 state->callback_valuestring = values.str();
03349
03350 }
03351 #undef CALLBACKDATA
03352 #endif
03353
03354 drudeBondTempAvg += drudeBondTemp;
03355
03356 temp_avg += temperature;
03357 pressure_avg += trace(pressure)/3.;
03358 groupPressure_avg += trace(groupPressure)/3.;
03359 avg_count += 1;
03360
03361 if ( simParams->outputPairlists && pairlistWarnings &&
03362 ! (step % simParams->outputPairlists) ) {
03363 iout << iINFO << pairlistWarnings <<
03364 " pairlist warnings in past " << simParams->outputPairlists <<
03365 " steps.\n" << endi;
03366 pairlistWarnings = 0;
03367 }
03368
03369 BigReal enthalpy;
03370 if (simParameters->multigratorOn && ((step % simParameters->outputEnergies) == 0)) {
03371 enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
03372 }
03373
03374
03375 if ( ! minimize && step % simParameters->outputEnergies ) return;
03376
03377
03378 if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
03379 IMDEnergies energies;
03380 energies.tstep = step;
03381 energies.T = temp_avg/avg_count;
03382 energies.Etot = totalEnergy;
03383 energies.Epot = potentialEnergy;
03384 energies.Evdw = ljEnergy;
03385 energies.Eelec = electEnergy + electEnergySlow;
03386 energies.Ebond = bondEnergy;
03387 energies.Eangle = angleEnergy;
03388 energies.Edihe = dihedralEnergy + crosstermEnergy;
03389 energies.Eimpr = improperEnergy;
03390 Node::Object()->imd->gather_energies(&energies);
03391 }
03392
03393 if ( marginViolations ) {
03394 iout << iERROR << marginViolations <<
03395 " margin violations detected since previous energy output.\n" << endi;
03396 }
03397 marginViolations = 0;
03398
03399 if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
03400 {
03401 iout << "ETITLE: TS";
03402 iout << FORMAT("BOND");
03403 iout << FORMAT("ANGLE");
03404 iout << FORMAT("DIHED");
03405 if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
03406 iout << FORMAT("IMPRP");
03407 iout << " ";
03408 iout << FORMAT("ELECT");
03409 iout << FORMAT("VDW");
03410 iout << FORMAT("BOUNDARY");
03411 iout << FORMAT("MISC");
03412 iout << FORMAT("KINETIC");
03413 iout << " ";
03414 iout << FORMAT("TOTAL");
03415 iout << FORMAT("TEMP");
03416 iout << FORMAT("POTENTIAL");
03417
03418 iout << FORMAT("TOTAL3");
03419 iout << FORMAT("TEMPAVG");
03420 if ( volume != 0. ) {
03421 iout << " ";
03422 iout << FORMAT("PRESSURE");
03423 iout << FORMAT("GPRESSURE");
03424 iout << FORMAT("VOLUME");
03425 iout << FORMAT("PRESSAVG");
03426 iout << FORMAT("GPRESSAVG");
03427 }
03428 if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03429 iout << " ";
03430 iout << FORMAT("HEAT");
03431 iout << FORMAT("WORK");
03432 }
03433 if (simParameters->drudeOn) {
03434 iout << " ";
03435 iout << FORMAT("DRUDEBOND");
03436 iout << FORMAT("DRBONDAVG");
03437 }
03438
03439 if (simParameters->goGroPair) {
03440 iout << " ";
03441 iout << FORMAT("GRO_PAIR_LJ");
03442 iout << FORMAT("GRO_PAIR_GAUSS");
03443 }
03444
03445 if (simParameters->goForcesOn) {
03446 iout << " ";
03447 iout << FORMAT("NATIVE");
03448 iout << FORMAT("NONNATIVE");
03449
03450
03451 iout << FORMAT("GOTOTAL");
03452
03453 }
03454
03455
03456 if (simParameters->alchOn) {
03457 if (simParameters->alchThermIntOn) {
03458 iout << "\nTITITLE: TS";
03459 iout << FORMAT("BOND1");
03460 iout << FORMAT("ELECT1");
03461 iout << FORMAT("VDW1");
03462 iout << FORMAT("BOND2");
03463 iout << " ";
03464 iout << FORMAT("ELECT2");
03465 iout << FORMAT("VDW2");
03466 if (simParameters->alchLambdaFreq > 0) {
03467 iout << FORMAT("LAMBDA");
03468 iout << FORMAT("ALCHWORK");
03469 iout << FORMAT("CUMALCHWORK");
03470 }
03471 } else if (simParameters->alchFepOn) {
03472 iout << "\nFEPTITLE: TS";
03473 iout << FORMAT("BOND2");
03474 iout << FORMAT("ELECT2");
03475 iout << FORMAT("VDW2");
03476 if (simParameters->alchLambdaFreq > 0) {
03477 iout << FORMAT("LAMBDA");
03478 }
03479 }
03480 }
03481
03482 iout << "\n\n" << endi;
03483
03484 if (simParameters->qmForcesOn) {
03485 iout << "QMETITLE: TS";
03486 iout << FORMAT("QMID");
03487 iout << FORMAT("ENERGY");
03488 if (simParameters->PMEOn) iout << FORMAT("PMECORRENERGY");
03489 iout << "\n\n" << endi;
03490 }
03491
03492 }
03493
03494
03495
03496 iout << ETITLE(step);
03497 iout << FORMAT(bondEnergy);
03498 iout << FORMAT(angleEnergy);
03499 if ( simParameters->mergeCrossterms ) {
03500 iout << FORMAT(dihedralEnergy+crosstermEnergy);
03501 } else {
03502 iout << FORMAT(dihedralEnergy);
03503 iout << FORMAT(crosstermEnergy);
03504 }
03505 iout << FORMAT(improperEnergy);
03506 iout << " ";
03507 iout << FORMAT(electEnergy+electEnergySlow);
03508 iout << FORMAT(ljEnergy);
03509 iout << FORMAT(boundaryEnergy);
03510 iout << FORMAT(miscEnergy);
03511 iout << FORMAT(kineticEnergy);
03512 iout << " ";
03513 iout << FORMAT(totalEnergy);
03514 iout << FORMAT(temperature);
03515 iout << FORMAT(potentialEnergy);
03516
03517 iout << FORMAT(smoothEnergy);
03518 iout << FORMAT(temp_avg/avg_count);
03519 if ( volume != 0. )
03520 {
03521 iout << " ";
03522 iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
03523 iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
03524 iout << FORMAT(volume);
03525 iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
03526 iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
03527 }
03528 if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03529 iout << " ";
03530 iout << FORMAT(heat);
03531 iout << FORMAT(work);
03532 }
03533 if (simParameters->drudeOn) {
03534 iout << " ";
03535 iout << FORMAT(drudeBondTemp);
03536 iout << FORMAT(drudeBondTempAvg/avg_count);
03537 }
03538
03539 if (simParameters->goGroPair) {
03540 iout << " ";
03541 iout << FORMAT(groLJEnergy);
03542 iout << FORMAT(groGaussEnergy);
03543 }
03544
03545 if (simParameters->goForcesOn) {
03546 iout << " ";
03547 iout << FORMAT(goNativeEnergy);
03548 iout << FORMAT(goNonnativeEnergy);
03549
03550
03551 iout << FORMAT(goTotalEnergy);
03552
03553 }
03554
03555 if (simParameters->alchOn) {
03556 if (simParameters->alchThermIntOn) {
03557 iout << "\n";
03558 iout << TITITLE(step);
03559 iout << FORMAT(bondedEnergy_ti_1);
03560 iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 +
03561 electEnergyPME_ti_1);
03562 iout << FORMAT(ljEnergy_ti_1);
03563 iout << FORMAT(bondedEnergy_ti_2);
03564 iout << " ";
03565 iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 +
03566 electEnergyPME_ti_2);
03567 iout << FORMAT(ljEnergy_ti_2);
03568 if (simParameters->alchLambdaFreq > 0) {
03569 iout << FORMAT(simParameters->getCurrentLambda(step));
03570 iout << FORMAT(alchWork);
03571 iout << FORMAT(cumAlchWork);
03572 }
03573 } else if (simParameters->alchFepOn) {
03574 iout << "\n";
03575 iout << FEPTITLE2(step);
03576 iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy
03577 + improperEnergy + bondedEnergyDiff_f);
03578 iout << FORMAT(electEnergy_f + electEnergySlow_f);
03579 iout << FORMAT(ljEnergy_f);
03580 if (simParameters->alchLambdaFreq > 0) {
03581 iout << FORMAT(simParameters->getCurrentLambda(step));
03582 }
03583 }
03584 }
03585
03586 iout << "\n\n" << endi;
03587
03588 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
03589 char webout[80];
03590 sprintf(webout,"%d %d %d %d",(int)totalEnergy,
03591 (int)(potentialEnergy),
03592 (int)kineticEnergy,(int)temperature);
03593 CApplicationDepositNode0Data(webout);
03594 #endif
03595
03596 if (simParameters->pairInteractionOn) {
03597 iout << "PAIR INTERACTION:";
03598 iout << " STEP: " << step;
03599 iout << " VDW_FORCE: ";
03600 iout << FORMAT(pairVDWForce.x);
03601 iout << FORMAT(pairVDWForce.y);
03602 iout << FORMAT(pairVDWForce.z);
03603 iout << " ELECT_FORCE: ";
03604 iout << FORMAT(pairElectForce.x);
03605 iout << FORMAT(pairElectForce.y);
03606 iout << FORMAT(pairElectForce.z);
03607 iout << "\n" << endi;
03608 }
03609 drudeBondTempAvg = 0;
03610 temp_avg = 0;
03611 pressure_avg = 0;
03612 groupPressure_avg = 0;
03613 avg_count = 0;
03614
03615 if ( fflush_count ) {
03616 --fflush_count;
03617 fflush(stdout);
03618 }
03619 }
03620
03621 void Controller::writeExtendedSystemLabels(ofstream_namd &file) {
03622 Lattice &lattice = state->lattice;
03623 file << "#$LABELS step";
03624 if ( lattice.a_p() ) file << " a_x a_y a_z";
03625 if ( lattice.b_p() ) file << " b_x b_y b_z";
03626 if ( lattice.c_p() ) file << " c_x c_y c_z";
03627 file << " o_x o_y o_z";
03628 if ( simParams->langevinPistonOn ) {
03629 file << " s_x s_y s_z s_u s_v s_w";
03630 }
03631 file << std::endl;
03632 }
03633
03634 void Controller::writeExtendedSystemData(int step, ofstream_namd &file) {
03635 Lattice &lattice = state->lattice;
03636 file.precision(12);
03637 file << step;
03638 if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
03639 if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
03640 if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
03641 file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
03642 if ( simParams->langevinPistonOn ) {
03643 Vector strainRate = diagonal(langevinPiston_strainRate);
03644 Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
03645 file << " " << strainRate.x;
03646 file << " " << strainRate.y;
03647 file << " " << strainRate.z;
03648 file << " " << strainRate2.x;
03649 file << " " << strainRate2.y;
03650 file << " " << strainRate2.z;
03651 }
03652 file << std::endl;
03653 }
03654
03655 void Controller::enqueueCollections(int timestep)
03656 {
03657 if ( Output::coordinateNeeded(timestep) ){
03658 collection->enqueuePositions(timestep,state->lattice);
03659 }
03660 if ( Output::velocityNeeded(timestep) )
03661 collection->enqueueVelocities(timestep);
03662 if ( Output::forceNeeded(timestep) )
03663 collection->enqueueForces(timestep);
03664 }
03665
03666
03667 void Controller::outputFepEnergy(int step) {
03668 if (simParams->alchOn && simParams->alchFepOn) {
03669 const int stepInRun = step - simParams->firstTimestep;
03670 const int alchEquilSteps = simParams->alchEquilSteps;
03671 const BigReal alchLambda = simParams->alchLambda;
03672 const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03673
03674 if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
03675 FepNo = 0;
03676 exp_dE_ByRT = 0.0;
03677 net_dE = 0.0;
03678 }
03679 dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
03680 electEnergy - electEnergySlow - ljEnergy;
03681 BigReal RT = BOLTZMANN * simParams->alchTemp;
03682
03683 if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
03684 (step / simParams->alchIDWSFreq) % 2 == 1 )) {
03685
03686 FepNo++;
03687 exp_dE_ByRT += exp(-dE/RT);
03688 net_dE += dE;
03689 }
03690
03691 if (simParams->alchOutFreq) {
03692 if (stepInRun == 0) {
03693 if (!fepFile.is_open()) {
03694 NAMD_backup_file(simParams->alchOutFile);
03695 fepFile.open(simParams->alchOutFile);
03696 iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
03697 if(alchEnsembleAvg){
03698 fepSum = 0.0;
03699 fepFile << "# STEP Elec "
03700 << "vdW dE dE_avg Temp dG\n"
03701 << "# l l+dl "
03702 << " l l+dl E(l+dl)-E(l)" << std::endl;
03703 }
03704 else{
03705 fepFile << "# STEP Elec "
03706 << "vdW dE Temp\n"
03707 << "# l l+dl "
03708 << " l l+dl E(l+dl)-E(l)" << std::endl;
03709 }
03710 }
03711 if(!step){
03712 fepFile << "#NEW FEP WINDOW: "
03713 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
03714 << simParams->alchLambda2;
03715 if ( simParams->alchLambdaIDWS >= 0. ) {
03716 fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
03717 }
03718 fepFile << std::endl;
03719 }
03720 }
03721 if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03722 fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03723 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03724 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03725 }
03726 if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
03727 writeFepEnergyData(step, fepFile);
03728 fepFile.flush();
03729 }
03730 if (alchEnsembleAvg && (step == simParams->N)) {
03731 fepSum = fepSum + dG;
03732 fepFile << "#Free energy change for lambda window [ " << alchLambda
03733 << " " << simParams->alchLambda2 << " ] is " << dG
03734 << " ; net change until now is " << fepSum << std::endl;
03735 fepFile.flush();
03736 }
03737 }
03738 }
03739 }
03740
03741 void Controller::outputTiEnergy(int step) {
03742 if (simParams->alchOn && simParams->alchThermIntOn) {
03743 const int stepInRun = step - simParams->firstTimestep;
03744 const int alchEquilSteps = simParams->alchEquilSteps;
03745 const int alchLambdaFreq = simParams->alchLambdaFreq;
03746
03747 if (stepInRun == 0 || stepInRun == alchEquilSteps) {
03748 TiNo = 0;
03749 net_dEdl_bond_1 = 0;
03750 net_dEdl_bond_2 = 0;
03751 net_dEdl_elec_1 = 0;
03752 net_dEdl_elec_2 = 0;
03753 net_dEdl_lj_1 = 0;
03754 net_dEdl_lj_2 = 0;
03755 }
03756 TiNo++;
03757 net_dEdl_bond_1 += bondedEnergy_ti_1;
03758 net_dEdl_bond_2 += bondedEnergy_ti_2;
03759 net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
03760 + electEnergyPME_ti_1);
03761 net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
03762 + electEnergyPME_ti_2);
03763 net_dEdl_lj_1 += ljEnergy_ti_1;
03764 net_dEdl_lj_2 += ljEnergy_ti_2;
03765
03766
03767
03768 if (simParams->alchOutFreq) {
03769 if (stepInRun == 0 || stepInRun == alchEquilSteps
03770 || (! ((step - 1) % simParams->alchOutFreq))) {
03771
03772
03773 recent_TiNo = 0;
03774 recent_dEdl_bond_1 = 0;
03775 recent_dEdl_bond_2 = 0;
03776 recent_dEdl_elec_1 = 0;
03777 recent_dEdl_elec_2 = 0;
03778 recent_dEdl_lj_1 = 0;
03779 recent_dEdl_lj_2 = 0;
03780 recent_alchWork = 0;
03781 }
03782 recent_TiNo++;
03783 recent_dEdl_bond_1 += bondedEnergy_ti_1;
03784 recent_dEdl_bond_2 += bondedEnergy_ti_2;
03785 recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
03786 + electEnergyPME_ti_1);
03787 recent_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
03788 + electEnergyPME_ti_2);
03789 recent_dEdl_lj_1 += ljEnergy_ti_1;
03790 recent_dEdl_lj_2 += ljEnergy_ti_2;
03791 recent_alchWork += alchWork;
03792
03793 if (stepInRun == 0) {
03794 if (!tiFile.is_open()) {
03795 NAMD_backup_file(simParams->alchOutFile);
03796 tiFile.open(simParams->alchOutFile);
03797
03798
03799
03800
03801
03802
03803 iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
03804 tiFile << "#TITITLE: TS";
03805 tiFile << FORMAT("BOND1");
03806 tiFile << FORMAT("AVGBOND1");
03807 tiFile << FORMAT("ELECT1");
03808 tiFile << FORMAT("AVGELECT1");
03809 tiFile << " ";
03810 tiFile << FORMAT("VDW1");
03811 tiFile << FORMAT("AVGVDW1");
03812 tiFile << FORMAT("BOND2");
03813 tiFile << FORMAT("AVGBOND2");
03814 tiFile << FORMAT("ELECT2");
03815 tiFile << " ";
03816 tiFile << FORMAT("AVGELECT2");
03817 tiFile << FORMAT("VDW2");
03818 tiFile << FORMAT("AVGVDW2");
03819 if (alchLambdaFreq > 0) {
03820 tiFile << FORMAT("ALCHWORK");
03821 tiFile << FORMAT("CUMALCHWORK");
03822 }
03823 tiFile << std::endl;
03824 }
03825
03826 if (alchLambdaFreq > 0) {
03827 tiFile << "#ALCHEMICAL SWITCHING ACTIVE "
03828 << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
03829 << "\n#LAMBDA SCHEDULE: "
03830 << "dL: " << simParams->getLambdaDelta()
03831 << " Freq: " << alchLambdaFreq;
03832 }
03833 else {
03834 const BigReal alchLambda = simParams->alchLambda;
03835 const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
03836 const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
03837 const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
03838 const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
03839 const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
03840 const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
03841 tiFile << "#NEW TI WINDOW: "
03842 << "LAMBDA " << alchLambda
03843 << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
03844 << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
03845 << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
03846 << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
03847 }
03848 tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
03849 << std::endl;
03850 }
03851
03852 if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03853 tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03854 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03855 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03856 }
03857 if ((step%simParams->alchOutFreq) == 0) {
03858 writeTiEnergyData(step, tiFile);
03859 tiFile.flush();
03860 }
03861 }
03862 }
03863 }
03864
03865
03866
03867
03868
03869
03870 BigReal Controller::computeAlchWork(const int step) {
03871
03872 const BigReal oldLambda = simParams->getCurrentLambda(step-1);
03873 const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
03874 const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
03875 const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
03876 const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
03877 const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
03878 const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
03879
03880 const BigReal alchLambda = simParams->getCurrentLambda(step);
03881 const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
03882 const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
03883 const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
03884 const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
03885 const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
03886 const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda);
03887
03888 return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
03889 (elec_lambda_12 - elec_lambda_1)*
03890 (electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1) +
03891 (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
03892 (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
03893 (elec_lambda_22 - elec_lambda_2)*
03894 (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) +
03895 (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
03896 );
03897 }
03898
03899 void Controller::writeFepEnergyData(int step, ofstream_namd &file) {
03900 BigReal eeng = electEnergy + electEnergySlow;
03901 BigReal eeng_f = electEnergy_f + electEnergySlow_f;
03902 BigReal RT = BOLTZMANN * simParams->alchTemp;
03903
03904 const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03905 const int stepInRun = step - simParams->firstTimestep;
03906 const BigReal alchLambda = simParams->alchLambda;
03907
03908 if(stepInRun){
03909 if ( simParams->alchLambdaIDWS >= 0. &&
03910 (step / simParams->alchIDWSFreq) % 2 == 0 ) {
03911
03912 fepFile << FEPTITLE_BACK(step);
03913 } else {
03914
03915 fepFile << FEPTITLE(step);
03916 }
03917 fepFile << FORMAT(eeng);
03918 fepFile << FORMAT(eeng_f);
03919 fepFile << FORMAT(ljEnergy);
03920 fepFile << FORMAT(ljEnergy_f);
03921 fepFile << FORMAT(dE);
03922 if(alchEnsembleAvg){
03923 BigReal dE_avg = net_dE / FepNo;
03924 fepFile << FORMAT(dE_avg);
03925 }
03926 fepFile << FORMAT(temperature);
03927 if(alchEnsembleAvg){
03928 dG = -(RT * log(exp_dE_ByRT / FepNo));
03929 fepFile << FORMAT(dG);
03930 }
03931 fepFile << std::endl;
03932 }
03933 }
03934
03935 void Controller::writeTiEnergyData(int step, ofstream_namd &file) {
03936 tiFile << TITITLE(step);
03937 tiFile << FORMAT(recent_dEdl_bond_1 / recent_TiNo);
03938 tiFile << FORMAT(net_dEdl_bond_1 / TiNo);
03939 tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
03940 tiFile << FORMAT(net_dEdl_elec_1/TiNo);
03941 tiFile << " ";
03942 tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
03943 tiFile << FORMAT(net_dEdl_lj_1/TiNo);
03944 tiFile << FORMAT(recent_dEdl_bond_2 / recent_TiNo);
03945 tiFile << FORMAT(net_dEdl_bond_2 / TiNo);
03946 tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
03947 tiFile << " ";
03948 tiFile << FORMAT(net_dEdl_elec_2/TiNo);
03949 tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
03950 tiFile << FORMAT(net_dEdl_lj_2/TiNo);
03951 if (simParams->alchLambdaFreq > 0) {
03952 tiFile << FORMAT(recent_alchWork / recent_TiNo);
03953 tiFile << FORMAT(cumAlchWork);
03954 }
03955 tiFile << std::endl;
03956 }
03957
03958
03959 void Controller::outputExtendedSystem(int step)
03960 {
03961
03962 if ( step >= 0 ) {
03963
03964
03965 if ( simParams->xstFrequency &&
03966 ((step % simParams->xstFrequency) == 0) )
03967 {
03968 if ( ! xstFile.is_open() )
03969 {
03970 iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
03971 NAMD_backup_file(simParams->xstFilename);
03972 xstFile.open(simParams->xstFilename);
03973 while (!xstFile) {
03974 if ( errno == EINTR ) {
03975 CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
03976 xstFile.clear();
03977 xstFile.open(simParams->xstFilename);
03978 continue;
03979 }
03980 char err_msg[257];
03981 sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
03982 NAMD_err(err_msg);
03983 }
03984 xstFile << "# NAMD extended system trajectory file" << std::endl;
03985 writeExtendedSystemLabels(xstFile);
03986 }
03987 writeExtendedSystemData(step,xstFile);
03988 xstFile.flush();
03989 }
03990
03991
03992
03993 if ( simParams->restartFrequency &&
03994 ((step % simParams->restartFrequency) == 0) &&
03995 (step != simParams->firstTimestep) )
03996 {
03997 iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
03998 << step << "\n" << endi;
03999 char fname[140];
04000 const char *bsuffix = ".old";
04001 strcpy(fname, simParams->restartFilename);
04002 if ( simParams->restartSave ) {
04003 char timestepstr[20];
04004 sprintf(timestepstr,".%d",step);
04005 strcat(fname, timestepstr);
04006 bsuffix = ".BAK";
04007 }
04008 strcat(fname, ".xsc");
04009 NAMD_backup_file(fname,bsuffix);
04010 ofstream_namd xscFile(fname);
04011 while (!xscFile) {
04012 if ( errno == EINTR ) {
04013 CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
04014 xscFile.clear();
04015 xscFile.open(fname);
04016 continue;
04017 }
04018 char err_msg[257];
04019 sprintf(err_msg, "Error opening XSC restart file %s",fname);
04020 NAMD_err(err_msg);
04021 }
04022 xscFile << "# NAMD extended system configuration restart file" << std::endl;
04023 writeExtendedSystemLabels(xscFile);
04024 writeExtendedSystemData(step,xscFile);
04025 if (!xscFile) {
04026 char err_msg[257];
04027 sprintf(err_msg, "Error writing XSC restart file %s",fname);
04028 NAMD_err(err_msg);
04029 }
04030 }
04031
04032 }
04033
04034
04035 if (step == FILE_OUTPUT || step == END_OF_RUN)
04036 {
04037 int realstep = ( step == FILE_OUTPUT ?
04038 simParams->firstTimestep : simParams->N );
04039 iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
04040 << realstep << "\n" << endi;
04041 static char fname[140];
04042 strcpy(fname, simParams->outputFilename);
04043 strcat(fname, ".xsc");
04044 NAMD_backup_file(fname);
04045 ofstream_namd xscFile(fname);
04046 while (!xscFile) {
04047 if ( errno == EINTR ) {
04048 CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
04049 xscFile.clear();
04050 xscFile.open(fname);
04051 continue;
04052 }
04053 char err_msg[257];
04054 sprintf(err_msg, "Error opening XSC output file %s",fname);
04055 NAMD_err(err_msg);
04056 }
04057 xscFile << "# NAMD extended system configuration output file" << std::endl;
04058 writeExtendedSystemLabels(xscFile);
04059 writeExtendedSystemData(realstep,xscFile);
04060 if (!xscFile) {
04061 char err_msg[257];
04062 sprintf(err_msg, "Error writing XSC output file %s",fname);
04063 NAMD_err(err_msg);
04064 }
04065 }
04066
04067
04068 if (step == END_OF_RUN) {
04069 if ( xstFile.is_open() ) {
04070 xstFile.close();
04071 iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
04072 }
04073 }
04074
04075 }
04076
04077
04078 void Controller::recvCheckpointReq(const char *key, int task, checkpoint &cp) {
04079 switch ( task ) {
04080 case SCRIPT_CHECKPOINT_STORE:
04081 if ( ! checkpoints.count(key) ) {
04082 checkpoints[key] = new checkpoint;
04083 }
04084 *checkpoints[key] = cp;
04085 break;
04086 case SCRIPT_CHECKPOINT_LOAD:
04087 if ( ! checkpoints.count(key) ) {
04088 NAMD_die("Unable to load checkpoint, requested key was never stored.");
04089 }
04090 cp = *checkpoints[key];
04091 break;
04092 case SCRIPT_CHECKPOINT_SWAP:
04093 if ( ! checkpoints.count(key) ) {
04094 NAMD_die("Unable to swap checkpoint, requested key was never stored.");
04095 }
04096 std::swap(cp,*checkpoints[key]);
04097 break;
04098 case SCRIPT_CHECKPOINT_FREE:
04099 if ( ! checkpoints.count(key) ) {
04100 NAMD_die("Unable to free checkpoint, requested key was never stored.");
04101 }
04102 delete checkpoints[key];
04103 checkpoints.erase(key);
04104 break;
04105 }
04106 }
04107
04108 void Controller::recvCheckpointAck(checkpoint &cp) {
04109 if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
04110 state->lattice = cp.lattice;
04111 *(ControllerState*)this = cp.state;
04112 }
04113 CkpvAccess(_qd)->process();
04114 }
04115
04116
04117 void Controller::rebalanceLoad(int step)
04118 {
04119 if ( ! ldbSteps ) {
04120 ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
04121 }
04122 if ( ! --ldbSteps ) {
04123 startBenchTime -= CmiWallTimer();
04124 Node::Object()->outputPatchComputeMaps("before_ldb", step);
04125 LdbCoordinator::Object()->rebalance(this);
04126 startBenchTime += CmiWallTimer();
04127 fflush_count = 3;
04128 }
04129 }
04130
04131 void Controller::cycleBarrier(int doBarrier, int step) {
04132 #if USE_BARRIER
04133 if (doBarrier) {
04134 broadcast->cycleBarrier.publish(step,1);
04135 CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
04136 CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04137 }
04138 #endif
04139 }
04140
04141 void Controller::traceBarrier(int turnOnTrace, int step) {
04142 CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04143 CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04144 nd.traceBarrier(turnOnTrace, step);
04145 CthSuspend();
04146 }
04147
04148 void Controller::resumeAfterTraceBarrier(int step){
04149 broadcast->traceBarrier.publish(step,1);
04150 CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04151 awaken();
04152 }
04153
04154 #ifdef MEASURE_NAMD_WITH_PAPI
04155 void Controller::papiMeasureBarrier(int turnOnMeasure, int step){
04156 CkPrintf("Cycle time at PAPI measurement sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04157 CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04158 nd.papiMeasureBarrier(turnOnMeasure, step);
04159 CthSuspend();
04160 }
04161
04162 void Controller::resumeAfterPapiMeasureBarrier(int step){
04163 broadcast->papiMeasureBarrier.publish(step,1);
04164 CkPrintf("Cycle time at PAPI measurement sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04165 awaken();
04166 }
04167 #endif
04168
04169 void Controller::terminate(void) {
04170 BackEnd::awaken();
04171 CthFree(thread);
04172 CthSuspend();
04173 }
04174
