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