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 <iomanip>
00040 #include <errno.h>
00041
00042 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
00043 extern "C" void CApplicationDepositNode0Data(char *);
00044 #endif
00045
00046 #ifndef cbrt
00047
00048 #define cbrt(x) pow(x,(double)(1.0/3.0))
00049 #endif
00050
00051
00052 #define MIN_DEBUG_LEVEL 3
00053
00054 #include "Debug.h"
00055
00056 #define XXXBIGREAL 1.0e32
00057
00058 class PressureProfileReduction {
00059 private:
00060 RequireReduction *reduction;
00061 const int nslabs;
00062 const int freq;
00063 int nelements;
00064 int count;
00065 char *name;
00066
00067 BigReal *average;
00068
00069 public:
00070 PressureProfileReduction(int rtag, int numslabs, int numpartitions,
00071 const char *myname, int outputfreq)
00072 : nslabs(numslabs), freq(outputfreq) {
00073 name = strdup(myname);
00074 nelements = 3*nslabs * numpartitions;
00075 reduction = ReductionMgr::Object()->willRequire(rtag,nelements);
00076
00077 average = new BigReal[nelements];
00078 count = 0;
00079 }
00080 ~PressureProfileReduction() {
00081 delete [] average;
00082 delete reduction;
00083 free(name);
00084 }
00085
00086 void getData(int firsttimestep, int step, const Lattice &lattice,
00087 BigReal *total) {
00088 reduction->require();
00089
00090 int i;
00091 double inv_volume = 1.0 / lattice.volume();
00092
00093 int arraysize = 3*nslabs;
00094 for (i=0; i<nelements; i++) {
00095 BigReal val = reduction->item(i) * inv_volume;
00096 average[i] += val;
00097 total[i % arraysize] += val;
00098 }
00099 count++;
00100 if (!(step % freq)) {
00101
00102 BigReal scalefac = PRESSUREFACTOR * nslabs;
00103
00104 iout << "PPROFILE" << name << ": " << step << " ";
00105 if (step == firsttimestep) {
00106
00107 for (i=0; i<nelements; i++)
00108 iout << reduction->item(i)*scalefac*inv_volume << " ";
00109 } else {
00110
00111 scalefac /= count;
00112 for (i=0; i<nelements; i++)
00113 iout << average[i]*scalefac << " ";
00114 }
00115 iout << "\n" << endi;
00116
00117 memset(average, 0, nelements*sizeof(BigReal));
00118 count = 0;
00119 }
00120 }
00121 };
00122
00123
00124 Controller::Controller(NamdState *s) :
00125 computeChecksum(0), marginViolations(0), pairlistWarnings(0),
00126 simParams(Node::Object()->simParameters),
00127 state(s),
00128 collection(CollectionMaster::Object()),
00129 startCTime(0),
00130 firstCTime(CmiTimer()),
00131 startWTime(0),
00132 firstWTime(CmiWallTimer()),
00133 startBenchTime(0),
00134 ldbSteps(0),
00135 fflush_count(3)
00136 {
00137 broadcast = new ControllerBroadcasts;
00138 reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
00139
00140 if (simParams->accelMDOn) {
00141 amd_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_AMD);
00142 } else {
00143 amd_reduction = NULL;
00144 }
00145
00146 pressureProfileSlabs = 0;
00147 pressureProfileCount = 0;
00148 ppbonded = ppnonbonded = ppint = NULL;
00149 if (simParams->pressureProfileOn || simParams->pressureProfileEwaldOn) {
00150 int ntypes = simParams->pressureProfileAtomTypes;
00151 int nslabs = pressureProfileSlabs = simParams->pressureProfileSlabs;
00152
00153 int npairs = (ntypes * (ntypes+1))/2;
00154 pressureProfileAverage = new BigReal[3*nslabs];
00155 memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
00156 int freq = simParams->pressureProfileFreq;
00157 if (simParams->pressureProfileOn) {
00158 ppbonded = new PressureProfileReduction(REDUCTIONS_PPROF_BONDED,
00159 nslabs, npairs, "BONDED", freq);
00160 ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED,
00161 nslabs, npairs, "NONBONDED", freq);
00162
00163 ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL,
00164 nslabs, ntypes, "INTERNAL", freq);
00165 } else {
00166
00167 ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED,
00168 nslabs, npairs, "NONBONDED", freq);
00169 }
00170 }
00171 random = new Random(simParams->randomSeed);
00172 random->split(0,PatchMap::Object()->numPatches()+1);
00173
00174 rescaleVelocities_sumTemps = 0; rescaleVelocities_numTemps = 0;
00175 berendsenPressure_avg = 0; berendsenPressure_count = 0;
00176
00177 langevinPiston_strainRate =
00178 Tensor::symmetric(simParams->strainRate,simParams->strainRate2);
00179 if ( ! simParams->useFlexibleCell ) {
00180 BigReal avg = trace(langevinPiston_strainRate) / 3.;
00181 langevinPiston_strainRate = Tensor::identity(avg);
00182 } else if ( simParams->useConstantRatio ) {
00183 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
00184 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
00185 AVGXY(langevinPiston_strainRate);
00186 #undef AVGXY
00187 }
00188 smooth2_avg = XXXBIGREAL;
00189 temp_avg = 0;
00190 pressure_avg = 0;
00191 groupPressure_avg = 0;
00192 avg_count = 0;
00193 pressure_tavg = 0;
00194 groupPressure_tavg = 0;
00195 tavg_count = 0;
00196 checkpoint_stored = 0;
00197 drudeComTemp = 0;
00198 drudeBondTemp = 0;
00199 drudeComTempAvg = 0;
00200 drudeBondTempAvg = 0;
00201 }
00202
00203 Controller::~Controller(void)
00204 {
00205 delete broadcast;
00206 delete reduction;
00207 delete amd_reduction;
00208 delete ppbonded;
00209 delete ppnonbonded;
00210 delete ppint;
00211 delete [] pressureProfileAverage;
00212 delete random;
00213 }
00214
00215 void Controller::threadRun(Controller* arg)
00216 {
00217 arg->algorithm();
00218 }
00219
00220 void Controller::run(void)
00221 {
00222
00223 DebugM(4, "Starting thread in controller on this=" << this << "\n");
00224 thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
00225 CthSetStrategyDefault(thread);
00226 #if CMK_BLUEGENE_CHARM
00227 BgAttach(thread);
00228 #endif
00229 awaken();
00230 }
00231
00232
00233 void Controller::algorithm(void)
00234 {
00235 int scriptTask;
00236 int scriptSeq = 0;
00237 BackEnd::awaken();
00238 while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
00239 switch ( scriptTask ) {
00240 case SCRIPT_OUTPUT:
00241 enqueueCollections(FILE_OUTPUT);
00242 outputExtendedSystem(FILE_OUTPUT);
00243 break;
00244 case SCRIPT_FORCEOUTPUT:
00245 enqueueCollections(FORCE_OUTPUT);
00246 break;
00247 case SCRIPT_MEASURE:
00248 enqueueCollections(EVAL_MEASURE);
00249 break;
00250 case SCRIPT_REINITVELS:
00251 iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
00252 << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
00253 break;
00254 case SCRIPT_RESCALEVELS:
00255 iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
00256 << " BY " << simParams->scriptArg1 << "\n" << endi;
00257 break;
00258 case SCRIPT_CHECKPOINT:
00259 iout << "CHECKPOINTING POSITIONS AT STEP " << simParams->firstTimestep
00260 << "\n" << endi;
00261 checkpoint_stored = 1;
00262 checkpoint_lattice = state->lattice;
00263 checkpoint_langevinPiston_strainRate = langevinPiston_strainRate;
00264 checkpoint_berendsenPressure_avg = berendsenPressure_avg;
00265 checkpoint_berendsenPressure_count = berendsenPressure_count;
00266 checkpoint_smooth2_avg = smooth2_avg;
00267 break;
00268 case SCRIPT_REVERT:
00269 iout << "REVERTING POSITIONS AT STEP " << simParams->firstTimestep
00270 << "\n" << endi;
00271 if ( ! checkpoint_stored )
00272 NAMD_die("Unable to revert, checkpoint was never called!");
00273 state->lattice = checkpoint_lattice;
00274 langevinPiston_strainRate = checkpoint_langevinPiston_strainRate;
00275 berendsenPressure_avg = checkpoint_berendsenPressure_avg;
00276 berendsenPressure_count = checkpoint_berendsenPressure_count;
00277 smooth2_avg = checkpoint_smooth2_avg;
00278 break;
00279 case SCRIPT_MINIMIZE:
00280 minimize();
00281 break;
00282 case SCRIPT_RUN:
00283 integrate();
00284 break;
00285 }
00286 BackEnd::awaken();
00287 }
00288 enqueueCollections(END_OF_RUN);
00289 outputExtendedSystem(END_OF_RUN);
00290
00291
00292 if(traceAvailable())
00293 traceClose();
00294 terminate();
00295 }
00296
00297
00298
00299
00300
00301
00302 extern int eventEndOfTimeStep;
00303
00304
00305 static int gotsigint = 0;
00306 static void my_sigint_handler(int sig) {
00307 if (sig == SIGINT) gotsigint = 1;
00308 }
00309 extern "C" {
00310 typedef void (*namd_sighandler_t)(int);
00311 }
00312
00313 void Controller::integrate() {
00314 char traceNote[24];
00315
00316 int step = simParams->firstTimestep;
00317
00318 const int numberOfSteps = simParams->N;
00319 const int stepsPerCycle = simParams->stepsPerCycle;
00320
00321 const int zeroMomentum = simParams->zeroMomentum;
00322
00323 nbondFreq = simParams->nonbondedFrequency;
00324 const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00325 if (dofull)
00326 slowFreq = simParams->fullElectFrequency;
00327 else
00328 slowFreq = simParams->nonbondedFrequency;
00329 if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
00330
00331 reassignVelocities(step);
00332 rescaleaccelMD(step);
00333 adaptTempUpdate(step);
00334
00335 receivePressure(step);
00336 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00337 correctMomentum(step);
00338 printFepMessage(step);
00339 printTiMessage(step);
00340 printDynamicsEnergies(step);
00341 outputFepEnergy(step);
00342 outputTiEnergy(step);
00343 if(traceIsOn()){
00344 traceUserEvent(eventEndOfTimeStep);
00345 sprintf(traceNote, "s:%d", step);
00346 traceUserSuppliedNote(traceNote);
00347 }
00348 outputExtendedSystem(step);
00349 rebalanceLoad(step);
00350
00351
00352
00353
00354
00355
00356 for ( ++step ; step <= numberOfSteps; ++step )
00357 {
00358 adaptTempUpdate(step);
00359 rescaleVelocities(step);
00360 tcoupleVelocities(step);
00361 berendsenPressure(step);
00362 langevinPiston1(step);
00363 rescaleaccelMD(step);
00364 enqueueCollections(step);
00365 receivePressure(step);
00366 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00367 correctMomentum(step);
00368 langevinPiston2(step);
00369 reassignVelocities(step);
00370 printDynamicsEnergies(step);
00371 outputFepEnergy(step);
00372 outputTiEnergy(step);
00373 if(traceIsOn()){
00374 traceUserEvent(eventEndOfTimeStep);
00375 sprintf(traceNote, "s:%d", step);
00376 traceUserSuppliedNote(traceNote);
00377 }
00378
00379
00380
00381
00382 outputExtendedSystem(step);
00383 #if CYCLE_BARRIER
00384 cycleBarrier(!((step+1) % stepsPerCycle),step);
00385 #elif PME_BARRIER
00386 cycleBarrier(dofull && !(step%slowFreq),step);
00387 #elif STEP_BARRIER
00388 cycleBarrier(1, step);
00389 #endif
00390
00391 if(Node::Object()->specialTracing || simParams->statsOn){
00392 int bstep = simParams->traceStartStep;
00393 int estep = bstep + simParams->numTraceSteps;
00394 if(step == bstep){
00395 traceBarrier(1, step);
00396 }
00397 if(step == estep){
00398 traceBarrier(0, step);
00399 }
00400 }
00401
00402 #ifdef MEASURE_NAMD_WITH_PAPI
00403 if(simParams->papiMeasure) {
00404 int bstep = simParams->papiMeasureStartStep;
00405 int estep = bstep + simParams->numPapiMeasureSteps;
00406 if(step == bstep) {
00407 papiMeasureBarrier(1, step);
00408 }
00409 if(step == estep) {
00410 papiMeasureBarrier(0, step);
00411 }
00412 }
00413 #endif
00414
00415 rebalanceLoad(step);
00416
00417 #if PME_BARRIER
00418 cycleBarrier(dofull && !((step+1)%slowFreq),step);
00419 #endif
00420 }
00421
00422 }
00423
00424
00425 #define CALCULATE \
00426 printMinimizeEnergies(step); \
00427 outputExtendedSystem(step); \
00428 rebalanceLoad(step); \
00429 if ( step == numberOfSteps ) return; \
00430 else ++step;
00431
00432 #define MOVETO(X) \
00433 if ( step == numberOfSteps ) { \
00434 if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \
00435 if ( newDir || (mid.x-last.x) ) { \
00436 broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \
00437 } else { \
00438 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00439 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00440 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00441 } \
00442 enqueueCollections(step); \
00443 CALCULATE \
00444 } else if ( (X)-last.x ) { \
00445 broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \
00446 newDir = 0; \
00447 last.x = (X); \
00448 enqueueCollections(step); \
00449 CALCULATE \
00450 last.u = min_energy; \
00451 last.dudx = -1. * min_f_dot_v; \
00452 last.noGradient = min_huge_count; \
00453 if ( minVerbose ) { \
00454 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \
00455 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \
00456 iout << "\n" << endi; \
00457 } \
00458 }
00459
00460 struct minpoint {
00461 BigReal x, u, dudx, noGradient;
00462 minpoint() : x(0), u(0), dudx(0), noGradient(1) { ; }
00463 };
00464
00465 void Controller::minimize() {
00466
00467
00468 const int minVerbose = simParams->minVerbose;
00469 const int numberOfSteps = simParams->N;
00470 int step = simParams->firstTimestep;
00471 slowFreq = nbondFreq = 1;
00472 BigReal tinystep = simParams->minTinyStep;
00473 BigReal babystep = simParams->minBabyStep;
00474 BigReal linegoal = simParams->minLineGoal;
00475 BigReal initstep = tinystep;
00476 const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
00477
00478 CALCULATE
00479
00480 int minSeq = 0;
00481
00482
00483 int old_min_huge_count = 2000000000;
00484 int steps_at_same_min_huge_count = 0;
00485 for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
00486 if ( min_huge_count >= old_min_huge_count ) {
00487 if ( ++steps_at_same_min_huge_count > 10 ) break;
00488 } else {
00489 old_min_huge_count = min_huge_count;
00490 steps_at_same_min_huge_count = 0;
00491 }
00492 iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
00493 broadcast->minimizeCoefficient.publish(minSeq++,1.);
00494 enqueueCollections(step);
00495 CALCULATE
00496 }
00497 if ( min_huge_count ) {
00498 iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
00499 }
00500 iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00501
00502 int atStart = 2;
00503 int errorFactor = 10;
00504 BigReal old_f_dot_f = min_f_dot_f;
00505 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00506 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00507 int newDir = 1;
00508 min_f_dot_v = min_f_dot_f;
00509 min_v_dot_v = min_f_dot_f;
00510 while ( 1 ) {
00511
00512
00513 minpoint lo,hi,mid,last;
00514 BigReal x = 0;
00515 lo.x = x;
00516 lo.u = min_energy;
00517 lo.dudx = -1. * min_f_dot_v;
00518 lo.noGradient = min_huge_count;
00519 mid = lo;
00520 last = mid;
00521 if ( minVerbose ) {
00522 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
00523 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
00524 iout << "\n" << endi;
00525 }
00526 BigReal tol = fabs( linegoal * min_f_dot_v );
00527 if ( initstep > babystep ) initstep = babystep;
00528 if ( initstep < 1.0e-300 ) initstep = 1.0e-300;
00529 iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
00530 fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
00531 int start_with_huge = last.noGradient;
00532 x = initstep;
00533 x *= sqrt( min_f_dot_f / min_v_dot_v ); MOVETO(x)
00534 if ( ! start_with_huge ) {
00535 for ( int i_huge = 0; last.noGradient && i_huge < 10; ++i_huge ) {
00536 x *= 0.25; MOVETO(x);
00537 initstep *= 0.25;
00538 }
00539 }
00540
00541 initstep *= 0.25;
00542 BigReal maxinitstep = initstep * 16.0;
00543 while ( last.u < mid.u ) {
00544 initstep *= 2.0;
00545 lo = mid; mid = last;
00546 x *= 2.0; MOVETO(x)
00547 }
00548 if ( initstep > maxinitstep ) initstep = maxinitstep;
00549 hi = last;
00550 #define PRINT_BRACKET \
00551 iout << "LINE MINIMIZER BRACKET: DX " \
00552 << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
00553 " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
00554 lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
00555 PRINT_BRACKET
00556
00557 int itcnt;
00558 for ( itcnt = 10; itcnt > 0; --itcnt ) {
00559 int progress = 1;
00560
00561 if ( mid.noGradient ) {
00562 if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) {
00563 x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
00564 MOVETO(x)
00565 if ( last.u <= mid.u ) {
00566 hi = mid; mid = last;
00567 } else {
00568 lo = last;
00569 }
00570 } else {
00571 x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
00572 MOVETO(x)
00573 if ( last.u <= mid.u ) {
00574 lo = mid; mid = last;
00575 } else {
00576 hi = last;
00577 }
00578 }
00579 } else {
00580 if ( mid.dudx > 0. ) {
00581 BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
00582 BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
00583 x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
00584 BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
00585 if ( xdiv ) x /= xdiv; else x = last.x;
00586 if ( x > altxhi ) x = altxhi;
00587 if ( x < altxlo ) x = altxlo;
00588 if ( x-last.x == 0 ) break;
00589 MOVETO(x)
00590 if ( last.u <= mid.u ) {
00591 hi = mid; mid = last;
00592 } else {
00593 if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
00594 lo = last;
00595 }
00596 } else {
00597 BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
00598 BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
00599 x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
00600 BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
00601 if ( xdiv ) x /= xdiv; else x = last.x;
00602 if ( x < altxlo ) x = altxlo;
00603 if ( x > altxhi ) x = altxhi;
00604 if ( x-last.x == 0 ) break;
00605 MOVETO(x)
00606 if ( last.u <= mid.u ) {
00607 lo = mid; mid = last;
00608 } else {
00609 if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
00610 hi = last;
00611 }
00612 }
00613 }
00614 PRINT_BRACKET
00615 if ( ! progress ) {
00616 MOVETO(mid.x);
00617 break;
00618 }
00619 if ( fabs(last.dudx) < tol ) break;
00620 if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
00621 if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
00622 }
00623
00624 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00625 BigReal c = min_f_dot_f / old_f_dot_f;
00626 c = ( c > 1.5 ? 1.5 : c );
00627 if ( atStart ) { c = 0; --atStart; }
00628 if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
00629 c = 0;
00630 if ( errorFactor < 100 ) errorFactor += 10;
00631 }
00632 if ( c == 0 ) {
00633 iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00634 }
00635 broadcast->minimizeCoefficient.publish(minSeq++,c);
00636 newDir = 1;
00637 old_f_dot_f = min_f_dot_f;
00638 min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
00639 min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
00640 }
00641
00642 }
00643
00644 #undef MOVETO
00645 #undef CALCULATE
00646
00647 void Controller::berendsenPressure(int step)
00648 {
00649 if ( simParams->berendsenPressureOn ) {
00650 berendsenPressure_count += 1;
00651 berendsenPressure_avg += controlPressure;
00652 const int freq = simParams->berendsenPressureFreq;
00653 if ( ! (berendsenPressure_count % freq) ) {
00654 Tensor factor = berendsenPressure_avg / berendsenPressure_count;
00655 berendsenPressure_avg = 0;
00656 berendsenPressure_count = 0;
00657
00658 factor = Tensor::diagonal(diagonal(factor));
00659 factor -= Tensor::identity(simParams->berendsenPressureTarget);
00660 factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
00661 simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
00662 factor += Tensor::identity(1.0);
00663 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
00664 if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
00665 if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
00666 int limited = 0;
00667 LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
00668 LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
00669 LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
00670 #undef LIMIT_SCALING
00671 if ( limited ) {
00672 iout << iERROR << "Step " << step <<
00673 " cell rescaling factor limited.\n" << endi;
00674 }
00675 broadcast->positionRescaleFactor.publish(step,factor);
00676 state->lattice.rescale(factor);
00677 }
00678 } else {
00679 berendsenPressure_avg = 0;
00680 berendsenPressure_count = 0;
00681 }
00682 }
00683
00684 void Controller::langevinPiston1(int step)
00685 {
00686 if ( simParams->langevinPistonOn )
00687 {
00688 Tensor &strainRate = langevinPiston_strainRate;
00689 int cellDims = simParams->useFlexibleCell ? 1 : 3;
00690 BigReal dt = simParams->dt;
00691 BigReal dt_long = slowFreq * dt;
00692 BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
00693 BigReal tau = simParams->langevinPistonPeriod;
00694 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00695
00696 #ifdef DEBUG_PRESSURE
00697 iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
00698 #endif
00699
00700 if ( ! ( (step-1) % slowFreq ) )
00701 {
00702 BigReal gamma = 1 / simParams->langevinPistonDecay;
00703 BigReal f1 = exp( -0.5 * dt_long * gamma );
00704 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00705 strainRate *= f1;
00706 if ( simParams->useFlexibleCell ) {
00707
00708 if ( simParams->useConstantRatio ) {
00709 BigReal r = f2 * random->gaussian();
00710 strainRate.xx += r;
00711 strainRate.yy += r;
00712 strainRate.zz += f2 * random->gaussian();
00713 } else {
00714 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00715 }
00716 } else {
00717 strainRate += f2 * Tensor::identity(random->gaussian());
00718 }
00719
00720 #ifdef DEBUG_PRESSURE
00721 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00722 #endif
00723 }
00724
00725
00726
00727
00728 Tensor ptarget;
00729 ptarget.zz = simParams->langevinPistonTarget;
00730 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
00731 simParams->surfaceTensionTarget / state->lattice.c().z;
00732
00733 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00734 ( controlPressure - ptarget );
00735
00736 if (simParams->fixCellDims) {
00737 if (simParams->fixCellDimX) strainRate.xx = 0;
00738 if (simParams->fixCellDimY) strainRate.yy = 0;
00739 if (simParams->fixCellDimZ) strainRate.zz = 0;
00740 }
00741
00742
00743 #ifdef DEBUG_PRESSURE
00744 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00745 #endif
00746
00747 if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
00748 {
00749
00750 Tensor factor;
00751 if ( !simParams->useConstantArea ) {
00752 factor.xx = exp( dt_long * strainRate.xx );
00753 factor.yy = exp( dt_long * strainRate.yy );
00754 } else {
00755 factor.xx = factor.yy = 1;
00756 }
00757 factor.zz = exp( dt_long * strainRate.zz );
00758 broadcast->positionRescaleFactor.publish(step,factor);
00759 state->lattice.rescale(factor);
00760 #ifdef DEBUG_PRESSURE
00761 iout << iINFO << "rescaling by: " << factor << "\n";
00762 #endif
00763 }
00764
00765
00766 if ( ! ( step % nbondFreq ) )
00767 {
00768 #ifdef DEBUG_PRESSURE
00769 iout << iINFO << "correcting strain rate for nbond, ";
00770 #endif
00771 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00772 ( (nbondFreq - 1) * controlPressure_nbond );
00773 #ifdef DEBUG_PRESSURE
00774 iout << "strain rate: " << strainRate << "\n";
00775 #endif
00776 }
00777 if ( ! ( step % slowFreq ) )
00778 {
00779 #ifdef DEBUG_PRESSURE
00780 iout << iINFO << "correcting strain rate for slow, ";
00781 #endif
00782 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00783 ( (slowFreq - 1) * controlPressure_slow );
00784 #ifdef DEBUG_PRESSURE
00785 iout << "strain rate: " << strainRate << "\n";
00786 #endif
00787 }
00788 if (simParams->fixCellDims) {
00789 if (simParams->fixCellDimX) strainRate.xx = 0;
00790 if (simParams->fixCellDimY) strainRate.yy = 0;
00791 if (simParams->fixCellDimZ) strainRate.zz = 0;
00792 }
00793
00794 }
00795
00796 }
00797
00798 void Controller::langevinPiston2(int step)
00799 {
00800 if ( simParams->langevinPistonOn )
00801 {
00802 Tensor &strainRate = langevinPiston_strainRate;
00803 int cellDims = simParams->useFlexibleCell ? 1 : 3;
00804 BigReal dt = simParams->dt;
00805 BigReal dt_long = slowFreq * dt;
00806 BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
00807 BigReal tau = simParams->langevinPistonPeriod;
00808 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00809
00810
00811 if ( ! ( step % nbondFreq ) )
00812 {
00813 #ifdef DEBUG_PRESSURE
00814 iout << iINFO << "correcting strain rate for nbond, ";
00815 #endif
00816 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00817 ( (nbondFreq - 1) * controlPressure_nbond );
00818 #ifdef DEBUG_PRESSURE
00819 iout << "strain rate: " << strainRate << "\n";
00820 #endif
00821 }
00822 if ( ! ( step % slowFreq ) )
00823 {
00824 #ifdef DEBUG_PRESSURE
00825 iout << iINFO << "correcting strain rate for slow, ";
00826 #endif
00827 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00828 ( (slowFreq - 1) * controlPressure_slow );
00829 #ifdef DEBUG_PRESSURE
00830 iout << "strain rate: " << strainRate << "\n";
00831 #endif
00832 }
00833
00834
00835
00836
00837 Tensor ptarget;
00838 ptarget.zz = simParams->langevinPistonTarget;
00839 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
00840 simParams->surfaceTensionTarget / state->lattice.c().z;
00841
00842 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00843 ( controlPressure - ptarget );
00844
00845
00846 #ifdef DEBUG_PRESSURE
00847 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00848 #endif
00849
00850 if ( ! ( step % slowFreq ) )
00851 {
00852 BigReal gamma = 1 / simParams->langevinPistonDecay;
00853 BigReal f1 = exp( -0.5 * dt_long * gamma );
00854 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00855 strainRate *= f1;
00856 if ( simParams->useFlexibleCell ) {
00857
00858 if ( simParams->useConstantRatio ) {
00859 BigReal r = f2 * random->gaussian();
00860 strainRate.xx += r;
00861 strainRate.yy += r;
00862 strainRate.zz += f2 * random->gaussian();
00863 } else {
00864 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00865 }
00866 } else {
00867 strainRate += f2 * Tensor::identity(random->gaussian());
00868 }
00869 #ifdef DEBUG_PRESSURE
00870 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00871 #endif
00872 }
00873
00874 #ifdef DEBUG_PRESSURE
00875 iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
00876 #endif
00877 if (simParams->fixCellDims) {
00878 if (simParams->fixCellDimX) strainRate.xx = 0;
00879 if (simParams->fixCellDimY) strainRate.yy = 0;
00880 if (simParams->fixCellDimZ) strainRate.zz = 0;
00881 }
00882 }
00883
00884
00885 }
00886
00887 void Controller::rescaleVelocities(int step)
00888 {
00889 const int rescaleFreq = simParams->rescaleFreq;
00890 if ( rescaleFreq > 0 ) {
00891 rescaleVelocities_sumTemps += temperature; ++rescaleVelocities_numTemps;
00892 if ( rescaleVelocities_numTemps == rescaleFreq ) {
00893 BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
00894 BigReal rescaleTemp = simParams->rescaleTemp;
00895 if ( simParams->adaptTempOn && simParams->adaptTempRescale && (step > simParams->adaptTempFirstStep ) &&
00896 (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
00897 rescaleTemp = adaptTempT;
00898 }
00899 BigReal factor = sqrt(rescaleTemp/avgTemp);
00900 broadcast->velocityRescaleFactor.publish(step,factor);
00901
00902
00903
00904 rescaleVelocities_sumTemps = 0; rescaleVelocities_numTemps = 0;
00905 }
00906 }
00907 }
00908
00909 void Controller::correctMomentum(int step) {
00910
00911 Vector momentum;
00912 momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
00913 momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
00914 momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
00915 const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
00916
00917 if ( momentum.length2() == 0. )
00918 iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
00919 if ( mass == 0. )
00920 NAMD_die("Total mass is zero in Controller::correctMomentum");
00921
00922 momentum *= (-1./mass);
00923
00924 broadcast->momentumCorrection.publish(step+slowFreq,momentum);
00925 }
00926
00927
00928 void Controller::printFepMessage(int step)
00929 {
00930 if (simParams->alchFepOn) {
00931 const BigReal alchLambda = simParams->alchLambda;
00932 const BigReal alchLambda2 = simParams->alchLambda2;
00933 const BigReal alchTemp = simParams->alchTemp;
00934 const int alchEquilSteps = simParams->alchEquilSteps;
00935 iout << "FEP: RESETTING FOR NEW FEP WINDOW "
00936 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2
00937 << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
00938 << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
00939 << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
00940 << " K FOR FEP CALCULATION\n" << endi;
00941 }
00942 }
00943 void Controller::printTiMessage(int step)
00944 {
00945 if (simParams->alchThermIntOn) {
00946 const BigReal alchLambda = simParams->alchLambda;
00947 const BigReal alchTemp = simParams->alchTemp;
00948 const int alchEquilSteps = simParams->alchEquilSteps;
00949 iout << "TI: RESETTING FOR NEW WINDOW "
00950 << "LAMBDA SET TO " << alchLambda
00951 << "\nTI: WINDOW TO HAVE " << alchEquilSteps
00952 << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
00953 }
00954 }
00955
00956
00957 void Controller::reassignVelocities(int step)
00958 {
00959 const int reassignFreq = simParams->reassignFreq;
00960 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
00961 BigReal newTemp = simParams->reassignTemp;
00962 newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
00963 if ( simParams->reassignIncr > 0.0 ) {
00964 if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
00965 newTemp = simParams->reassignHold;
00966 } else {
00967 if ( newTemp < simParams->reassignHold )
00968 newTemp = simParams->reassignHold;
00969 }
00970 iout << "REASSIGNING VELOCITIES AT STEP " << step
00971 << " TO " << newTemp << " KELVIN.\n" << endi;
00972 }
00973 }
00974
00975 void Controller::tcoupleVelocities(int step)
00976 {
00977 if ( simParams->tCoupleOn )
00978 {
00979 const BigReal tCoupleTemp = simParams->tCoupleTemp;
00980 BigReal coefficient = 1.;
00981 if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
00982 broadcast->tcoupleCoefficient.publish(step,coefficient);
00983 }
00984 }
00985
00986 static char *FORMAT(BigReal X)
00987 {
00988 static char tmp_string[25];
00989 const double maxnum = 9999999999.9999;
00990 if ( X > maxnum ) X = maxnum;
00991 if ( X < -maxnum ) X = -maxnum;
00992 sprintf(tmp_string," %14.4f",X);
00993 return tmp_string;
00994 }
00995
00996 static char *FORMAT(const char *X)
00997 {
00998 static char tmp_string[25];
00999 sprintf(tmp_string," %14s",X);
01000 return tmp_string;
01001 }
01002
01003 static char *ETITLE(int X)
01004 {
01005 static char tmp_string[21];
01006 sprintf(tmp_string,"ENERGY: %7d",X);
01007 return tmp_string;
01008 }
01009
01010 void Controller::receivePressure(int step, int minimize)
01011 {
01012 Node *node = Node::Object();
01013 Molecule *molecule = node->molecule;
01014 SimParameters *simParameters = node->simParameters;
01015 Lattice &lattice = state->lattice;
01016
01017 reduction->require();
01018
01019 Tensor virial;
01020 Tensor virial_normal;
01021 Tensor virial_nbond;
01022 Tensor virial_slow;
01023 Tensor pressure_amd;
01024 #ifdef ALTVIRIAL
01025 Tensor altVirial_normal;
01026 Tensor altVirial_nbond;
01027 Tensor altVirial_slow;
01028 #endif
01029 Tensor intVirial;
01030 Tensor intVirial_normal;
01031 Tensor intVirial_nbond;
01032 Tensor intVirial_slow;
01033 Vector extForce_normal;
01034 Vector extForce_nbond;
01035 Vector extForce_slow;
01036 BigReal volume;
01037
01038 #if 1
01039 numDegFreedom = molecule->num_deg_freedom();
01040 int numGroupDegFreedom = molecule->num_group_deg_freedom();
01041 int numFixedGroups = molecule->num_fixed_groups();
01042 int numFixedAtoms = molecule->num_fixed_atoms();
01043 #endif
01044 #if 0
01045 int numAtoms = molecule->numAtoms;
01046 numDegFreedom = 3 * numAtoms;
01047 int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
01048 int numFixedAtoms =
01049 ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
01050 int numLonepairs = molecule->numLonepairs;
01051 int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
01052 if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
01053 if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01054 if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01055 if ( ! ( numFixedAtoms || molecule->numConstraints
01056 || simParameters->comMove || simParameters->langevinOn ) ) {
01057 numDegFreedom -= 3;
01058 numGroupDegFreedom -= 3;
01059 }
01060 if (simParameters->pairInteractionOn) {
01061
01062 numDegFreedom = 3 * molecule->numFepInitial;
01063 }
01064 int numRigidBonds = molecule->numRigidBonds;
01065 int numFixedRigidBonds =
01066 ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01067
01068
01069
01070 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01071 #endif
01072
01073 kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01074 kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01075
01076 BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01077 reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01078 BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01079 reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01080
01081 BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01082 / ( numDegFreedom * BOLTZMANN );
01083 BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01084 / ( numDegFreedom * BOLTZMANN );
01085 BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01086 / ( numGroupDegFreedom * BOLTZMANN );
01087 BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01088 / ( numGroupDegFreedom * BOLTZMANN );
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100 GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01101 GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01102 GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01103
01104 #ifdef ALTVIRIAL
01105 GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01106 GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01107 GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01108 #endif
01109
01110 GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01111 GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01112 GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01113
01114 GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01115 GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01116 GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01117 Vector extPosition = lattice.origin();
01118 virial_normal -= outer(extForce_normal,extPosition);
01119 virial_nbond -= outer(extForce_nbond,extPosition);
01120 virial_slow -= outer(extForce_slow,extPosition);
01121
01122
01123 kineticEnergy = kineticEnergyCentered;
01124 temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );
01125
01126 if (simParameters->drudeOn) {
01127 BigReal drudeComKE
01128 = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01129 BigReal drudeBondKE
01130 = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01131 int g_bond = 3 * molecule->numDrudeAtoms;
01132 int g_com = numDegFreedom - g_bond;
01133 drudeComTemp = 2.0 * drudeComKE / (g_com * BOLTZMANN);
01134 drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
01135 }
01136
01137 if ( (volume=lattice.volume()) != 0. )
01138 {
01139
01140 if (simParameters->LJcorrection && volume) {
01141
01142
01143
01144 virial_normal += Tensor::identity(molecule->tail_corr_virial / volume);
01145 }
01146
01147
01148 pressure_normal = virial_normal / volume;
01149 groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01150
01151 if (simParameters->accelMDOn) {
01152 pressure_amd = virial_amd / volume;
01153 pressure_normal += pressure_amd;
01154 groupPressure_normal += pressure_amd;
01155 }
01156
01157 if ( minimize || ! ( step % nbondFreq ) )
01158 {
01159 pressure_nbond = virial_nbond / volume;
01160 groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01161 }
01162
01163 if ( minimize || ! ( step % slowFreq ) )
01164 {
01165 pressure_slow = virial_slow / volume;
01166 groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01167 }
01168
01169
01170
01171
01172
01173
01174
01175
01176 pressure = pressure_normal + pressure_nbond + pressure_slow;
01177 groupPressure = groupPressure_normal + groupPressure_nbond +
01178 groupPressure_slow;
01179 }
01180 else
01181 {
01182 pressure = Tensor();
01183 groupPressure = Tensor();
01184 }
01185
01186 if ( simParameters->useGroupPressure )
01187 {
01188 controlPressure_normal = groupPressure_normal;
01189 controlPressure_nbond = groupPressure_nbond;
01190 controlPressure_slow = groupPressure_slow;
01191 controlPressure = groupPressure;
01192 controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01193 if ( ! ( numFixedAtoms || molecule->numConstraints
01194 || simParameters->comMove || simParameters->langevinOn ) ) {
01195 controlNumDegFreedom -= 1;
01196 }
01197 }
01198 else
01199 {
01200 controlPressure_normal = pressure_normal;
01201 controlPressure_nbond = pressure_nbond;
01202 controlPressure_slow = pressure_slow;
01203 controlPressure = pressure;
01204 controlNumDegFreedom = numDegFreedom / 3;
01205 }
01206
01207 if (simParameters->fixCellDims) {
01208 if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01209 if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01210 if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01211 }
01212
01213 if ( simParameters->useFlexibleCell ) {
01214
01215
01216
01217
01218
01219
01220 controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01221 controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01222 controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01223 controlPressure = Tensor::diagonal(diagonal(controlPressure));
01224 if ( simParameters->useConstantRatio ) {
01225 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01226 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01227 AVGXY(controlPressure_normal);
01228 AVGXY(controlPressure_nbond);
01229 AVGXY(controlPressure_slow);
01230 AVGXY(controlPressure);
01231 #undef AVGXY
01232 }
01233 } else {
01234 controlPressure_normal =
01235 Tensor::identity(trace(controlPressure_normal)/3.);
01236 controlPressure_nbond =
01237 Tensor::identity(trace(controlPressure_nbond)/3.);
01238 controlPressure_slow =
01239 Tensor::identity(trace(controlPressure_slow)/3.);
01240 controlPressure =
01241 Tensor::identity(trace(controlPressure)/3.);
01242 }
01243
01244 #ifdef DEBUG_PRESSURE
01245 iout << iINFO << "Control pressure = " << controlPressure <<
01246 " = " << controlPressure_normal << " + " <<
01247 controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01248 #endif
01249 if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
01250 iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
01251 << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
01252 << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
01253 << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
01254 << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
01255 << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
01256 << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
01257 << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
01258 << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
01259 << endi;
01260 }
01261 }
01262
01263 void Controller::rescaleaccelMD(int step, int minimize)
01264 {
01265 if ( !simParams->accelMDOn ) return;
01266
01267 amd_reduction->require();
01268
01269 if (step == simParams->firstTimestep) accelMDdVAverage = 0;
01270
01271 if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
01272
01273 Node *node = Node::Object();
01274 Molecule *molecule = node->molecule;
01275 Lattice &lattice = state->lattice;
01276
01277 const BigReal accelMDE = simParams->accelMDE;
01278 const BigReal accelMDalpha = simParams->accelMDalpha;
01279 const BigReal accelMDTE = simParams->accelMDTE;
01280 const BigReal accelMDTalpha = simParams->accelMDTalpha;
01281 const int accelMDOutFreq = simParams->accelMDOutFreq;
01282
01283 BigReal bondEnergy;
01284 BigReal angleEnergy;
01285 BigReal dihedralEnergy;
01286 BigReal improperEnergy;
01287 BigReal crosstermEnergy;
01288 BigReal amd_electEnergy;
01289 BigReal amd_ljEnergy;
01290 BigReal amd_electEnergySlow;
01291 BigReal potentialEnergy;
01292 BigReal factor_dihe = 1;
01293 BigReal factor_tot = 1;
01294 BigReal testV;
01295 BigReal dV = 0;
01296 Vector accelMDfactor;
01297 Tensor vir;
01298 Tensor vir_dihe;
01299 Tensor vir_normal;
01300 Tensor vir_nbond;
01301 Tensor vir_slow;
01302
01303 bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
01304 angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
01305 dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
01306 improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
01307 crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
01308
01309 GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
01310 GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
01311 GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
01312 GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
01313
01314 if ( !( step % nbondFreq ) ) {
01315 amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
01316 amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
01317 } else {
01318 amd_electEnergy = electEnergy;
01319 amd_ljEnergy = ljEnergy;
01320 }
01321
01322 if ( !( step % slowFreq ) ) {
01323 amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
01324 } else {
01325 amd_electEnergySlow = electEnergySlow;
01326 }
01327
01328 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy + improperEnergy
01329 + crosstermEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy;
01330
01331 if (simParams->accelMDdihe) {
01332
01333 testV = dihedralEnergy + crosstermEnergy;
01334 if ( testV < accelMDE ) {
01335 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
01336 factor_dihe *= factor_dihe;
01337 vir = vir_dihe * (factor_dihe - 1.0);
01338 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
01339 accelMDdVAverage += dV;
01340 }
01341
01342 } else if (simParams->accelMDdual) {
01343
01344 testV = dihedralEnergy + crosstermEnergy;
01345 if ( testV < accelMDE ) {
01346 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
01347 factor_dihe *= factor_dihe;
01348 vir = vir_dihe * (factor_dihe - 1.0) ;
01349 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
01350 }
01351
01352 testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
01353 if ( testV < accelMDTE ) {
01354 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
01355 factor_tot *= factor_tot;
01356 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
01357 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
01358 }
01359 accelMDdVAverage += dV;
01360
01361 } else {
01362
01363 testV = potentialEnergy;
01364 if ( testV < accelMDE ) {
01365 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
01366 factor_tot *= factor_tot;
01367 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
01368 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
01369 accelMDdVAverage += dV;
01370 }
01371 }
01372
01373 accelMDfactor[0]=factor_dihe;
01374 accelMDfactor[1]=factor_tot;
01375 accelMDfactor[2]=1;
01376 broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
01377 virial_amd = vir;
01378
01379 if ( factor_tot < 0.001 ) {
01380 iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
01381 << "\n" << endi;
01382 }
01383
01384 if ( ! (step % accelMDOutFreq) ) {
01385 if ( !(step == simParams->firstTimestep) ) {
01386 accelMDdVAverage = accelMDdVAverage/accelMDOutFreq;
01387 }
01388 iout << "ACCELERATED MD: STEP " << step
01389 << " dV " << dV
01390 << " dVAVG " << accelMDdVAverage
01391 << " BOND " << bondEnergy
01392 << " ANGLE " << angleEnergy
01393 << " DIHED " << dihedralEnergy+crosstermEnergy
01394 << " IMPRP " << improperEnergy
01395 << " ELECT " << amd_electEnergy+amd_electEnergySlow
01396 << " VDW " << amd_ljEnergy
01397 << " POTENTIAL " << potentialEnergy << "\n"
01398 << endi;
01399
01400 accelMDdVAverage = 0;
01401
01402 if (simParams->accelMDDebugOn) {
01403 BigReal volume;
01404 Tensor p_normal;
01405 Tensor p_nbond;
01406 Tensor p_slow;
01407 Tensor p;
01408 if ( (volume=lattice.volume()) != 0. ) {
01409 p_normal = vir_normal/volume;
01410 p_nbond = vir_nbond/volume;
01411 p_slow = vir_slow/volume;
01412 p = vir/volume;
01413 }
01414 iout << " accelMD Scaling Factor: " << accelMDfactor << "\n"
01415 << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
01416 << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
01417 << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
01418 << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n"
01419 << endi;
01420 }
01421 }
01422 }
01423
01424 void Controller::adaptTempInit(int step) {
01425 if (!simParams->adaptTempOn) return;
01426 iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
01427 adaptTempDtMin = 0;
01428 adaptTempDtMax = 0;
01429 adaptTempAutoDt = false;
01430 if (simParams->adaptTempBins == 0) {
01431 iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
01432 std::ifstream adaptTempRead(simParams->adaptTempInFile);
01433 if (adaptTempRead) {
01434 int readInt;
01435 BigReal readReal;
01436 bool readBool;
01437
01438 adaptTempRead >> readInt;
01439
01440 adaptTempRead >> adaptTempT;
01441 adaptTempRead >> adaptTempBetaMin;
01442 adaptTempRead >> adaptTempBetaMax;
01443 adaptTempBetaMin = 1./adaptTempBetaMin;
01444 adaptTempBetaMax = 1./adaptTempBetaMax;
01445
01446 if (adaptTempBetaMin > adaptTempBetaMax){
01447 readReal = adaptTempBetaMax;
01448 adaptTempBetaMax = adaptTempBetaMin;
01449 adaptTempBetaMin = adaptTempBetaMax;
01450 }
01451 adaptTempRead >> adaptTempBins;
01452 adaptTempRead >> adaptTempCg;
01453 adaptTempRead >> adaptTempDt;
01454 adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
01455 adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
01456 adaptTempPotEnergySamples = new int[adaptTempBins];
01457 adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
01458 adaptTempPotEnergyVar = new BigReal[adaptTempBins];
01459 adaptTempPotEnergyAve = new BigReal[adaptTempBins];
01460 adaptTempBetaN = new BigReal[adaptTempBins];
01461 adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
01462 for(int j = 0; j < adaptTempBins; ++j) {
01463 adaptTempRead >> adaptTempPotEnergyAve[j];
01464 adaptTempRead >> adaptTempPotEnergyVar[j];
01465 adaptTempRead >> adaptTempPotEnergySamples[j];
01466 adaptTempRead >> adaptTempPotEnergyAveNum[j];
01467 adaptTempRead >> adaptTempPotEnergyVarNum[j];
01468 adaptTempRead >> adaptTempPotEnergyAveDen[j];
01469 adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
01470 }
01471 adaptTempRead.close();
01472 }
01473 else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
01474 }
01475 else {
01476 adaptTempBins = simParams->adaptTempBins;
01477 adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
01478 adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
01479 adaptTempPotEnergySamples = new int[adaptTempBins];
01480 adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
01481 adaptTempPotEnergyVar = new BigReal[adaptTempBins];
01482 adaptTempPotEnergyAve = new BigReal[adaptTempBins];
01483 adaptTempBetaN = new BigReal[adaptTempBins];
01484 adaptTempBetaMax = 1./simParams->adaptTempTmin;
01485 adaptTempBetaMin = 1./simParams->adaptTempTmax;
01486 adaptTempCg = simParams->adaptTempCgamma;
01487 adaptTempDt = simParams->adaptTempDt;
01488 adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
01489 adaptTempT = simParams->initialTemp;
01490 if (simParams->langevinOn)
01491 adaptTempT = simParams->langevinTemp;
01492 else if (simParams->rescaleFreq > 0)
01493 adaptTempT = simParams->rescaleTemp;
01494 for(int j = 0; j < adaptTempBins; ++j){
01495 adaptTempPotEnergyAveNum[j] = 0.;
01496 adaptTempPotEnergyAveDen[j] = 0.;
01497 adaptTempPotEnergySamples[j] = 0;
01498 adaptTempPotEnergyVarNum[j] = 0.;
01499 adaptTempPotEnergyVar[j] = 0.;
01500 adaptTempPotEnergyAve[j] = 0.;
01501 adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
01502 }
01503 }
01504 if (simParams->adaptTempAutoDt > 0.0) {
01505 adaptTempAutoDt = true;
01506 adaptTempDtMin = simParams->adaptTempAutoDt - 0.01;
01507 adaptTempDtMax = simParams->adaptTempAutoDt + 0.01;
01508 }
01509 adaptTempDTave = 0;
01510 adaptTempDTavenum = 0;
01511 iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
01512 iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
01513 iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
01514 iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
01515 iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
01516 if (simParams->adaptTempRestartFreq > 0) {
01517 NAMD_backup_file(simParams->adaptTempRestartFile);
01518 adaptTempRestartFile.open(simParams->adaptTempRestartFile);
01519 if(!adaptTempRestartFile)
01520 NAMD_die("Error opening restart file");
01521 }
01522 }
01523
01524 void Controller::adaptTempWriteRestart(int step) {
01525 if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
01526 adaptTempRestartFile.seekp(std::ios::beg);
01527 iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
01528 adaptTempRestartFile << step << " ";
01529
01530 adaptTempRestartFile << adaptTempT<< " ";
01531 adaptTempRestartFile << 1./adaptTempBetaMin << " ";
01532 adaptTempRestartFile << 1./adaptTempBetaMax << " ";
01533 adaptTempRestartFile << adaptTempBins << " ";
01534 adaptTempRestartFile << adaptTempCg << " ";
01535 adaptTempRestartFile << adaptTempDt ;
01536 adaptTempRestartFile << "\n" ;
01537 for(int j = 0; j < adaptTempBins; ++j) {
01538 adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
01539 adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
01540 adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
01541 adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
01542 adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
01543 adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
01544 adaptTempRestartFile << "\n";
01545 }
01546 adaptTempRestartFile.flush();
01547 }
01548 }
01549
01550 void Controller::adaptTempUpdate(int step, int minimize)
01551 {
01552
01553 if ( !simParams->adaptTempOn ) return;
01554 int j = 0;
01555 if (step == simParams->firstTimestep) {
01556 adaptTempInit(step);
01557 return;
01558 }
01559 if ( minimize || (step < simParams->adaptTempFirstStep ) ||
01560 ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
01561 const int adaptTempOutFreq = simParams->adaptTempOutFreq;
01562 const bool adaptTempDebug = simParams->adaptTempDebug;
01563
01564 BigReal adaptTempBeta = 1./adaptTempT;
01565 adaptTempBin = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
01566
01567 if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
01568 iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin
01569 << " adaptTempBeta: " << adaptTempBeta
01570 << " adaptTempDBeta: " << adaptTempDBeta
01571 << " betaMin:" << adaptTempBetaMin
01572 << " betaMax: " << adaptTempBetaMax << "\n";
01573 adaptTempPotEnergySamples[adaptTempBin] += 1;
01574 BigReal gammaAve = 1.-adaptTempCg/adaptTempPotEnergySamples[adaptTempBin];
01575
01576 BigReal potentialEnergy;
01577 BigReal potEnergyAverage;
01578 BigReal potEnergyVariance;
01579 potentialEnergy = totalEnergy-kineticEnergy;
01580
01581
01582 adaptTempPotEnergyAveNum[adaptTempBin] = adaptTempPotEnergyAveNum[adaptTempBin]*gammaAve + potentialEnergy;
01583 adaptTempPotEnergyAveDen[adaptTempBin] = adaptTempPotEnergyAveDen[adaptTempBin]*gammaAve + 1;
01584 adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
01585
01586 potEnergyAverage = adaptTempPotEnergyAveNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
01587 potEnergyVariance = adaptTempPotEnergyVarNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
01588 potEnergyVariance -= potEnergyAverage*potEnergyAverage;
01589
01590 adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
01591 adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
01592
01593
01594
01595
01596
01597 if ( ! ( step % simParams->adaptTempFreq ) ) {
01598
01599 if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
01600
01601 BigReal deltaBeta = 0.04*adaptTempBeta;
01602 BigReal betaPlus;
01603 BigReal betaMinus;
01604 int nMinus =0;
01605 int nPlus = 0;
01606 if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
01607 if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
01608 betaMinus = adaptTempBeta - deltaBeta;
01609 betaPlus = adaptTempBeta + deltaBeta;
01610 BigReal betaMinus2 = betaMinus*betaMinus;
01611 BigReal betaPlus2 = betaPlus*betaPlus;
01612 nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
01613 nPlus = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
01614
01615 BigReal potEnergyAve0 = 0.0;
01616 BigReal potEnergyAve1 = 0.0;
01617
01618 BigReal A0 = 0;
01619 BigReal A1 = 0;
01620 BigReal A2 = 0;
01621
01622 BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1];
01623 BigReal DeltaE2Ave;
01624 BigReal DeltaE2AveSum = 0;
01625 BigReal betaj;
01626 for (j = nMinus+1; j <= adaptTempBin; ++j) {
01627 potEnergyAve0 += adaptTempPotEnergyAve[j];
01628 DeltaE2Ave = adaptTempPotEnergyVar[j];
01629 DeltaE2AveSum += DeltaE2Ave;
01630 A0 += j*DeltaE2Ave;
01631 }
01632 A0 *= adaptTempDBeta;
01633 A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
01634 A0 *= adaptTempDBeta;
01635 betaj = adaptTempBetaN[nMinus+1]-betaMinus;
01636 betaj *= betaj;
01637 A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
01638 A0 /= (betaNp1-betaMinus);
01639
01640
01641 DeltaE2AveSum = 0;
01642 for (j = adaptTempBin+1; j < nPlus; j++) {
01643 potEnergyAve1 += adaptTempPotEnergyAve[j];
01644 DeltaE2Ave = adaptTempPotEnergyVar[j];
01645 DeltaE2AveSum += DeltaE2Ave;
01646 A1 += j*DeltaE2Ave;
01647 }
01648 A1 *= adaptTempDBeta;
01649 A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
01650 A1 *= adaptTempDBeta;
01651 if ( nPlus < adaptTempBins ) {
01652 betaj = betaPlus-adaptTempBetaN[nPlus];
01653 betaj *= betaj;
01654 A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
01655 }
01656 A1 /= (betaPlus-betaNp1);
01657
01658
01659 A2 = 0.5*adaptTempDBeta*potEnergyVariance;
01660
01661
01662 DeltaE2Ave = A0-A1;
01663 BigReal aplus = 0;
01664 BigReal aminus = 0;
01665 if (DeltaE2Ave != 0) {
01666 aplus = (A0-A2)/(A0-A1);
01667 if (aplus < 0) {
01668 aplus = 0;
01669 }
01670 if (aplus > 1) {
01671 aplus = 1;
01672 }
01673 aminus = 1-aplus;
01674 potEnergyAve0 *= adaptTempDBeta;
01675 potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
01676 potEnergyAve0 /= (betaNp1-betaMinus);
01677 potEnergyAve1 *= adaptTempDBeta;
01678 if ( nPlus < adaptTempBins ) {
01679 potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
01680 }
01681 potEnergyAve1 /= (betaPlus-betaNp1);
01682 potEnergyAverage = aminus*potEnergyAve0;
01683 potEnergyAverage += aplus*potEnergyAve1;
01684 }
01685 if (simParams->adaptTempDebug) {
01686 iout << "ADAPTEMP DEBUG:" << "\n"
01687 << " adaptTempBin: " << adaptTempBin << "\n"
01688 << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
01689 << " adaptTempBeta: " << adaptTempBeta << "\n"
01690 << " adaptTemp: " << adaptTempT<< "\n"
01691 << " betaMin: " << adaptTempBetaMin << "\n"
01692 << " betaMax: " << adaptTempBetaMax << "\n"
01693 << " gammaAve: " << gammaAve << "\n"
01694 << " deltaBeta: " << deltaBeta << "\n"
01695 << " betaMinus: " << betaMinus << "\n"
01696 << " betaPlus: " << betaPlus << "\n"
01697 << " nMinus: " << nMinus << "\n"
01698 << " nPlus: " << nPlus << "\n"
01699 << " A0: " << A0 << "\n"
01700 << " A1: " << A1 << "\n"
01701 << " A2: " << A2 << "\n"
01702 << " a+: " << aplus << "\n"
01703 << " a-: " << aminus << "\n"
01704 << endi;
01705 }
01706 }
01707 else {
01708 if (simParams->adaptTempDebug) {
01709 iout << "ADAPTEMP DEBUG:" << "\n"
01710 << " adaptTempBin: " << adaptTempBin << "\n"
01711 << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
01712 << " adaptTempBeta: " << adaptTempBeta << "\n"
01713 << " adaptTemp: " << adaptTempT<< "\n"
01714 << " betaMin: " << adaptTempBetaMin << "\n"
01715 << " betaMax: " << adaptTempBetaMax << "\n"
01716 << " gammaAve: " << gammaAve << "\n"
01717 << " deltaBeta: N/A\n"
01718 << " betaMinus: N/A\n"
01719 << " betaPlus: N/A\n"
01720 << " nMinus: N/A\n"
01721 << " nPlus: N/A\n"
01722 << " A0: N/A\n"
01723 << " A1: N/A\n"
01724 << " A2: N/A\n"
01725 << " a+: N/A\n"
01726 << " a-: N/A\n"
01727 << endi;
01728 }
01729 }
01730
01731
01732 BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
01733 dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
01734 dT += adaptTempT;
01735
01736
01737
01738 if ( dT > 1./adaptTempBetaMin || dT < 1./adaptTempBetaMax ) {
01739 dT = ((potentialEnergy-adaptTempPotEnergyAve[adaptTempBin])/BOLTZMANN+adaptTempT)*adaptTempDt;
01740 dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
01741 dT += adaptTempT;
01742
01743 if ( dT > 1./adaptTempBetaMin ) {
01744 if (!simParams->adaptTempRandom) {
01745
01746
01747
01748
01749 dT = adaptTempT;
01750 }
01751 else {
01752
01753
01754
01755 dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
01756 dT = 1./dT;
01757 }
01758 }
01759 else if ( dT < 1./adaptTempBetaMax ) {
01760 if (!simParams->adaptTempRandom) {
01761
01762
01763
01764 dT = adaptTempT;
01765 }
01766 else {
01767
01768
01769
01770 dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
01771 dT = 1./dT;
01772 }
01773 }
01774 else if (adaptTempAutoDt) {
01775
01776
01777 BigReal adaptTempTdiff = fabs(dT-adaptTempT);
01778 if (adaptTempTdiff > 0) {
01779 adaptTempDTave += adaptTempTdiff;
01780 adaptTempDTavenum++;
01781
01782 }
01783 if(adaptTempDTavenum == 100){
01784 BigReal Frac;
01785 adaptTempDTave /= adaptTempDTavenum;
01786 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
01787 Frac = adaptTempDTave/Frac;
01788
01789
01790 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
01791 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
01792 Frac = adaptTempDtMax/Frac;
01793 iout << "ADAPTEMP: Updating adaptTempDt to ";
01794 adaptTempDt *= Frac;
01795 iout << adaptTempDt << "\n" << endi;
01796 }
01797 adaptTempDTave = 0;
01798 adaptTempDTavenum = 0;
01799 }
01800 }
01801 }
01802 else if (adaptTempAutoDt) {
01803
01804
01805 BigReal adaptTempTdiff = fabs(dT-adaptTempT);
01806 if (adaptTempTdiff > 0) {
01807 adaptTempDTave += adaptTempTdiff;
01808 adaptTempDTavenum++;
01809
01810 }
01811 if(adaptTempDTavenum == 100){
01812 BigReal Frac;
01813 adaptTempDTave /= adaptTempDTavenum;
01814 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
01815 Frac = adaptTempDTave/Frac;
01816
01817
01818 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
01819 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
01820 Frac = adaptTempDtMax/Frac;
01821 iout << "ADAPTEMP: Updating adaptTempDt to ";
01822 adaptTempDt *= Frac;
01823 iout << adaptTempDt << "\n" << endi;
01824 }
01825 adaptTempDTave = 0;
01826 adaptTempDTavenum = 0;
01827
01828 }
01829
01830 }
01831
01832 adaptTempT = dT;
01833 broadcast->adaptTemperature.publish(step,adaptTempT);
01834 }
01835 adaptTempWriteRestart(step);
01836 if ( ! (step % adaptTempOutFreq) ) {
01837 iout << "ADAPTEMP: STEP " << step
01838 << " TEMP " << adaptTempT
01839 << " ENERGY " << std::setprecision(10) << potentialEnergy
01840 << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
01841 << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
01842 iout << "\n" << endi;
01843 }
01844
01845 }
01846
01847
01848 void Controller::compareChecksums(int step, int forgiving) {
01849 Node *node = Node::Object();
01850 Molecule *molecule = node->molecule;
01851
01852
01853 BigReal checksum, checksum_b;
01854
01855 checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
01856 if ( ((int)checksum) != molecule->numAtoms ) {
01857 char errmsg[256];
01858 sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
01859 (int)checksum, molecule->numAtoms);
01860 NAMD_bug(errmsg);
01861 }
01862
01863 checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
01864 if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
01865 if ( ! computeChecksum ) {
01866 computesPartitioned = 0;
01867 } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
01868 computesPartitioned = 1;
01869 } else {
01870 NAMD_bug("Bad global compute count!\n");
01871 }
01872 computeChecksum = ((int)checksum);
01873 }
01874
01875 checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
01876 if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
01877 if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
01878 iout << iWARN << "Bad global bond count!\n" << endi;
01879 else NAMD_bug("Bad global bond count!\n");
01880 }
01881
01882 checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
01883 if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
01884 if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
01885 iout << iWARN << "Bad global angle count!\n" << endi;
01886 else NAMD_bug("Bad global angle count!\n");
01887 }
01888
01889 checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
01890 if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
01891 if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
01892 iout << iWARN << "Bad global dihedral count!\n" << endi;
01893 else NAMD_bug("Bad global dihedral count!\n");
01894 }
01895
01896 checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
01897 if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
01898 if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
01899 iout << iWARN << "Bad global improper count!\n" << endi;
01900 else NAMD_bug("Bad global improper count!\n");
01901 }
01902
01903 checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
01904 if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
01905 if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
01906 iout << iWARN << "Bad global Thole count!\n" << endi;
01907 else NAMD_bug("Bad global Thole count!\n");
01908 }
01909
01910 checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
01911 if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
01912 if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
01913 iout << iWARN << "Bad global Aniso count!\n" << endi;
01914 else NAMD_bug("Bad global Aniso count!\n");
01915 }
01916
01917 checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
01918 if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
01919 if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
01920 iout << iWARN << "Bad global crossterm count!\n" << endi;
01921 else NAMD_bug("Bad global crossterm count!\n");
01922 }
01923
01924 checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01925 if ( ((int)checksum) > molecule->numCalcExclusions &&
01926 ( ! simParams->mollyOn || step % slowFreq ) ) {
01927 if ( forgiving )
01928 iout << iWARN << "High global exclusion count ("
01929 << ((int)checksum) << " vs "
01930 << molecule->numCalcExclusions << "), possible error!\n"
01931 << iWARN << "This warning is not unusual during minimization.\n"
01932 << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
01933 else {
01934 char errmsg[256];
01935 sprintf(errmsg, "High global exclusion count! (%d vs %d) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
01936 (int)checksum, molecule->numCalcExclusions);
01937 NAMD_bug(errmsg);
01938 }
01939 }
01940 if ( ((int)checksum) && ((int)checksum) < molecule->numCalcExclusions ) {
01941 if ( forgiving || ! simParams->fullElectFrequency ) {
01942 iout << iWARN << "Low global exclusion count! ("
01943 << ((int)checksum) << " vs " << molecule->numCalcExclusions << ")";
01944 if ( forgiving ) {
01945 iout << "\n"
01946 << iWARN << "This warning is not unusual during minimization.\n"
01947 << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
01948 } else {
01949 iout << " System unstable or pairlistdist or cutoff too small.\n";
01950 }
01951 iout << endi;
01952 } else {
01953 char errmsg[256];
01954 sprintf(errmsg, "Low global exclusion count! (%d vs %d) System unstable or pairlistdist or cutoff too small.\n",
01955 (int)checksum, molecule->numCalcExclusions);
01956 NAMD_bug(errmsg);
01957 }
01958 }
01959
01960 checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
01961 if ( ((int)checksum) && ! marginViolations ) {
01962 iout << iERROR << "Margin is too small for " << ((int)checksum) <<
01963 " atoms during timestep " << step << ".\n" << iERROR <<
01964 "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
01965 }
01966 marginViolations += (int)checksum;
01967
01968 checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
01969 if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
01970 iout << iINFO <<
01971 "Pairlistdist is too small for " << ((int)checksum) <<
01972 " computes during timestep " << step << ".\n" << endi;
01973 }
01974 if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
01975
01976 checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
01977 if ( checksum ) {
01978 if ( forgiving )
01979 iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
01980 else NAMD_bug("Stray PME grid charges detected!\n");
01981 }
01982 }
01983
01984 void Controller::printTiming(int step) {
01985
01986 if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
01987 {
01988 const double endWTime = CmiWallTimer() - firstWTime;
01989 const double endCTime = CmiTimer() - firstCTime;
01990
01991
01992 if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
01993
01994 const double elapsedW =
01995 (endWTime - startWTime) / simParams->outputTiming;
01996 const double elapsedC =
01997 (endCTime - startCTime) / simParams->outputTiming;
01998
01999 const double remainingW = elapsedW * (simParams->N - step);
02000 const double remainingW_hours = remainingW / 3600;
02001
02002 startWTime = endWTime;
02003 startCTime = endCTime;
02004
02005 #ifdef NAMD_CUDA
02006 if ( simParams->outputEnergies < 60 &&
02007 step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
02008 iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
02009 }
02010 #endif
02011 if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
02012 CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
02013 ", %g hours remaining, %f MB of memory in use.\n",
02014 step, endCTime, elapsedC, endWTime, elapsedW,
02015 remainingW_hours, memusage_MB());
02016 if ( fflush_count ) { --fflush_count; fflush(stdout); }
02017 }
02018 }
02019 }
02020
02021 void Controller::printMinimizeEnergies(int step) {
02022
02023 rescaleaccelMD(step,1);
02024 receivePressure(step,1);
02025
02026 Node *node = Node::Object();
02027 Molecule *molecule = node->molecule;
02028 compareChecksums(step,1);
02029
02030 printEnergies(step,1);
02031
02032 min_energy = totalEnergy;
02033 min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
02034 min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
02035 min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
02036 min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
02037
02038 }
02039
02040 void Controller::printDynamicsEnergies(int step) {
02041
02042 Node *node = Node::Object();
02043 Molecule *molecule = node->molecule;
02044 SimParameters *simParameters = node->simParameters;
02045 Lattice &lattice = state->lattice;
02046
02047 compareChecksums(step);
02048
02049 printEnergies(step,0);
02050 }
02051
02052 void Controller::printEnergies(int step, int minimize)
02053 {
02054 Node *node = Node::Object();
02055 Molecule *molecule = node->molecule;
02056 SimParameters *simParameters = node->simParameters;
02057 Lattice &lattice = state->lattice;
02058
02059
02060
02061
02062 BigReal bondEnergy;
02063 BigReal angleEnergy;
02064 BigReal dihedralEnergy;
02065 BigReal improperEnergy;
02066 BigReal crosstermEnergy;
02067 BigReal boundaryEnergy;
02068 BigReal miscEnergy;
02069 BigReal potentialEnergy;
02070 BigReal flatEnergy;
02071 BigReal smoothEnergy;
02072
02073 Vector momentum;
02074 Vector angularMomentum;
02075 BigReal volume = lattice.volume();
02076
02077 bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
02078 angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
02079 dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
02080 improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
02081 crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
02082 boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
02083 miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
02084
02085 if ( minimize || ! ( step % nbondFreq ) )
02086 {
02087 electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
02088 ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
02089
02090
02091 goNativeEnergy = reduction->item(REDUCTION_GO_NATIVE_ENERGY);
02092 goNonnativeEnergy = reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
02093 goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
02094
02095
02096 electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
02097 ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
02098
02099 electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
02100 ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
02101 electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
02102 ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
02103
02104 }
02105
02106 if ( minimize || ! ( step % slowFreq ) )
02107 {
02108 electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
02109
02110 electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
02111
02112 electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
02113 electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
02114 electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
02115 electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
02116
02117 }
02118
02119 if (simParameters->LJcorrection && volume) {
02120
02121
02122
02123 ljEnergy += molecule->tail_corr_ener / volume;
02124 }
02125
02126
02127 momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
02128 momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
02129 momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
02130 angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
02131 angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
02132 angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
02133
02134
02135 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
02136 improperEnergy + electEnergy + electEnergySlow + ljEnergy +
02137 crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy;
02138
02139 totalEnergy = potentialEnergy + kineticEnergy;
02140 flatEnergy = totalEnergy +
02141 (1.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
02142 if ( !(step%slowFreq) ) {
02143
02144 BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
02145 if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
02146 if ( step != simParams->firstTimestep ) {
02147 smooth2_avg *= 0.9375;
02148 smooth2_avg += 0.0625 * s;
02149 }
02150 }
02151 smoothEnergy = flatEnergy + smooth2_avg -
02152 (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
02153
02154 if ( simParameters->outputMomenta && ! minimize &&
02155 ! ( step % simParameters->outputMomenta ) )
02156 {
02157 iout << "MOMENTUM: " << step
02158 << " P: " << momentum
02159 << " L: " << angularMomentum
02160 << "\n" << endi;
02161 }
02162
02163 if ( simParameters->outputPressure ) {
02164 pressure_tavg += pressure;
02165 groupPressure_tavg += groupPressure;
02166 tavg_count += 1;
02167 if ( minimize || ! ( step % simParameters->outputPressure ) ) {
02168 iout << "PRESSURE: " << step << " "
02169 << PRESSUREFACTOR * pressure << "\n"
02170 << "GPRESSURE: " << step << " "
02171 << PRESSUREFACTOR * groupPressure << "\n";
02172 if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
02173 << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
02174 << "GPRESSAVG: " << step << " "
02175 << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
02176 iout << endi;
02177 pressure_tavg = 0;
02178 groupPressure_tavg = 0;
02179 tavg_count = 0;
02180 }
02181 }
02182
02183
02184 if (pressureProfileSlabs) {
02185 const int freq = simParams->pressureProfileFreq;
02186 const int arraysize = 3*pressureProfileSlabs;
02187
02188 BigReal *total = new BigReal[arraysize];
02189 memset(total, 0, arraysize*sizeof(BigReal));
02190 const int first = simParams->firstTimestep;
02191
02192 if (ppbonded) ppbonded->getData(first, step, lattice, total);
02193 if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
02194 if (ppint) ppint->getData(first, step, lattice, total);
02195 for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
02196 pressureProfileCount++;
02197
02198 if (!(step % freq)) {
02199
02200 BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
02201
02202 iout << "PRESSUREPROFILE: " << step << " ";
02203 if (step == first) {
02204
02205 for (int i=0; i<arraysize; i++) {
02206 iout << total[i] * scalefac << " ";
02207 }
02208 } else {
02209
02210 scalefac /= pressureProfileCount;
02211 for (int i=0; i<arraysize; i++)
02212 iout << pressureProfileAverage[i]*scalefac << " ";
02213 }
02214 iout << "\n" << endi;
02215
02216
02217 memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
02218 pressureProfileCount = 0;
02219 }
02220 delete [] total;
02221 }
02222
02223 int stepInRun = step - simParams->firstTimestep;
02224 if ( stepInRun % simParams->firstLdbStep == 0 ) {
02225 int benchPhase = stepInRun / simParams->firstLdbStep;
02226 if ( benchPhase > 0 && benchPhase < 10 ) {
02227 #ifdef NAMD_CUDA
02228 if ( simParams->outputEnergies < 60 ) {
02229 iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
02230 }
02231 #endif
02232 iout << iINFO;
02233 if ( benchPhase < 4 ) iout << "Initial time: ";
02234 else iout << "Benchmark time: ";
02235 iout << CkNumPes() << " CPUs ";
02236 {
02237 BigReal wallPerStep =
02238 (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
02239 BigReal ns = simParams->dt / 1000000.0;
02240 BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
02241 BigReal daysPerNano = wallPerStep * days / ns;
02242 iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
02243 iout << memusage_MB() << " MB memory\n" << endi;
02244 }
02245 }
02246 startBenchTime = CmiWallTimer();
02247 }
02248
02249 printTiming(step);
02250
02251 Vector pairVDWForce, pairElectForce;
02252 if ( simParameters->pairInteractionOn ) {
02253 GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
02254 GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
02255 }
02256
02257
02258 #ifdef NAMD_TCL
02259 #define CALLBACKDATA(LABEL,VALUE) \
02260 labels << (LABEL) << " "; values << (VALUE) << " ";
02261 #define CALLBACKLIST(LABEL,VALUE) \
02262 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
02263 if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
02264 std::ostringstream labels, values;
02265 #if CMK_BLUEGENEL
02266
02267 values.precision(16);
02268 #else
02269 values << std::setprecision(16);
02270 #endif
02271 CALLBACKDATA("TS",step);
02272 CALLBACKDATA("BOND",bondEnergy);
02273 CALLBACKDATA("ANGLE",angleEnergy);
02274 CALLBACKDATA("DIHED",dihedralEnergy);
02275 CALLBACKDATA("CROSS",crosstermEnergy);
02276 CALLBACKDATA("IMPRP",improperEnergy);
02277 CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
02278 CALLBACKDATA("VDW",ljEnergy);
02279 CALLBACKDATA("BOUNDARY",boundaryEnergy);
02280 CALLBACKDATA("MISC",miscEnergy);
02281 CALLBACKDATA("POTENTIAL",potentialEnergy);
02282 CALLBACKDATA("KINETIC",kineticEnergy);
02283 CALLBACKDATA("TOTAL",totalEnergy);
02284 CALLBACKDATA("TEMP",temperature);
02285 CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
02286 CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
02287 CALLBACKDATA("VOLUME",lattice.volume());
02288 CALLBACKLIST("CELL_A",lattice.a());
02289 CALLBACKLIST("CELL_B",lattice.b());
02290 CALLBACKLIST("CELL_C",lattice.c());
02291 CALLBACKLIST("CELL_O",lattice.origin());
02292 labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
02293 << lattice.b_p() << " " << lattice.c_p() << "} ";
02294 if ( simParameters->pairInteractionOn ) {
02295 CALLBACKLIST("VDW_FORCE",pairVDWForce);
02296 CALLBACKLIST("ELECT_FORCE",pairElectForce);
02297 }
02298
02299 labels << '\0'; values << '\0';
02300 state->callback_labelstring = labels.str();
02301 state->callback_valuestring = values.str();
02302
02303 }
02304 #undef CALLBACKDATA
02305 #endif
02306
02307 drudeComTempAvg += drudeComTemp;
02308 drudeBondTempAvg += drudeBondTemp;
02309
02310 temp_avg += temperature;
02311 pressure_avg += trace(pressure)/3.;
02312 groupPressure_avg += trace(groupPressure)/3.;
02313 avg_count += 1;
02314
02315 if ( simParams->outputPairlists && pairlistWarnings &&
02316 ! (step % simParams->outputPairlists) ) {
02317 iout << iINFO << pairlistWarnings <<
02318 " pairlist warnings in past " << simParams->outputPairlists <<
02319 " steps.\n" << endi;
02320 pairlistWarnings = 0;
02321 }
02322
02323
02324 if ( ! minimize && step % simParameters->outputEnergies ) return;
02325
02326
02327 if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
02328 IMDEnergies energies;
02329 energies.tstep = step;
02330 energies.T = temp_avg/avg_count;
02331 energies.Etot = totalEnergy;
02332 energies.Epot = potentialEnergy;
02333 energies.Evdw = ljEnergy;
02334 energies.Eelec = electEnergy + electEnergySlow;
02335 energies.Ebond = bondEnergy;
02336 energies.Eangle = angleEnergy;
02337 energies.Edihe = dihedralEnergy + crosstermEnergy;
02338 energies.Eimpr = improperEnergy;
02339 Node::Object()->imd->gather_energies(&energies);
02340 }
02341
02342 if ( marginViolations ) {
02343 iout << iERROR << marginViolations <<
02344 " margin violations detected since previous energy output.\n" << endi;
02345 }
02346 marginViolations = 0;
02347
02348 if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
02349 {
02350 iout << "ETITLE: TS";
02351 iout << FORMAT("BOND");
02352 iout << FORMAT("ANGLE");
02353 iout << FORMAT("DIHED");
02354 if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
02355 iout << FORMAT("IMPRP");
02356 iout << " ";
02357 iout << FORMAT("ELECT");
02358 iout << FORMAT("VDW");
02359 iout << FORMAT("BOUNDARY");
02360 iout << FORMAT("MISC");
02361 iout << FORMAT("KINETIC");
02362 iout << " ";
02363 iout << FORMAT("TOTAL");
02364 iout << FORMAT("TEMP");
02365 iout << FORMAT("POTENTIAL");
02366
02367 iout << FORMAT("TOTAL3");
02368 iout << FORMAT("TEMPAVG");
02369 if ( volume != 0. ) {
02370 iout << " ";
02371 iout << FORMAT("PRESSURE");
02372 iout << FORMAT("GPRESSURE");
02373 iout << FORMAT("VOLUME");
02374 iout << FORMAT("PRESSAVG");
02375 iout << FORMAT("GPRESSAVG");
02376 }
02377 if (simParameters->drudeOn) {
02378 iout << " ";
02379 iout << FORMAT("DRUDECOM");
02380 iout << FORMAT("DRUDEBOND");
02381 iout << FORMAT("DRCOMAVG");
02382 iout << FORMAT("DRBONDAVG");
02383 }
02384
02385 if (simParameters->goForcesOn) {
02386 iout << " ";
02387 iout << FORMAT("NATIVE");
02388 iout << FORMAT("NONNATIVE");
02389
02390
02391 iout << FORMAT("GOTOTAL");
02392
02393 }
02394
02395 iout << "\n\n" << endi;
02396 }
02397
02398
02399
02400
02401 iout << ETITLE(step);
02402 iout << FORMAT(bondEnergy);
02403 iout << FORMAT(angleEnergy);
02404 if ( simParameters->mergeCrossterms ) {
02405 iout << FORMAT(dihedralEnergy+crosstermEnergy);
02406 } else {
02407 iout << FORMAT(dihedralEnergy);
02408 iout << FORMAT(crosstermEnergy);
02409 }
02410 iout << FORMAT(improperEnergy);
02411 iout << " ";
02412 iout << FORMAT(electEnergy+electEnergySlow);
02413 iout << FORMAT(ljEnergy);
02414 iout << FORMAT(boundaryEnergy);
02415 iout << FORMAT(miscEnergy);
02416 iout << FORMAT(kineticEnergy);
02417 iout << " ";
02418 iout << FORMAT(totalEnergy);
02419 iout << FORMAT(temperature);
02420 iout << FORMAT(potentialEnergy);
02421
02422 iout << FORMAT(smoothEnergy);
02423 iout << FORMAT(temp_avg/avg_count);
02424 if ( volume != 0. )
02425 {
02426 iout << " ";
02427 iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
02428 iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
02429 iout << FORMAT(volume);
02430 iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
02431 iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
02432 }
02433 if (simParameters->drudeOn) {
02434 iout << " ";
02435 iout << FORMAT(drudeComTemp);
02436 iout << FORMAT(drudeBondTemp);
02437 iout << FORMAT(drudeComTempAvg/avg_count);
02438 iout << FORMAT(drudeBondTempAvg/avg_count);
02439 }
02440
02441 if (simParameters->goForcesOn) {
02442 iout << " ";
02443 iout << FORMAT(goNativeEnergy);
02444 iout << FORMAT(goNonnativeEnergy);
02445
02446
02447 iout << FORMAT(goTotalEnergy);
02448
02449 }
02450 iout << "\n\n" << endi;
02451
02452 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
02453 char webout[80];
02454 sprintf(webout,"%d %d %d %d",(int)totalEnergy,
02455 (int)(potentialEnergy),
02456 (int)kineticEnergy,(int)temperature);
02457 CApplicationDepositNode0Data(webout);
02458 #endif
02459
02460 if (simParameters->pairInteractionOn) {
02461 iout << "PAIR INTERACTION:";
02462 iout << " STEP: " << step;
02463 iout << " VDW_FORCE: ";
02464 iout << FORMAT(pairVDWForce.x);
02465 iout << FORMAT(pairVDWForce.y);
02466 iout << FORMAT(pairVDWForce.z);
02467 iout << " ELECT_FORCE: ";
02468 iout << FORMAT(pairElectForce.x);
02469 iout << FORMAT(pairElectForce.y);
02470 iout << FORMAT(pairElectForce.z);
02471 iout << "\n" << endi;
02472 }
02473 drudeComTempAvg = 0;
02474 drudeBondTempAvg = 0;
02475 temp_avg = 0;
02476 pressure_avg = 0;
02477 groupPressure_avg = 0;
02478 avg_count = 0;
02479
02480 if ( fflush_count ) {
02481 --fflush_count;
02482 fflush(stdout);
02483 }
02484 }
02485
02486 void Controller::writeExtendedSystemLabels(std::ofstream &file) {
02487 Lattice &lattice = state->lattice;
02488 file << "#$LABELS step";
02489 if ( lattice.a_p() ) file << " a_x a_y a_z";
02490 if ( lattice.b_p() ) file << " b_x b_y b_z";
02491 if ( lattice.c_p() ) file << " c_x c_y c_z";
02492 file << " o_x o_y o_z";
02493 if ( simParams->langevinPistonOn ) {
02494 file << " s_x s_y s_z s_u s_v s_w";
02495 }
02496 file << std::endl;
02497 }
02498
02499 void Controller::writeExtendedSystemData(int step, std::ofstream &file) {
02500 Lattice &lattice = state->lattice;
02501 file << step;
02502 if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
02503 if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
02504 if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
02505 file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
02506 if ( simParams->langevinPistonOn ) {
02507 Vector strainRate = diagonal(langevinPiston_strainRate);
02508 Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
02509 file << " " << strainRate.x;
02510 file << " " << strainRate.y;
02511 file << " " << strainRate.z;
02512 file << " " << strainRate2.x;
02513 file << " " << strainRate2.y;
02514 file << " " << strainRate2.z;
02515 }
02516 file << std::endl;
02517 }
02518
02519 void Controller::enqueueCollections(int timestep)
02520 {
02521 if ( Output::coordinateNeeded(timestep) ){
02522 collection->enqueuePositions(timestep,state->lattice);
02523 }
02524 if ( Output::velocityNeeded(timestep) )
02525 collection->enqueueVelocities(timestep);
02526 if ( Output::forceNeeded(timestep) )
02527 collection->enqueueForces(timestep);
02528 }
02529
02530
02531 static char *FEPTITLE(int X)
02532 {
02533 static char tmp_string[21];
02534 sprintf(tmp_string, "FepEnergy: %6d ",X);
02535 return tmp_string;
02536 }
02537
02538 static char *TITITLE(int X)
02539 {
02540 static char tmp_string[21];
02541 sprintf(tmp_string, "TI %6d ",X);
02542 return tmp_string;
02543 }
02544
02545 void Controller::outputFepEnergy(int step) {
02546 if (simParams->alchFepOn) {
02547 const int stepInRun = step - simParams->firstTimestep;
02548 const int alchEquilSteps = simParams->alchEquilSteps;
02549 const BigReal alchLambda = simParams->alchLambda;
02550 const BigReal alchLambda2 = simParams->alchLambda2;
02551 const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
02552 if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
02553 FepNo = 0;
02554 exp_dE_ByRT = 0.0;
02555 net_dE = 0.0;
02556 }
02557 BigReal dE = electEnergy_f + electEnergySlow_f + ljEnergy_f
02558 - (electEnergy + electEnergySlow + ljEnergy);
02559 BigReal RT = BOLTZMANN * simParams->alchTemp;
02560
02561 if (alchEnsembleAvg){
02562 FepNo++;
02563 exp_dE_ByRT += exp(-dE/RT);
02564 net_dE += dE;
02565 }
02566
02567 if (stepInRun == 0) {
02568 if (!fepFile.rdbuf()->is_open()) {
02569 NAMD_backup_file(simParams->alchOutFile);
02570 fepFile.open(simParams->alchOutFile);
02571 iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
02572 if(alchEnsembleAvg){
02573 fepSum = 0.0;
02574 fepFile << "# STEP Elec "
02575 << "vdW dE dE_avg Temp dG\n"
02576 << "# l l+dl "
02577 << " l l+dl E(l+dl)-E(l)" << std::endl;
02578 }
02579 else{
02580 fepFile << "# STEP Elec "
02581 << "vdW dE Temp\n"
02582 << "# l l+dl "
02583 << " l l+dl E(l+dl)-E(l)" << std::endl;
02584 }
02585 }
02586 if(!step){
02587 fepFile << "#NEW FEP WINDOW: "
02588 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
02589 << alchLambda2 << std::endl;
02590 }
02591 }
02592 if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
02593 fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
02594 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
02595 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
02596 }
02597 if ((simParams->N) && (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0))) {
02598 writeFepEnergyData(step, fepFile);
02599 fepFile.flush();
02600 }
02601 if (alchEnsembleAvg && (step == simParams->N)) {
02602 fepSum = fepSum + dG;
02603 fepFile << "#Free energy change for lambda window [ " << alchLambda
02604 << " " << alchLambda2 << " ] is " << dG << " ; net change until now is " << fepSum << std::endl;
02605 }
02606 }
02607 }
02608
02609 void Controller::outputTiEnergy(int step) {
02610 if (simParams->alchThermIntOn) {
02611 const int stepInRun = step - simParams->firstTimestep;
02612 const int alchEquilSteps = simParams->alchEquilSteps;
02613 const BigReal alchLambda = simParams->alchLambda;
02614
02615 if (stepInRun == 0 || stepInRun == alchEquilSteps) {
02616 TiNo = 0;
02617 net_dEdl_elec_1 = 0;
02618 net_dEdl_elec_2 = 0;
02619 net_dEdl_lj_1 = 0;
02620 net_dEdl_lj_2 = 0;
02621 }
02622 if (stepInRun == 0 || (! ((step - 1) % simParams->alchOutFreq))) {
02623
02624
02625 recent_TiNo = 0;
02626 recent_dEdl_elec_1 = 0;
02627 recent_dEdl_elec_2 = 0;
02628 recent_dEdl_lj_1 = 0;
02629 recent_dEdl_lj_2 = 0;
02630 }
02631 TiNo++;
02632 recent_TiNo++;
02633
02634
02635 net_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1;
02636 net_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2;
02637 net_dEdl_lj_1 += ljEnergy_ti_1;
02638 net_dEdl_lj_2 += ljEnergy_ti_2;
02639 recent_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1;
02640 recent_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2;
02641 recent_dEdl_lj_1 += ljEnergy_ti_1;
02642 recent_dEdl_lj_2 += ljEnergy_ti_2;
02643
02644 if (stepInRun == 0) {
02645 BigReal alchElecLambdaStart = simParams->alchElecLambdaStart;
02646 BigReal alchVdwLambdaEnd = simParams->alchVdwLambdaEnd;
02647 BigReal elec_lambda_1 = (alchLambda <= alchElecLambdaStart)? 0. : \
02648 (alchLambda - alchElecLambdaStart) / (1. - alchElecLambdaStart);
02649 BigReal elec_lambda_2 = ((1-alchLambda) <= alchElecLambdaStart)? 0. : \
02650 ((1-alchLambda) - alchElecLambdaStart) / (1. - alchElecLambdaStart);
02651 BigReal vdw_lambda_1 = (alchLambda >= alchVdwLambdaEnd)? 1. : \
02652 alchLambda / alchVdwLambdaEnd;
02653 BigReal vdw_lambda_2 = ((1-alchLambda) >= alchVdwLambdaEnd)? 1. : \
02654 (1-alchLambda) / alchVdwLambdaEnd;
02655 if (!tiFile.rdbuf()->is_open()) {
02656
02657 NAMD_backup_file(simParams->alchOutFile);
02658 tiFile.open(simParams->alchOutFile);
02659 iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
02660 tiFile << "# STEP Elec_dU/dl Elec_avg vdW_dU/dl vdw_avg Elec_dU/dl Elec_avg vdW_dU/dl vdw_avg PME_dU/dl PME_avg\n"
02661 << "# <---------------------PARTITION 1------------------------> <---------------------PARTITION 2--------------------->"
02662 << std::endl;
02663 }
02664 tiFile << "#NEW TI WINDOW: "
02665 << "LAMBDA " << alchLambda
02666 << "\n#PARTITION 1 VDW LAMBDA " << vdw_lambda_1
02667 << "\n#PARTITION 1 ELEC LAMBDA " << elec_lambda_1
02668 << "\n#PARTITION 2 VDW LAMBDA " << vdw_lambda_2
02669 << "\n#PARTITION 2 ELEC LAMBDA " << elec_lambda_2
02670 << "\n" << std::endl;
02671 }
02672 if (stepInRun == alchEquilSteps) {
02673 tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
02674 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
02675 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
02676 }
02677 if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) {
02678 writeTiEnergyData(step, tiFile);
02679 tiFile.flush();
02680 }
02681 }
02682 }
02683
02684 void Controller::writeFepEnergyData(int step, std::ofstream &file) {
02685 BigReal eeng = electEnergy+electEnergySlow;
02686 BigReal eeng_f = electEnergy_f + electEnergySlow_f;
02687 BigReal dE = eeng_f + ljEnergy_f - eeng - ljEnergy;
02688 BigReal RT = BOLTZMANN * simParams->alchTemp;
02689 const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
02690 const int stepInRun = step - simParams->firstTimestep;
02691 if(stepInRun){
02692 fepFile << FEPTITLE(step);
02693 fepFile << FORMAT(eeng);
02694 fepFile << FORMAT(eeng_f);
02695 fepFile << FORMAT(ljEnergy);
02696 fepFile << FORMAT(ljEnergy_f);
02697 fepFile << FORMAT(dE);
02698 if(alchEnsembleAvg){
02699 BigReal dE_avg = net_dE/FepNo;
02700 fepFile << FORMAT(dE_avg);
02701 }
02702 fepFile << FORMAT(temperature);
02703 if(alchEnsembleAvg){
02704 dG = -(RT * log(exp_dE_ByRT/FepNo));
02705 fepFile << FORMAT(dG);
02706 }
02707 fepFile << std::endl;
02708 }
02709 }
02710
02711 void Controller::writeTiEnergyData(int step, std::ofstream &file) {
02712 tiFile << TITITLE(step);
02713 tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
02714 tiFile << FORMAT(net_dEdl_elec_1/TiNo);
02715 tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
02716 tiFile << FORMAT(net_dEdl_lj_1/TiNo);
02717 tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
02718 tiFile << FORMAT(net_dEdl_elec_2/TiNo);
02719 tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
02720 tiFile << FORMAT(net_dEdl_lj_2/TiNo);
02721 tiFile << std::endl;
02722 }
02723
02724 void Controller::outputExtendedSystem(int step)
02725 {
02726
02727 if ( step >= 0 ) {
02728
02729
02730 if ( simParams->xstFrequency &&
02731 ((step % simParams->xstFrequency) == 0) )
02732 {
02733 if ( ! xstFile.rdbuf()->is_open() )
02734 {
02735 iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
02736 NAMD_backup_file(simParams->xstFilename);
02737 xstFile.open(simParams->xstFilename);
02738 while (!xstFile) {
02739 if ( errno == EINTR ) {
02740 CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
02741 xstFile.clear();
02742 xstFile.open(simParams->xstFilename);
02743 continue;
02744 }
02745 char err_msg[257];
02746 sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
02747 NAMD_err(err_msg);
02748 }
02749 xstFile << "# NAMD extended system trajectory file" << std::endl;
02750 writeExtendedSystemLabels(xstFile);
02751 }
02752 writeExtendedSystemData(step,xstFile);
02753 xstFile.flush();
02754 }
02755
02756
02757
02758 if ( simParams->restartFrequency &&
02759 ((step % simParams->restartFrequency) == 0) &&
02760 (step != simParams->firstTimestep) )
02761 {
02762 iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
02763 << step << "\n" << endi;
02764 char fname[140];
02765 strcpy(fname, simParams->restartFilename);
02766 if ( simParams->restartSave ) {
02767 char timestepstr[20];
02768 sprintf(timestepstr,".%d",step);
02769 strcat(fname, timestepstr);
02770 }
02771 strcat(fname, ".xsc");
02772 NAMD_backup_file(fname,".old");
02773 std::ofstream xscFile(fname);
02774 while (!xscFile) {
02775 if ( errno == EINTR ) {
02776 CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
02777 xscFile.clear();
02778 xscFile.open(fname);
02779 continue;
02780 }
02781 char err_msg[257];
02782 sprintf(err_msg, "Error opening XSC restart file %s",fname);
02783 NAMD_err(err_msg);
02784 }
02785 xscFile << "# NAMD extended system configuration restart file" << std::endl;
02786 writeExtendedSystemLabels(xscFile);
02787 writeExtendedSystemData(step,xscFile);
02788 if (!xscFile) {
02789 char err_msg[257];
02790 sprintf(err_msg, "Error writing XSC restart file %s",fname);
02791 NAMD_err(err_msg);
02792 }
02793 }
02794
02795 }
02796
02797
02798 if (step == FILE_OUTPUT || step == END_OF_RUN)
02799 {
02800 int realstep = ( step == FILE_OUTPUT ?
02801 simParams->firstTimestep : simParams->N );
02802 iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
02803 << realstep << "\n" << endi;
02804 static char fname[140];
02805 strcpy(fname, simParams->outputFilename);
02806 strcat(fname, ".xsc");
02807 NAMD_backup_file(fname);
02808 std::ofstream xscFile(fname);
02809 while (!xscFile) {
02810 if ( errno == EINTR ) {
02811 CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
02812 xscFile.clear();
02813 xscFile.open(fname);
02814 continue;
02815 }
02816 char err_msg[257];
02817 sprintf(err_msg, "Error opening XSC output file %s",fname);
02818 NAMD_err(err_msg);
02819 }
02820 xscFile << "# NAMD extended system configuration output file" << std::endl;
02821 writeExtendedSystemLabels(xscFile);
02822 writeExtendedSystemData(realstep,xscFile);
02823 if (!xscFile) {
02824 char err_msg[257];
02825 sprintf(err_msg, "Error writing XSC output file %s",fname);
02826 NAMD_err(err_msg);
02827 }
02828 }
02829
02830
02831 if (step == END_OF_RUN) {
02832 if ( xstFile.rdbuf()->is_open() ) {
02833 xstFile.close();
02834 iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
02835 }
02836 }
02837
02838 }
02839
02840 void Controller::rebalanceLoad(int step)
02841 {
02842 if ( ! ldbSteps ) {
02843 ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
02844 }
02845 if ( ! --ldbSteps ) {
02846 startBenchTime -= CmiWallTimer();
02847 Node::Object()->outputPatchComputeMaps("before_ldb", step);
02848 LdbCoordinator::Object()->rebalance(this);
02849 startBenchTime += CmiWallTimer();
02850 fflush_count = 3;
02851 }
02852 }
02853
02854 void Controller::cycleBarrier(int doBarrier, int step) {
02855 #if USE_BARRIER
02856 if (doBarrier) {
02857 broadcast->cycleBarrier.publish(step,1);
02858 CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
02859 CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02860 }
02861 #endif
02862 }
02863
02864 void Controller::traceBarrier(int turnOnTrace, int step) {
02865 CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02866 CProxy_Node nd(CkpvAccess(BOCclass_group).node);
02867 nd.traceBarrier(turnOnTrace, step);
02868 CthSuspend();
02869 }
02870
02871 void Controller::resumeAfterTraceBarrier(int step){
02872 broadcast->traceBarrier.publish(step,1);
02873 CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02874 awaken();
02875 }
02876
02877 #ifdef MEASURE_NAMD_WITH_PAPI
02878 void Controller::papiMeasureBarrier(int turnOnMeasure, int step){
02879 CkPrintf("Cycle time at PAPI measurement sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02880 CProxy_Node nd(CkpvAccess(BOCclass_group).node);
02881 nd.papiMeasureBarrier(turnOnMeasure, step);
02882 CthSuspend();
02883 }
02884
02885 void Controller::resumeAfterPapiMeasureBarrier(int step){
02886 broadcast->papiMeasureBarrier.publish(step,1);
02887 CkPrintf("Cycle time at PAPI measurement sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02888 awaken();
02889 }
02890 #endif
02891
02892 void Controller::terminate(void) {
02893 BackEnd::awaken();
02894 CthFree(thread);
02895 CthSuspend();
02896 }
02897
