NAMD: Controller.C Source File

00001 
00007 /*****************************************************************************
00008  * $Source: /home/cvs/namd/cvsroot/namd2/src/Controller.C,v $
00009  * $Author: jim $
00010  * $Date: 2017/03/29 21:29:03 $
00011  * $Revision: 1.1326 $
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   // cbrt() not in math.h on goneril
00052   #define cbrt(x)  pow(x,(double)(1.0/3.0))
00053 #endif
00054 
00055 //#define DEBUG_PRESSURE
00056 #define MIN_DEBUG_LEVEL 3
00057 //#define DEBUGM
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     // accumulate the pressure profile computed for this step into the average.
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       // convert NAMD internal virial to pressure in units of bar 
00106       BigReal scalefac = PRESSUREFACTOR * nslabs;
00107  
00108       iout << "PPROFILE" << name << ": " << step << " ";
00109       if (step == firsttimestep) {
00110         // output pressure profile for this step
00111         for (i=0; i<nelements; i++) 
00112           iout << reduction->item(i)*scalefac*inv_volume << " ";
00113       } else {
00114         // output pressure profile averaged over the last Freq steps.
00115         scalefac /= count;
00116         for (i=0; i<nelements; i++) 
00117           iout << average[i]*scalefac << " ";
00118       }
00119       iout << "\n" << endi; 
00120       // Clear the average for the next block
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     // for accelMD
00145     if (simParams->accelMDOn) {
00146        // REDUCTIONS_BASIC wil contain all potential energies and arrive first
00147        amd_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
00148        // contents of amd_reduction be submitted to REDUCTIONS_AMD
00149        submit_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_AMD);
00150        // REDUCTIONS_AMD will contain Sequencer contributions and arrive second
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     // pressure profile reductions
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       // number of partitions for pairwise interactions
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         // internal partitions by atom type, but not in a pairwise fashion
00175         ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL, 
00176             nslabs, ntypes, "INTERNAL", freq);
00177       } else {
00178         // just doing Ewald, so only need nonbonded
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     // strainRate tensor is symmetric to avoid rotation
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     // create a Thread and invoke it
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         // Parameter setting already reported in ScriptTcl
00308         // Nothing to do!
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 // XXX static and global variables are unsafe for shared memory builds.
00416 // The use of global and static vars should be eliminated.
00417 //
00418 extern int eventEndOfTimeStep;
00419 
00420 // Handle SIGINT so that restart files get written completely.
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);  // only for full-step velecities
00450     rescaleaccelMD(step);
00451     adaptTempUpdate(step); // Init adaptive tempering;
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     // Handling SIGINT doesn't seem to be working on Lemieux, and it
00472     // sometimes causes the net-xxx versions of NAMD to segfault on exit, 
00473     // so disable it for now.
00474     // namd_sighandler_t oldhandler = signal(SIGINT, 
00475     //  (namd_sighandler_t)my_sigint_handler);
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);  // after lattice scaling!
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   // if (gotsigint) {
00507   //   iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
00508   //   NAMD_quit();
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);   // step before PME
00547 #endif
00548     }
00549     // signal(SIGINT, oldhandler);
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   // iout << "Controller::minimize() called.\n" << endi;
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;  // 1.0e-3
00602   const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
00603 
00604   CALCULATE
00605 
00606   int minSeq = 0;
00607 
00608   // just move downhill until initial bad contacts go away or 100 steps
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.); // v = f
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     // line minimization
00638     // bracket minimum on line
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     // bracket minimum on line
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.); // v = f
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     // converge on minimum on line
00687     int itcnt;
00688     for ( itcnt = 10; itcnt > 0; --itcnt ) {
00689       int progress = 1;
00690       // select new position
00691       if ( mid.noGradient ) {
00692        if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) {  // subdivide left side
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 {  // subdivide right side
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. ) {  // subdivide left side
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 {  // subdivide right side
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     // new direction
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); // v = c*v+f
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 // NOTE: Only isotropic case implemented here!
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       // Compute new scaling factors and send them to Sequencer
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       // fprintf(stderr, "%d | T %lf P %lf V %1.3lf\n", step, temperature, Pinst*PRESSUREFACTOR, V);
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       // Wait here for Sequencer to finish scaling and force computation
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         // Update temperature for the Temperature Cycle that is coming next
00836         BigReal kineticEnergy = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
00837         temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
00838       }
00839     }
00840 
00841     {
00842       // Update pressure integrator
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       // Wait here for Sequencer to finish scaling and re-calculating kinetic energy
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         // If this is pressure cycle, receive new momentum product
00892         GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
00893       }
00894     }
00895 
00896     // fprintf(stderr, "%d | T %lf scale %lf T' %lf\n", step, T1, v, T2);
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 // Calculates Enthalpy for multigrator
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 //  if (!(step % 100))
00944     // fprintf(stderr, "enthalpy %lf %lf %lf %lf %lf %lf %lf\n", enthalpy,
00945     //   kineticEnergy, potentialEnergy, eta*eta/(2.0*W), P0*V, nuOmegaSum, kT0*(Nf*multigratorZeta[0] + sumZeta));
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     // We only use on-diagonal terms (for now)
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         // We only use on-diagonal terms (for now)
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     // Apply surface tension.  If surfaceTensionTarget is zero, we get
01029     // the default (isotropic pressure) case.
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       // We only use on-diagonal terms (for now)
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 { // langevinPistonBarrier
01071     if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
01072     {
01073       if ( ! positionRescaleFactor.xx ) {  // first pass without MTS
01074       // We only use on-diagonal terms (for now)
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       // We only use on-diagonal terms (for now)
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     // corrections to integrator
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     // corrections to integrator
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     // Apply surface tension.  If surfaceTensionTarget is zero, we get
01179     // the default (isotropic pressure) case.
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         // We only use on-diagonal terms (for now)
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       //iout << "RESCALING VELOCITIES AT STEP " << step
01248       //     << " FROM AVERAGE TEMPERATURE OF " << avgTemp
01249       //     << " TO " << rescaleTemp << " KELVIN.\n" << endi;
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 //Modifications for alchemical fep
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 //fepe
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   {
01352     ++stochRescale_count;
01353     if ( stochRescale_count == simParams->stochRescaleFreq )
01354     { 
01355       const BigReal stochRescaleTemp = simParams->stochRescaleTemp;
01356 
01357       BigReal coefficient = 1.;
01358       if ( temperature > 0. ) 
01359       {
01360         BigReal R1 = random->gaussian();
01361         // BigReal gammaShape = 0.5*(numDegFreedom - 1);
01362         // R2sum is the sum of (numDegFreedom - 1) squared normal variables, which is
01363         // chi-squared distributed. This is in turn a special case of the Gamma
01364         // distribution, which converges to a normal distribution in the limit of a
01365         // large shape parameter.
01366         // BigReal R2sum = 2*(gammaShape + sqrt(gammaShape)*random->gaussian()) + R1*R1;
01367         BigReal R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
01368         BigReal tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
01369 
01370         coefficient = sqrt(stochRescaleTimefactor + (1 - stochRescaleTimefactor)*tempfactor*R2sum
01371                   + 2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*stochRescaleTimefactor)); 
01372       }
01373       broadcast->stochRescaleCoefficient.publish(step,coefficient);
01374       heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient - 1.0);
01375       stochRescale_count = 0;
01376     }
01377   }
01378 }
01379 
01380 static char *FORMAT(BigReal X)
01381 {
01382   static char tmp_string[25];
01383   const double maxnum = 99999999999.9999;
01384   if ( X > maxnum ) X = maxnum;
01385   if ( X < -maxnum ) X = -maxnum;
01386   sprintf(tmp_string," %14.4f",X); 
01387   return tmp_string;
01388 }
01389 
01390 static char *FORMAT(const char *X)
01391 {
01392   static char tmp_string[25];
01393   sprintf(tmp_string," %14s",X); 
01394   return tmp_string;
01395 }
01396 
01397 static char *ETITLE(int X)
01398 {
01399   static char tmp_string[21];
01400   sprintf(tmp_string,"ENERGY: %7d",X); 
01401   return  tmp_string;
01402 }
01403 
01404 void Controller::receivePressure(int step, int minimize)
01405 {
01406     Node *node = Node::Object();
01407     Molecule *molecule = node->molecule;
01408     SimParameters *simParameters = node->simParameters;
01409     Lattice &lattice = state->lattice;
01410 
01411     reduction->require();
01412 
01413     Tensor virial_normal;
01414     Tensor virial_nbond;
01415     Tensor virial_slow;
01416 #ifdef ALTVIRIAL
01417     Tensor altVirial_normal;
01418     Tensor altVirial_nbond;
01419     Tensor altVirial_slow;
01420 #endif
01421     Tensor intVirial_normal;
01422     Tensor intVirial_nbond;
01423     Tensor intVirial_slow;
01424     Vector extForce_normal;
01425     Vector extForce_nbond;
01426     Vector extForce_slow;
01427 
01428 #if 1
01429     numDegFreedom = molecule->num_deg_freedom();
01430     int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
01431     int numFixedGroups = molecule->num_fixed_groups();
01432     int numFixedAtoms = molecule->num_fixed_atoms();
01433 #endif
01434 #if 0
01435     int numAtoms = molecule->numAtoms;
01436     numDegFreedom = 3 * numAtoms;
01437     int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
01438     int numFixedAtoms =
01439         ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
01440     int numLonepairs = molecule->numLonepairs;
01441     int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
01442     if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
01443     if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01444     if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01445     if ( ! ( numFixedAtoms || molecule->numConstraints
01446         || simParameters->comMove || simParameters->langevinOn ) ) {
01447       numDegFreedom -= 3;
01448       numGroupDegFreedom -= 3;
01449     }
01450     if (simParameters->pairInteractionOn) {
01451       // this doesn't attempt to deal with fixed atoms or constraints
01452       numDegFreedom = 3 * molecule->numFepInitial;
01453     }
01454     int numRigidBonds = molecule->numRigidBonds;
01455     int numFixedRigidBonds =
01456         ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01457 
01458     // numLonepairs is subtracted here because all lonepairs have a rigid bond
01459     // to oxygen, but all of the LP degrees of freedom are dealt with above
01460     numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01461 #endif
01462 
01463     kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01464     kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01465 
01466     BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01467         reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01468     BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01469         reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01470 
01471     BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01472                                         / ( numDegFreedom * BOLTZMANN );
01473     BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01474                                         / ( numDegFreedom * BOLTZMANN );
01475     BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01476                                         / ( numGroupDegFreedom * BOLTZMANN );
01477     BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01478                                         / ( numGroupDegFreedom * BOLTZMANN );
01479 
01480     /*  test code for comparing different temperatures  
01481     iout << "TEMPTEST: " << step << " " << 
01482         atomTempHalfstep << " " <<
01483         atomTempCentered << " " <<
01484         groupTempHalfstep << " " <<
01485         groupTempCentered << "\n" << endi;
01486   iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
01487     numGroupDegFreedom << "\n" << endi;
01488      */
01489 
01490     GET_TENSOR(momentumSqrSum,reduction,REDUCTION_MOMENTUM_SQUARED);
01491 
01492     GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01493     GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01494     GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01495 
01496 #ifdef ALTVIRIAL
01497     GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01498     GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01499     GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01500 #endif
01501 
01502     GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01503     GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01504     GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01505 
01506     GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01507     GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01508     GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01509     // APH NOTE: These four lines are now done in calcPressure()
01510     // Vector extPosition = lattice.origin();
01511     // virial_normal -= outer(extForce_normal,extPosition);
01512     // virial_nbond -= outer(extForce_nbond,extPosition);
01513     // virial_slow -= outer(extForce_slow,extPosition);
01514 
01515     kineticEnergy = kineticEnergyCentered;
01516     temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );
01517 
01518     if (simParameters->drudeOn) {
01519       BigReal drudeComKE
01520         = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01521       BigReal drudeBondKE
01522         = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01523       int g_bond = 3 * molecule->numDrudeAtoms;
01524       int g_com = numDegFreedom - g_bond;
01525       temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
01526       drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
01527     }
01528 
01529     // Calculate number of degrees of freedom (controlNumDegFreedom)
01530     if ( simParameters->useGroupPressure )
01531     {
01532       controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01533       if ( ! ( numFixedAtoms || molecule->numConstraints
01534   || simParameters->comMove || simParameters->langevinOn ) ) {
01535         controlNumDegFreedom -= 1;
01536       }
01537     }
01538     else
01539     {
01540       controlNumDegFreedom = numDegFreedom / 3;
01541     }
01542     if (simParameters->fixCellDims) {
01543       if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01544       if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01545       if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01546     }
01547 
01548     // Calculate pressure tensors using virials
01549     calcPressure(step, minimize,
01550       virial_normal, virial_nbond, virial_slow,
01551       intVirial_normal, intVirial_nbond, intVirial_slow,
01552       extForce_normal, extForce_nbond, extForce_slow);
01553 
01554 #ifdef DEBUG_PRESSURE
01555     iout << iINFO << "Control pressure = " << controlPressure <<
01556       " = " << controlPressure_normal << " + " <<
01557       controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01558 #endif
01559     if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
01560         iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
01561              << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
01562              << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
01563              << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
01564              << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
01565              << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
01566              << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
01567              << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
01568              << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
01569              << endi;
01570    }
01571 }
01572 
01573 //
01574 // Calculates all pressure tensors using virials
01575 //
01576 // Sets variables:
01577 // pressure, pressure_normal, pressure_nbond, pressure_slow
01578 // groupPressure, groupPressure_normal, groupPressure_nbond, groupPressure_slow
01579 // controlPressure, controlPressure_normal, controlPressure_nbond, controlPressure_slow
01580 // pressure_amd
01581 // 
01582 void Controller::calcPressure(int step, int minimize,
01583   const Tensor& virial_normal_in, const Tensor& virial_nbond_in, const Tensor& virial_slow_in,
01584   const Tensor& intVirial_normal, const Tensor& intVirial_nbond, const Tensor& intVirial_slow,
01585   const Vector& extForce_normal, const Vector& extForce_nbond, const Vector& extForce_slow) {
01586 
01587   Tensor virial_normal = virial_normal_in;
01588   Tensor virial_nbond = virial_nbond_in;
01589   Tensor virial_slow = virial_slow_in;
01590 
01591   // Tensor tmp = virial_normal;
01592   // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
01593   //   tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
01594 
01595   Node *node = Node::Object();
01596   Molecule *molecule = node->molecule;
01597   SimParameters *simParameters = node->simParameters;
01598   Lattice &lattice = state->lattice;
01599 
01600   BigReal volume;
01601 
01602   Vector extPosition = lattice.origin();
01603   virial_normal -= outer(extForce_normal,extPosition);
01604   virial_nbond -= outer(extForce_nbond,extPosition);
01605   virial_slow -= outer(extForce_slow,extPosition);
01606 
01607   if ( (volume=lattice.volume()) != 0. )
01608   {
01609 
01610     if (simParameters->LJcorrection && volume) {
01611 #ifdef MEM_OPT_VERSION
01612       NAMD_bug("LJcorrection not supported in memory optimized build.");
01613 #else
01614       // Apply tail correction to pressure
01615       BigReal alchLambda = simParameters->getCurrentLambda(step);
01616       virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
01617 #endif
01618     }
01619 
01620     // kinetic energy component included in virials
01621     pressure_normal = virial_normal / volume;
01622     groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01623 
01624     if (simParameters->accelMDOn) {
01625       pressure_amd = virial_amd / volume;
01626       pressure_normal += pressure_amd;
01627       groupPressure_normal +=  pressure_amd;
01628     }
01629 
01630     if ( minimize || ! ( step % nbondFreq ) )
01631     {
01632       pressure_nbond = virial_nbond / volume;
01633       groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01634     }
01635 
01636     if ( minimize || ! ( step % slowFreq ) )
01637     {
01638       pressure_slow = virial_slow / volume;
01639       groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01640     }
01641 
01642     pressure = pressure_normal + pressure_nbond + pressure_slow; 
01643     groupPressure = groupPressure_normal + groupPressure_nbond +
01644           groupPressure_slow;
01645   }
01646   else
01647   {
01648     pressure = Tensor();
01649     groupPressure = Tensor();
01650   }
01651 
01652   if ( simParameters->useGroupPressure )
01653   {
01654     controlPressure_normal = groupPressure_normal;
01655     controlPressure_nbond = groupPressure_nbond;
01656     controlPressure_slow = groupPressure_slow;
01657     controlPressure = groupPressure;
01658   }
01659   else
01660   {
01661     controlPressure_normal = pressure_normal;
01662     controlPressure_nbond = pressure_nbond;
01663     controlPressure_slow = pressure_slow;
01664     controlPressure = pressure;
01665   }
01666 
01667   if ( simParameters->useFlexibleCell ) {
01668     // use symmetric pressure to control rotation
01669     // controlPressure_normal = symmetric(controlPressure_normal);
01670     // controlPressure_nbond = symmetric(controlPressure_nbond);
01671     // controlPressure_slow = symmetric(controlPressure_slow);
01672     // controlPressure = symmetric(controlPressure);
01673     // only use on-diagonal components for now
01674     controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01675     controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01676     controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01677     controlPressure = Tensor::diagonal(diagonal(controlPressure));
01678     if ( simParameters->useConstantRatio ) {
01679 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01680    T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01681       AVGXY(controlPressure_normal);
01682       AVGXY(controlPressure_nbond);
01683       AVGXY(controlPressure_slow);
01684       AVGXY(controlPressure);
01685 #undef AVGXY
01686     }
01687   } else {
01688     controlPressure_normal =
01689   Tensor::identity(trace(controlPressure_normal)/3.);
01690     controlPressure_nbond =
01691   Tensor::identity(trace(controlPressure_nbond)/3.);
01692     controlPressure_slow =
01693   Tensor::identity(trace(controlPressure_slow)/3.);
01694     controlPressure =
01695   Tensor::identity(trace(controlPressure)/3.);
01696   }
01697 }
01698 
01699 inline void Controller :: calc_accelMDG_mean_std
01700 (BigReal testV, int n, 
01701  BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV){
01702    BigReal delta; 
01703 
01704    if(testV > *Vmax){
01705       *Vmax = testV;
01706    }
01707    else if(testV < *Vmin){
01708       *Vmin = testV;
01709    }
01710 
01711    //mean and std calculated by Online Algorithm
01712    delta = testV - *Vavg;
01713    *Vavg += delta / (BigReal)n;
01714    *M2 += delta * (testV - (*Vavg));
01715 
01716    *sigmaV = sqrt(*M2/n);
01717 }
01718 
01719 inline void Controller :: calc_accelMDG_E_k
01720 (int iE, int step, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, 
01721  BigReal* k0, BigReal* k, BigReal* E, int *iEused, char *warn){
01722    switch(iE){
01723       case 2:
01724          *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
01725          // if k0 <=0 OR k0 > 1, switch to iE=1 mode
01726          if(*k0 > 1.)
01727             *warn |= 1;
01728          else if(*k0 <= 0.)
01729             *warn |= 2;
01730          //else stay in iE=2 mode
01731          else{
01732             *E = Vmin + (Vmax-Vmin)/(*k0);
01733             *iEused = 2;
01734             break;
01735          }
01736       case 1:
01737          *E = Vmax;
01738          *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
01739          if(*k0 > 1.)
01740             *k0 = 1.;
01741          *iEused = 1;
01742          break;
01743    }
01744 
01745    *k = *k0 / (Vmax - Vmin);
01746 }
01747 
01748 inline void Controller :: calc_accelMDG_force_factor
01749 (BigReal k, BigReal E, BigReal testV, Tensor vir_orig, 
01750  BigReal *dV, BigReal *factor, Tensor *vir){
01751    BigReal VE = testV - E;
01752    //if V < E, apply boost
01753    if(VE < 0){
01754       *factor = k * VE;
01755       *vir += vir_orig * (*factor);
01756       *dV += (*factor) * VE/2.;
01757       *factor += 1.;
01758    }
01759    else{
01760       *factor = 1.;
01761    }
01762 }
01763 
01764 void Controller :: write_accelMDG_rest_file
01765 (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){
01766    FILE *rest;
01767    char timestepstr[20];
01768    char *filename;
01769    int baselen;
01770    char *restart_name;
01771    const char *bsuffix;
01772 
01773    if(lasttime || simParams->restartFrequency == 0){
01774            filename = simParams->outputFilename;
01775            bsuffix = ".BAK";
01776    }
01777    else{
01778            filename = simParams->restartFilename;
01779            bsuffix = ".old";
01780    }
01781 
01782    baselen = strlen(filename);
01783    restart_name = new char[baselen+26];
01784 
01785    strcpy(restart_name, filename);
01786    if ( simParams->restartSave && !lasttime) {
01787       sprintf(timestepstr,".%d",step_n);
01788       strcat(restart_name, timestepstr);
01789       bsuffix = ".BAK";
01790    }
01791    strcat(restart_name, ".gamd");
01792 
01793    if(write_topic){
01794       NAMD_backup_file(restart_name,bsuffix);
01795 
01796       rest = fopen(restart_name, "w");
01797       if(rest){
01798          fprintf(rest, "# NAMD accelMDG restart file\n"
01799                "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
01800                "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01801                type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01802          fclose(rest);
01803          iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
01804       }
01805       else
01806          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01807    }
01808    else{
01809       rest = fopen(restart_name, "a");
01810       if(rest){
01811          fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01812                           type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01813          fclose(rest);
01814       }
01815       else
01816          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01817    }
01818 
01819    delete[] restart_name;
01820 }
01821 
01822 
01823 void Controller::rescaleaccelMD(int step, int minimize)
01824 {
01825     if ( !simParams->accelMDOn ) return;
01826 
01827     amd_reduction->require();
01828 
01829     // copy all to submit_reduction
01830     for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
01831       submit_reduction->item(i) += amd_reduction->item(i);
01832     }
01833     submit_reduction->submit();
01834 
01835     if (step == simParams->firstTimestep) accelMDdVAverage = 0;
01836 //    if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
01837     if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
01838 
01839     Node *node = Node::Object();
01840     Molecule *molecule = node->molecule;
01841     Lattice &lattice = state->lattice;
01842 
01843     const BigReal accelMDE = simParams->accelMDE;
01844     const BigReal accelMDalpha = simParams->accelMDalpha;
01845     const BigReal accelMDTE = simParams->accelMDTE;
01846     const BigReal accelMDTalpha = simParams->accelMDTalpha;
01847     const int accelMDOutFreq = simParams->accelMDOutFreq;
01848 
01849     //GaMD
01850     static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
01851     static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
01852     static BigReal k0P, kP, EP;
01853     static BigReal k0D, kD, ED;
01854     static int V_n = 1, iEusedD, iEusedP;
01855     static char warnD, warnP;
01856     static int restfreq;
01857 
01858     const int ntcmdprep = simParams->accelMDGcMDPrepSteps;
01859     const int ntcmd     = simParams->accelMDGcMDSteps;
01860     const int ntebprep  = ntcmd + simParams->accelMDGEquiPrepSteps;
01861     const int nteb      = ntcmd + simParams->accelMDGEquiSteps;
01862 
01863     BigReal volume;
01864     BigReal bondEnergy;
01865     BigReal angleEnergy;
01866     BigReal dihedralEnergy;
01867     BigReal improperEnergy;
01868     BigReal crosstermEnergy;
01869     BigReal boundaryEnergy;
01870     BigReal miscEnergy;
01871     BigReal amd_electEnergy;
01872     BigReal amd_ljEnergy;
01873     BigReal amd_ljEnergy_Corr = 0.;
01874     BigReal amd_electEnergySlow;
01875     BigReal amd_groLJEnergy;
01876     BigReal amd_groGaussEnergy;
01877     BigReal amd_goTotalEnergy;
01878     BigReal potentialEnergy;
01879     BigReal factor_dihe = 1;
01880     BigReal factor_tot = 1;
01881     BigReal testV;
01882     BigReal dV = 0.;
01883     Vector  accelMDfactor;
01884     Tensor vir; //auto initialized to 0
01885     Tensor vir_dihe;
01886     Tensor vir_normal;
01887     Tensor vir_nbond;
01888     Tensor vir_slow;
01889 
01890     volume = lattice.volume();
01891 
01892     bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
01893     angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
01894     dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
01895     improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
01896     crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
01897     boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
01898     miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);
01899 
01900     GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
01901     GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
01902     GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
01903     GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
01904 
01905     if ( !( step % nbondFreq ) ) {
01906       amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
01907       amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
01908       amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
01909       amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
01910       BigReal goNativeEnergy = amd_reduction->item(REDUCTION_GO_NATIVE_ENERGY);
01911       BigReal goNonnativeEnergy = amd_reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
01912       amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
01913     } else {
01914       amd_electEnergy = electEnergy;
01915       amd_ljEnergy = ljEnergy;
01916       amd_groLJEnergy = groLJEnergy;
01917       amd_groGaussEnergy = groGaussEnergy;
01918       amd_goTotalEnergy = goTotalEnergy;
01919     }
01920 
01921     if( simParams->LJcorrection && volume ) {
01922         // Obtain tail correction (copied from printEnergies())
01923         // This value is only printed for sanity check
01924         // accelMD currently does not 'boost' tail correction
01925         amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
01926     }
01927 
01928     if ( !( step % slowFreq ) ) {
01929       amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
01930     } else {
01931       amd_electEnergySlow = electEnergySlow;
01932     }
01933 
01934     potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
01935         improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
01936         crosstermEnergy + boundaryEnergy + miscEnergy +
01937         amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
01938 
01939     //GaMD
01940     if(simParams->accelMDG){
01941        //fix step when there is accelMDFirstStep
01942        step -= simParams->accelMDFirstStep;
01943 
01944        if(step == simParams->firstTimestep) {
01945           iEusedD = iEusedP = simParams->accelMDGiE;
01946           warnD   = warnP   = '\0';
01947 
01948           //restart file reading
01949           if(simParams->accelMDGRestart){
01950              FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
01951              char line[256];
01952              int dihe_n=0, tot_n=0;
01953              if(!rest){
01954                 sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
01955                 NAMD_die(line);
01956              }
01957 
01958              while(fgets(line, 256, rest)){
01959                 if(line[0] == 'D'){
01960                    dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
01961                                    &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
01962                 }
01963                 else if(line[0] == 'T'){
01964                    tot_n  = sscanf(line+1, " %d %la %la %la %la %la %la %la",
01965                                    &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
01966                 }
01967              }
01968 
01969              fclose(rest);
01970 
01971              V_n++;
01972 
01973              //check dihe settings
01974              if(simParams->accelMDdihe || simParams->accelMDdual){
01975                 if(dihe_n !=8)
01976                    NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
01977                 k0D = kD * (VmaxD - VminD);
01978                 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
01979                    << " Vmax " << VmaxD << " Vmin " << VminD 
01980                    << " Vavg " << VavgD << " sigmaV " << sigmaVD
01981                    << " E " << ED << " k " << kD
01982                    << "\n" << endi;
01983              }
01984 
01985              //check tot settings
01986              if(!simParams->accelMDdihe || simParams->accelMDdual){
01987                 if(tot_n !=8)
01988                    NAMD_die("Invalid total potential energy format in accelMDG restart file");
01989                 k0P = kP * (VmaxP - VminP);
01990                 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
01991                    << " Vmax " << VmaxP << " Vmin " << VminP 
01992                    << " Vavg " << VavgP << " sigmaV " << sigmaVP
01993                    << " E " << EP << " k " << kP
01994                    << "\n" << endi;
01995             }
01996 
01997             iEusedD = (ED == VmaxD) ? 1 : 2;
01998             iEusedP = (EP == VmaxP) ? 1 : 2;
01999           }
02000           //local restfreq follows NAMD restartfreq (default: 500)
02001           restfreq = simParams->restartFrequency ? simParams->restartFrequency : 500;
02002        }
02003 
02004        //for dihedral potential
02005        if(simParams->accelMDdihe || simParams->accelMDdual){
02006           testV = dihedralEnergy + crosstermEnergy;
02007 
02008           //write restart file every restartfreq steps
02009           if(step > simParams->firstTimestep - simParams->accelMDFirstStep && step % restfreq == 0)
02010              write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
02011                    true, false);
02012           //write restart file at the end of the simulation
02013           if( simParams->accelMDLastStep > 0 ){
02014               if( step == simParams->accelMDLastStep )
02015                   write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
02016                           true, true);
02017           }
02018           else if(step == simParams->N - simParams->accelMDFirstStep)
02019               write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
02020                       true, true);
02021 
02022           //conventional MD
02023           if(step < ntcmd){
02024              //very first step
02025              if(step == 0 && !simParams->accelMDGRestart){
02026                 //initialize stat
02027                 VmaxD = VminD = VavgD = testV;
02028                 M2D = sigmaVD = 0.;
02029              }
02030              //first step after cmdprep
02031              else if(step == ntcmdprep && ntcmdprep != 0){
02032                 //reset stat
02033                 VmaxD = VminD = VavgD = testV;
02034                 M2D = sigmaVD = 0.;
02035                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
02036              }
02037              //normal steps
02038              else
02039                 calc_accelMDG_mean_std(testV, V_n,
02040                       &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
02041 
02042              //last cmd step
02043              if(step == ntcmd - 1){
02044                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D, 
02045                       VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
02046              }
02047           }
02048           //equilibration
02049           else if(step < nteb){
02050              calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
02051                    &dV, &factor_dihe, &vir);
02052 
02053              //first step after cmd
02054              if(step == ntcmd && simParams->accelMDGresetVaftercmd){
02055                 //reset stat
02056                 VmaxD = VminD = VavgD = testV;
02057                 M2D = sigmaVD = 0.;
02058                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
02059              }
02060              else
02061                 calc_accelMDG_mean_std(testV, V_n, 
02062                       &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
02063 
02064              //steps after ebprep
02065              if(step >= ntebprep)
02066                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D, 
02067                       VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
02068           }
02069           //production
02070           else{
02071              calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
02072                    &dV, &factor_dihe, &vir);
02073           }
02074 
02075        }
02076        //for total potential
02077        if(!simParams->accelMDdihe || simParams->accelMDdual){
02078           Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
02079           testV = potentialEnergy;
02080           if(simParams->accelMDdual){
02081              testV -= dihedralEnergy + crosstermEnergy;
02082              vir_tot -= vir_dihe;
02083           }
02084 
02085           //write restart file every restartfreq steps
02086           if(step > simParams->firstTimestep - simParams->accelMDFirstStep && step % restfreq == 0)
02087              write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
02088                    !simParams->accelMDdual, false);
02089           //write restart file at the end of simulation
02090           if( simParams->accelMDLastStep > 0 ){
02091               if( step == simParams->accelMDLastStep )
02092                   write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
02093                           !simParams->accelMDdual, true);
02094           }
02095           else if(step == simParams->N - simParams->accelMDFirstStep)
02096              write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
02097                    !simParams->accelMDdual, true);
02098 
02099           //conventional MD
02100           if(step < ntcmd){
02101              //very first step
02102              if(step == 0 && !simParams->accelMDGRestart){
02103                 //initialize stat
02104                 VmaxP = VminP = VavgP = testV;
02105                 M2P = sigmaVP = 0.;
02106              }
02107              //first step after cmdprep
02108              else if(step == ntcmdprep && ntcmdprep != 0){
02109                 //reset stat
02110                 VmaxP = VminP = VavgP = testV;
02111                 M2P = sigmaVP = 0.;
02112                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
02113              }
02114              //normal steps
02115              else
02116                 calc_accelMDG_mean_std(testV, V_n,
02117                       &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
02118              //last cmd step
02119              if(step == ntcmd - 1){
02120                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P, 
02121                       VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
02122              }
02123           }
02124           //equilibration
02125           else if(step < nteb){
02126              calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
02127                    &dV, &factor_tot, &vir);
02128 
02129              //first step after cmd
02130              if(step == ntcmd && simParams->accelMDGresetVaftercmd){
02131                 //reset stat
02132                 VmaxP = VminP = VavgP = testV;
02133                 M2P = sigmaVP = 0.;
02134                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
02135              }
02136              else
02137                 calc_accelMDG_mean_std(testV, V_n,
02138                       &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
02139 
02140              //steps after ebprep
02141              if(step >= ntebprep)
02142                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P, 
02143                       VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
02144           }
02145           //production
02146           else{
02147              calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
02148                    &dV, &factor_tot, &vir);
02149           }
02150 
02151        }
02152        accelMDdVAverage += dV;
02153 
02154        //first step after ntcmdprep
02155        if((ntcmdprep > 0 && step == ntcmdprep) || 
02156           (simParams->accelMDGresetVaftercmd && step == ntcmd)){
02157           V_n = 1;
02158        }
02159 
02160        if(step < nteb)
02161            V_n++;
02162 
02163        // restore step
02164        step += simParams->accelMDFirstStep;
02165     }
02166     //aMD
02167     else{
02168         if (simParams->accelMDdihe) {
02169 
02170             testV = dihedralEnergy + crosstermEnergy;
02171             if ( testV < accelMDE ) {
02172                 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
02173                 factor_dihe *= factor_dihe;
02174                 vir = vir_dihe * (factor_dihe - 1.0);
02175                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02176                 accelMDdVAverage += dV;
02177             }  
02178 
02179         } else if (simParams->accelMDdual) {
02180 
02181             testV = dihedralEnergy + crosstermEnergy;
02182             if ( testV < accelMDE ) {
02183                 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
02184                 factor_dihe *= factor_dihe;
02185                 vir = vir_dihe * (factor_dihe - 1.0) ;
02186                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02187             }
02188 
02189             testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
02190             if ( testV < accelMDTE ) {
02191                 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
02192                 factor_tot *= factor_tot;
02193                 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
02194                 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
02195             }
02196             accelMDdVAverage += dV;
02197 
02198         } else {
02199 
02200             testV = potentialEnergy;
02201             if ( testV < accelMDE ) {
02202                 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
02203                 factor_tot *= factor_tot;
02204                 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
02205                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02206                 accelMDdVAverage += dV;
02207             }
02208         } 
02209     }
02210 
02211     accelMDfactor[0]=factor_dihe;
02212     accelMDfactor[1]=factor_tot;
02213     accelMDfactor[2]=1;
02214     broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
02215     virial_amd = vir; 
02216 
02217     if ( factor_tot < 0.001 ) {
02218        iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
02219             << "\n" << endi;
02220     }
02221 
02222     if ( ! (step % accelMDOutFreq) ) {
02223         if ( !(step == simParams->firstTimestep) ) {
02224              accelMDdVAverage = accelMDdVAverage/accelMDOutFreq; 
02225         }
02226         iout << "ACCELERATED MD: STEP " << step
02227              << " dV "   << dV
02228              << " dVAVG " << accelMDdVAverage 
02229              << " BOND " << bondEnergy
02230              << " ANGLE " << angleEnergy
02231              << " DIHED " << dihedralEnergy+crosstermEnergy
02232              << " IMPRP " << improperEnergy
02233              << " ELECT " << amd_electEnergy+amd_electEnergySlow
02234              << " VDW "  << amd_ljEnergy
02235              << " POTENTIAL "  << potentialEnergy;
02236         if(amd_ljEnergy_Corr)
02237             iout << " LJcorr " << amd_ljEnergy_Corr;
02238         iout << "\n" << endi;
02239         if(simParams->accelMDG){
02240             if(simParams->accelMDdihe || simParams->accelMDdual)
02241                 iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD 
02242                     << " Vmax " << VmaxD << " Vmin " << VminD 
02243                     << " Vavg " << VavgD << " sigmaV " << sigmaVD
02244                     << " E " << ED << " k0 " << k0D << " k " << kD
02245                     << "\n" << endi;
02246             if(warnD & 1)
02247                 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
02248                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02249             if(warnD & 2)
02250                 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
02251                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02252             if(!simParams->accelMDdihe || simParams->accelMDdual)
02253                 iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP 
02254                     << " Vmax " << VmaxP << " Vmin " << VminP 
02255                     << " Vavg " << VavgP << " sigmaV " << sigmaVP
02256                     << " E " << EP << " k0 " << k0P << " k " << kP
02257                     << "\n" << endi;
02258             if(warnP & 1)
02259                 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
02260                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02261             if(warnP & 2)
02262                 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
02263                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02264             warnD = warnP = '\0';
02265         }
02266 
02267         accelMDdVAverage = 0;
02268 
02269         if (simParams->accelMDDebugOn) {
02270            Tensor p_normal;
02271            Tensor p_nbond;
02272            Tensor p_slow;
02273            Tensor p;
02274            if ( volume != 0. )  {
02275                  p_normal = vir_normal/volume;
02276                  p_nbond  = vir_nbond/volume;
02277                  p_slow = vir_slow/volume;
02278                  p = vir/volume;
02279            }    
02280            iout << " accelMD Scaling Factor: " << accelMDfactor << "\n" 
02281                 << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
02282                 << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
02283                 << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
02284                 << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n" 
02285                 << endi;
02286         }
02287    }
02288 }
02289 
02290 void Controller::adaptTempInit(int step) {
02291     if (!simParams->adaptTempOn) return;
02292     iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
02293     adaptTempDtMin = 0;
02294     adaptTempDtMax = 0;
02295     adaptTempAutoDt = false;
02296     if (simParams->adaptTempBins == 0) {
02297       iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
02298       std::ifstream adaptTempRead(simParams->adaptTempInFile);
02299       if (adaptTempRead) {
02300         int readInt;
02301         BigReal readReal;
02302         bool readBool;
02303         // step
02304         adaptTempRead >> readInt;
02305         // Start with min and max temperatures
02306         adaptTempRead >> adaptTempT;     // KELVIN
02307         adaptTempRead >> adaptTempBetaMin;  // KELVIN
02308         adaptTempRead >> adaptTempBetaMax;  // KELVIN
02309         adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
02310         adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
02311         // In case file is manually edited
02312         if (adaptTempBetaMin > adaptTempBetaMax){
02313             readReal = adaptTempBetaMax;
02314             adaptTempBetaMax = adaptTempBetaMin;
02315             adaptTempBetaMin = adaptTempBetaMax;
02316         }
02317         adaptTempRead >> adaptTempBins;     
02318         adaptTempRead >> adaptTempCg;
02319         adaptTempRead >> adaptTempDt;
02320         adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02321         adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02322         adaptTempPotEnergySamples = new int[adaptTempBins];
02323         adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02324         adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02325         adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02326         adaptTempBetaN           = new BigReal[adaptTempBins];
02327         adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02328         for(int j = 0; j < adaptTempBins; ++j) {
02329           adaptTempRead >> adaptTempPotEnergyAve[j];
02330           adaptTempRead >> adaptTempPotEnergyVar[j];
02331           adaptTempRead >> adaptTempPotEnergySamples[j];
02332           adaptTempRead >> adaptTempPotEnergyAveNum[j];
02333           adaptTempRead >> adaptTempPotEnergyVarNum[j];
02334           adaptTempRead >> adaptTempPotEnergyAveDen[j];
02335           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02336         } 
02337         adaptTempRead.close();
02338       }
02339       else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
02340     } 
02341     else {
02342       adaptTempBins = simParams->adaptTempBins;
02343       adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02344       adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02345       adaptTempPotEnergySamples = new int[adaptTempBins];
02346       adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02347       adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02348       adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02349       adaptTempBetaN           = new BigReal[adaptTempBins];
02350       adaptTempBetaMax = 1./simParams->adaptTempTmin;
02351       adaptTempBetaMin = 1./simParams->adaptTempTmax;
02352       adaptTempCg = simParams->adaptTempCgamma;   
02353       adaptTempDt  = simParams->adaptTempDt;
02354       adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02355       adaptTempT = simParams->initialTemp; 
02356       if (simParams->langevinOn)
02357         adaptTempT = simParams->langevinTemp;
02358       else if (simParams->rescaleFreq > 0)
02359         adaptTempT = simParams->rescaleTemp;
02360       for(int j = 0; j < adaptTempBins; ++j){
02361           adaptTempPotEnergyAveNum[j] = 0.;
02362           adaptTempPotEnergyAveDen[j] = 0.;
02363           adaptTempPotEnergySamples[j] = 0;
02364           adaptTempPotEnergyVarNum[j] = 0.;
02365           adaptTempPotEnergyVar[j] = 0.;
02366           adaptTempPotEnergyAve[j] = 0.;
02367           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02368       }
02369     }
02370     if (simParams->adaptTempAutoDt > 0.0) {
02371        adaptTempAutoDt = true;
02372        adaptTempDtMin =  simParams->adaptTempAutoDt - 0.01;
02373        adaptTempDtMax =  simParams->adaptTempAutoDt + 0.01;
02374     }
02375     adaptTempDTave = 0;
02376     adaptTempDTavenum = 0;
02377     iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
02378     iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
02379     iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
02380     iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
02381     iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
02382     if (simParams->adaptTempRestartFreq > 0) {
02383       NAMD_backup_file(simParams->adaptTempRestartFile);
02384       adaptTempRestartFile.open(simParams->adaptTempRestartFile);
02385        if(!adaptTempRestartFile)
02386         NAMD_die("Error opening restart file");
02387     }
02388 }
02389 
02390 void Controller::adaptTempWriteRestart(int step) {
02391     if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
02392         adaptTempRestartFile.seekp(std::ios::beg);        
02393         iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
02394         adaptTempRestartFile << step << " ";
02395         // Start with min and max temperatures
02396         adaptTempRestartFile << adaptTempT<< " ";     // KELVIN
02397         adaptTempRestartFile << 1./adaptTempBetaMin << " ";  // KELVIN
02398         adaptTempRestartFile << 1./adaptTempBetaMax << " ";  // KELVIN
02399         adaptTempRestartFile << adaptTempBins << " ";     
02400         adaptTempRestartFile << adaptTempCg << " ";
02401         adaptTempRestartFile << adaptTempDt ;
02402         adaptTempRestartFile << "\n" ;
02403         for(int j = 0; j < adaptTempBins; ++j) {
02404           adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
02405           adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
02406           adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
02407           adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
02408           adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
02409           adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
02410           adaptTempRestartFile << "\n";          
02411         }
02412         adaptTempRestartFile.flush(); 
02413     }
02414 }    
02415 
02416 void Controller::adaptTempUpdate(int step, int minimize)
02417 {
02418     //Beta = 1./T
02419     if ( !simParams->adaptTempOn ) return;
02420     int j = 0;
02421     if (step == simParams->firstTimestep) {
02422         adaptTempInit(step);
02423         return;
02424     }
02425     if ( minimize || (step < simParams->adaptTempFirstStep ) || 
02426         ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
02427     const int adaptTempOutFreq  = simParams->adaptTempOutFreq;
02428     const bool adaptTempDebug  = simParams->adaptTempDebug;
02429     //Calculate Current inverse temperature and bin 
02430     BigReal adaptTempBeta = 1./adaptTempT;
02431     adaptTempBin   = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
02432 
02433     if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
02434         iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin  
02435                                << " adaptTempBeta: " << adaptTempBeta 
02436                               << " adaptTempDBeta: " << adaptTempDBeta 
02437                                << " betaMin:" << adaptTempBetaMin 
02438                                << " betaMax: " << adaptTempBetaMax << "\n";
02439     adaptTempPotEnergySamples[adaptTempBin] += 1;
02440     BigReal gammaAve = 1.-adaptTempCg/adaptTempPotEnergySamples[adaptTempBin];
02441 
02442     BigReal potentialEnergy;
02443     BigReal potEnergyAverage;
02444     BigReal potEnergyVariance;
02445     potentialEnergy = totalEnergy-kineticEnergy;
02446 
02447     //calculate new bin average and variance using adaptive averaging
02448     adaptTempPotEnergyAveNum[adaptTempBin] = adaptTempPotEnergyAveNum[adaptTempBin]*gammaAve + potentialEnergy;
02449     adaptTempPotEnergyAveDen[adaptTempBin] = adaptTempPotEnergyAveDen[adaptTempBin]*gammaAve + 1;
02450     adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
02451     
02452     potEnergyAverage = adaptTempPotEnergyAveNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
02453     potEnergyVariance = adaptTempPotEnergyVarNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
02454     potEnergyVariance -= potEnergyAverage*potEnergyAverage;
02455 
02456     adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
02457     adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
02458     
02459     // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
02460     // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
02461     // is constant for each bin. This is to estimate <E(beta)> where beta \in
02462     // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
02463     if ( ! ( step % simParams->adaptTempFreq ) ) {
02464       // If adaptTempBin not at the edge, calculate improved average:
02465       if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
02466           // Get Averaging Limits:
02467           BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
02468           BigReal betaPlus;
02469           BigReal betaMinus;
02470           int     nMinus =0;
02471           int     nPlus = 0;
02472           if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
02473           if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
02474           betaMinus = adaptTempBeta - deltaBeta;
02475           betaPlus =  adaptTempBeta + deltaBeta;
02476           BigReal betaMinus2 = betaMinus*betaMinus;
02477           BigReal betaPlus2  = betaPlus*betaPlus;
02478           nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
02479           nPlus  = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
02480           // Variables for <E(beta)> estimate:
02481           BigReal potEnergyAve0 = 0.0;
02482           BigReal potEnergyAve1 = 0.0;
02483           // Integral terms
02484           BigReal A0 = 0;
02485           BigReal A1 = 0;
02486           BigReal A2 = 0;
02487           //A0 phi_s integral for beta_minus < beta < beta_{i+1}
02488           BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1]; 
02489           BigReal DeltaE2Ave;
02490           BigReal DeltaE2AveSum = 0;
02491           BigReal betaj;
02492           for (j = nMinus+1; j <= adaptTempBin; ++j) {
02493               potEnergyAve0 += adaptTempPotEnergyAve[j]; 
02494               DeltaE2Ave = adaptTempPotEnergyVar[j];
02495               DeltaE2AveSum += DeltaE2Ave;
02496               A0 += j*DeltaE2Ave;
02497           }
02498           A0 *= adaptTempDBeta;
02499           A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
02500           A0 *= adaptTempDBeta;
02501           betaj = adaptTempBetaN[nMinus+1]-betaMinus; 
02502           betaj *= betaj;
02503           A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
02504           A0 /= (betaNp1-betaMinus);
02505 
02506           //A1 phi_s integral for beta_{i+1} < beta < beta_plus
02507           DeltaE2AveSum = 0;
02508           for (j = adaptTempBin+1; j < nPlus; j++) {
02509               potEnergyAve1 += adaptTempPotEnergyAve[j];
02510               DeltaE2Ave = adaptTempPotEnergyVar[j];
02511               DeltaE2AveSum += DeltaE2Ave;
02512               A1 += j*DeltaE2Ave;
02513           }
02514           A1 *= adaptTempDBeta;   
02515           A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
02516           A1 *= adaptTempDBeta;
02517           if ( nPlus < adaptTempBins ) {
02518             betaj = betaPlus-adaptTempBetaN[nPlus];
02519             betaj *= betaj;
02520             A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
02521           }
02522           A1 /= (betaPlus-betaNp1);
02523 
02524           //A2 phi_t integral for beta_i
02525           A2 = 0.5*adaptTempDBeta*potEnergyVariance;
02526 
02527           // Now calculate a+ and a-
02528           DeltaE2Ave = A0-A1;
02529           BigReal aplus = 0;
02530           BigReal aminus = 0;
02531           if (DeltaE2Ave != 0) {
02532             aplus = (A0-A2)/(A0-A1);
02533             if (aplus < 0) {
02534                     aplus = 0;
02535             }
02536             if (aplus > 1)  {
02537                     aplus = 1;
02538             }
02539             aminus = 1-aplus;
02540             potEnergyAve0 *= adaptTempDBeta;
02541             potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
02542             potEnergyAve0 /= (betaNp1-betaMinus);
02543             potEnergyAve1 *= adaptTempDBeta;
02544             if ( nPlus < adaptTempBins ) {
02545                 potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
02546             }
02547             potEnergyAve1 /= (betaPlus-betaNp1);
02548             potEnergyAverage = aminus*potEnergyAve0;
02549             potEnergyAverage += aplus*potEnergyAve1;
02550           }
02551           if (simParams->adaptTempDebug) {
02552        iout << "ADAPTEMP DEBUG:"  << "\n"
02553             << "     adaptTempBin:    " << adaptTempBin << "\n"
02554             << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
02555             << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
02556             << "     adaptTemp:   " << adaptTempT<< "\n"
02557             << "     betaMin:   " << adaptTempBetaMin << "\n"
02558             << "     betaMax:   " << adaptTempBetaMax << "\n"
02559             << "     gammaAve:  " << gammaAve << "\n"
02560             << "     deltaBeta: " << deltaBeta << "\n"
02561             << "     betaMinus: " << betaMinus << "\n"
02562             << "     betaPlus:  " << betaPlus << "\n"
02563             << "     nMinus:    " << nMinus << "\n"
02564             << "     nPlus:     " << nPlus << "\n"
02565             << "     A0:        " << A0 << "\n"
02566             << "     A1:        " << A1 << "\n"
02567             << "     A2:        " << A2 << "\n"
02568             << "     a+:        " << aplus << "\n"
02569             << "     a-:        " << aminus << "\n"
02570             << endi;
02571           }
02572       }
02573       else {
02574           if (simParams->adaptTempDebug) {
02575        iout << "ADAPTEMP DEBUG:"  << "\n"
02576             << "     adaptTempBin:    " << adaptTempBin << "\n"
02577             << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
02578             << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
02579             << "     adaptTemp:   " << adaptTempT<< "\n"
02580             << "     betaMin:   " << adaptTempBetaMin << "\n"
02581             << "     betaMax:   " << adaptTempBetaMax << "\n"
02582             << "     gammaAve:  " << gammaAve << "\n"
02583             << "     deltaBeta:  N/A\n"
02584             << "     betaMinus:  N/A\n"
02585             << "     betaPlus:   N/A\n"
02586             << "     nMinus:     N/A\n"
02587             << "     nPlus:      N/A\n"
02588             << "     A0:         N/A\n"
02589             << "     A1:         N/A\n"
02590             << "     A2:         N/A\n"
02591             << "     a+:         N/A\n"
02592             << "     a-:         N/A\n"
02593             << endi;
02594           }
02595       }
02596       
02597       //dT is new temperature
02598       BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
02599       dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
02600       dT += adaptTempT;
02601       
02602      // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
02603      // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
02604       if ( dT > 1./adaptTempBetaMin || dT  < 1./adaptTempBetaMax ) {
02605         dT = ((potentialEnergy-adaptTempPotEnergyAve[adaptTempBin])/BOLTZMANN+adaptTempT)*adaptTempDt;
02606         dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
02607         dT += adaptTempT;
02608         // Check again, if not then keep original adaptTempTor assign random.
02609         if ( dT > 1./adaptTempBetaMin ) {
02610           if (!simParams->adaptTempRandom) {             
02611              //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
02612              //     << " K higher than adaptTempTmax."
02613              //     << " Keeping temperature at " 
02614              //     << adaptTempT<< "\n"<< endi;             
02615              dT = adaptTempT;
02616           }
02617           else {             
02618              //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
02619              //     << " K higher than adaptTempTmax."
02620              //     << " Assigning random temperature in range\n" << endi;
02621              dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);             
02622              dT = 1./dT;
02623           }
02624         } 
02625         else if ( dT  < 1./adaptTempBetaMax ) {
02626           if (!simParams->adaptTempRandom) {            
02627             //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
02628             //     << " K lower than adaptTempTmin."
02629             //     << " Keeping temperature at " << adaptTempT<< "\n" << endi; 
02630             dT = adaptTempT;
02631           }
02632           else {
02633             //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
02634             //     << " K lower than adaptTempTmin."
02635             //     << " Assigning random temperature in range\n" << endi;
02636             dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
02637             dT = 1./dT;
02638           }
02639         }
02640         else if (adaptTempAutoDt) {
02641           //update temperature step size counter
02642           //FOR "TRUE" ADAPTIVE TEMPERING 
02643           BigReal adaptTempTdiff = fabs(dT-adaptTempT);
02644           if (adaptTempTdiff > 0) {
02645             adaptTempDTave += adaptTempTdiff; 
02646             adaptTempDTavenum++;
02647 //            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
02648           }
02649           if(adaptTempDTavenum == 100){
02650                 BigReal Frac;
02651                 adaptTempDTave /= adaptTempDTavenum;
02652                 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
02653                 Frac = adaptTempDTave/Frac;
02654                 //if average temperature jump is > 3% of temperature range,
02655                 //modify jump size to match 3%
02656                 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
02657                 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
02658                     Frac = adaptTempDtMax/Frac;
02659                     iout << "ADAPTEMP: Updating adaptTempDt to ";
02660                     adaptTempDt *= Frac;
02661                     iout << adaptTempDt << "\n" << endi;
02662                 }
02663                 adaptTempDTave = 0;
02664                 adaptTempDTavenum = 0;
02665           }
02666         }
02667       }
02668       else if (adaptTempAutoDt) {
02669           //update temperature step size counter
02670           // FOR "TRUE" ADAPTIVE TEMPERING
02671           BigReal adaptTempTdiff = fabs(dT-adaptTempT);
02672           if (adaptTempTdiff > 0) {
02673             adaptTempDTave += adaptTempTdiff; 
02674             adaptTempDTavenum++;
02675 //            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
02676           }
02677           if(adaptTempDTavenum == 100){
02678                 BigReal Frac;
02679                 adaptTempDTave /= adaptTempDTavenum;
02680                 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
02681                 Frac = adaptTempDTave/Frac;
02682                 //if average temperature jump is > 3% of temperature range,
02683                 //modify jump size to match 3%
02684                 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
02685                 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
02686                     Frac = adaptTempDtMax/Frac;
02687                     iout << "ADAPTEMP: Updating adaptTempDt to ";
02688                     adaptTempDt *= Frac;
02689                     iout << adaptTempDt << "\n" << endi;
02690                 }
02691                 adaptTempDTave = 0;
02692                 adaptTempDTavenum = 0;
02693 
02694           }
02695           
02696       }
02697       
02698       adaptTempT = dT; 
02699       broadcast->adaptTemperature.publish(step,adaptTempT);
02700     }
02701     adaptTempWriteRestart(step);
02702     if ( ! (step % adaptTempOutFreq) ) {
02703         iout << "ADAPTEMP: STEP " << step
02704              << " TEMP "   << adaptTempT
02705              << " ENERGY " << std::setprecision(10) << potentialEnergy   
02706              << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
02707              << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
02708         iout << "\n" << endi;
02709    }
02710    
02711 }
02712 
02713 
02714 void Controller::compareChecksums(int step, int forgiving) {
02715     Node *node = Node::Object();
02716     Molecule *molecule = node->molecule;
02717 
02718     // Some basic correctness checking
02719     BigReal checksum, checksum_b;
02720     char errmsg[256];
02721 
02722     checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
02723     if ( ((int)checksum) != molecule->numAtoms ) {
02724       sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
02725               (int)checksum, molecule->numAtoms);
02726       NAMD_bug(errmsg);
02727     }
02728 
02729     checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
02730     if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
02731       if ( ! computeChecksum ) {
02732         computesPartitioned = 0;
02733       } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
02734         computesPartitioned = 1;
02735       } else {
02736         NAMD_bug("Bad global compute count!\n");
02737       }
02738       computeChecksum = ((int)checksum);
02739     }
02740 
02741     checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
02742     if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
02743       sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
02744               (int)checksum, molecule->numCalcBonds);
02745       if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
02746         iout << iWARN << errmsg << endi;
02747       else NAMD_bug(errmsg);
02748     }
02749 
02750     checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
02751     if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
02752       sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
02753               (int)checksum, molecule->numCalcAngles);
02754       if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
02755         iout << iWARN << errmsg << endi;
02756       else NAMD_bug(errmsg);
02757     }
02758 
02759     checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
02760     if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
02761       sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
02762               (int)checksum, molecule->numCalcDihedrals);
02763       if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
02764         iout << iWARN << errmsg << endi;
02765       else NAMD_bug(errmsg);
02766     }
02767 
02768     checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
02769     if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
02770       sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
02771               (int)checksum, molecule->numCalcImpropers);
02772       if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
02773         iout << iWARN << errmsg << endi;
02774       else NAMD_bug(errmsg);
02775     }
02776 
02777     checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
02778     if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
02779       sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
02780               (int)checksum, molecule->numCalcTholes);
02781       if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
02782         iout << iWARN << errmsg << endi;
02783       else NAMD_bug(errmsg);
02784     }
02785 
02786     checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
02787     if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
02788       sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
02789               (int)checksum, molecule->numCalcAnisos);
02790       if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
02791         iout << iWARN << errmsg << endi;
02792       else NAMD_bug(errmsg);
02793     }
02794 
02795     checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
02796     if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
02797       sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
02798               (int)checksum, molecule->numCalcCrossterms);
02799       if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
02800         iout << iWARN << errmsg << endi;
02801       else NAMD_bug(errmsg);
02802     }
02803 
02804     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
02805     if ( ((int64)checksum) > molecule->numCalcExclusions &&
02806          ( ! simParams->mollyOn || step % slowFreq ) ) {
02807       if ( forgiving )
02808         iout << iWARN << "High global exclusion count ("
02809                       << ((int64)checksum) << " vs "
02810                       << molecule->numCalcExclusions << "), possible error!\n"
02811           << iWARN << "This warning is not unusual during minimization.\n"
02812           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02813       else {
02814         char errmsg[256];
02815         sprintf(errmsg, "High global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02816                 (long long)checksum, (long long)molecule->numCalcExclusions);
02817         NAMD_bug(errmsg);
02818       }
02819     }
02820 
02821     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
02822         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02823       if ( forgiving || ! simParams->fullElectFrequency ) {
02824         iout << iWARN << "Low global exclusion count!  ("
02825           << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
02826         if ( forgiving ) {
02827           iout << "\n"
02828             << iWARN << "This warning is not unusual during minimization.\n"
02829             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02830         } else {
02831           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02832         }
02833         iout << endi;
02834       } else {
02835         char errmsg[256];
02836         sprintf(errmsg, "Low global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
02837                 (long long)checksum, (long long)molecule->numCalcExclusions);
02838         NAMD_bug(errmsg);
02839       }
02840     }
02841 
02842 #ifdef NAMD_CUDA
02843     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
02844     if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
02845          ( ! simParams->mollyOn || step % slowFreq ) ) {
02846       if ( forgiving )
02847         iout << iWARN << "High global CUDA exclusion count ("
02848                       << ((int64)checksum) << " vs "
02849                       << molecule->numCalcFullExclusions << "), possible error!\n"
02850           << iWARN << "This warning is not unusual during minimization.\n"
02851           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02852       else {
02853         char errmsg[256];
02854         sprintf(errmsg, "High global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02855                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02856         NAMD_bug(errmsg);
02857       }
02858     }
02859 
02860     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
02861         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02862       if ( forgiving || ! simParams->fullElectFrequency ) {
02863         iout << iWARN << "Low global CUDA exclusion count!  ("
02864           << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
02865         if ( forgiving ) {
02866           iout << "\n"
02867             << iWARN << "This warning is not unusual during minimization.\n"
02868             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02869         } else {
02870           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02871         }
02872         iout << endi;
02873       } else {
02874         char errmsg[256];
02875         sprintf(errmsg, "Low global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
02876                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02877         NAMD_bug(errmsg);
02878       }
02879     }
02880 #endif
02881 
02882     checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
02883     if ( ((int)checksum) && ! marginViolations ) {
02884       iout << iERROR << "Margin is too small for " << ((int)checksum) <<
02885         " atoms during timestep " << step << ".\n" << iERROR <<
02886       "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
02887     }
02888     marginViolations += (int)checksum;
02889 
02890     checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
02891     if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
02892       iout << iINFO <<
02893         "Pairlistdist is too small for " << ((int)checksum) <<
02894         " computes during timestep " << step << ".\n" << endi;
02895     }
02896     if ( simParams->outputPairlists )  pairlistWarnings += (int)checksum;
02897 
02898     checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
02899     if ( checksum ) {
02900       if ( forgiving )
02901         iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
02902       else NAMD_bug("Stray PME grid charges detected!\n");
02903     }
02904 }
02905 
02906 void Controller::printTiming(int step) {
02907 
02908     if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
02909     {
02910       const double endWTime = CmiWallTimer() - firstWTime;
02911       const double endCTime = CmiTimer() - firstCTime;
02912 
02913       // fflush about once per minute
02914       if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
02915 
02916       const double elapsedW = 
02917         (endWTime - startWTime) / simParams->outputTiming;
02918       const double elapsedC = 
02919         (endCTime - startCTime) / simParams->outputTiming;
02920 
02921       const double remainingW = elapsedW * (simParams->N - step);
02922       const double remainingW_hours = remainingW / 3600;
02923 
02924       startWTime = endWTime;
02925       startCTime = endCTime;
02926 
02927 #ifdef NAMD_CUDA
02928       if ( simParams->outputEnergies < 60 &&
02929            step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
02930         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
02931       }
02932 #endif
02933       if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
02934         CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
02935                   ", %g hours remaining, %f MB of memory in use.\n",
02936                   step, endCTime, elapsedC, endWTime, elapsedW,
02937                   remainingW_hours, memusage_MB());
02938         if ( fflush_count ) { --fflush_count; fflush(stdout); }
02939       }
02940     }
02941 }
02942 
02943 void Controller::printMinimizeEnergies(int step) {
02944 
02945     rescaleaccelMD(step,1);
02946     receivePressure(step,1);
02947 
02948     Node *node = Node::Object();
02949     Molecule *molecule = node->molecule;
02950     compareChecksums(step,1);
02951 
02952     printEnergies(step,1);
02953 
02954     min_energy = totalEnergy;
02955     min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
02956     min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
02957     min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
02958     min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
02959 
02960 }
02961 
02962 void Controller::printDynamicsEnergies(int step) {
02963 
02964     Node *node = Node::Object();
02965     Molecule *molecule = node->molecule;
02966     SimParameters *simParameters = node->simParameters;
02967     Lattice &lattice = state->lattice;
02968 
02969     compareChecksums(step);
02970 
02971     printEnergies(step,0);
02972 }
02973 
02974 void Controller::printEnergies(int step, int minimize)
02975 {
02976     Node *node = Node::Object();
02977     Molecule *molecule = node->molecule;
02978     SimParameters *simParameters = node->simParameters;
02979     Lattice &lattice = state->lattice;
02980 
02981     // Drude model ANISO energy is added into BOND energy
02982     // and THOLE energy is added into ELECT energy
02983 
02984     BigReal bondEnergy;
02985     BigReal angleEnergy;
02986     BigReal dihedralEnergy;
02987     BigReal improperEnergy;
02988     BigReal crosstermEnergy;
02989     BigReal boundaryEnergy;
02990     BigReal miscEnergy;
02991     BigReal potentialEnergy;
02992     BigReal flatEnergy;
02993     BigReal smoothEnergy;
02994     BigReal work;
02995 
02996     Vector momentum;
02997     Vector angularMomentum;
02998     BigReal volume = lattice.volume();
02999 
03000     bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
03001     angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
03002     dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
03003     improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
03004     crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
03005     boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
03006     miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
03007 
03008     if ( minimize || ! ( step % nbondFreq ) )
03009     {
03010       electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
03011       ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
03012 
03013       // JLai
03014       groLJEnergy = reduction->item(REDUCTION_GRO_LJ_ENERGY);
03015       groGaussEnergy = reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
03016 
03017       goNativeEnergy = reduction->item(REDUCTION_GO_NATIVE_ENERGY);
03018       goNonnativeEnergy = reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
03019       goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
03020 
03021 //fepb
03022       bondedEnergyDiff_f = reduction->item(REDUCTION_BONDED_ENERGY_F);
03023       electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
03024       ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
03025 
03026       bondedEnergy_ti_1 = reduction->item(REDUCTION_BONDED_ENERGY_TI_1);
03027       electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
03028       ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
03029       bondedEnergy_ti_2 = reduction->item(REDUCTION_BONDED_ENERGY_TI_2);
03030       electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
03031       ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
03032 //fepe
03033     }
03034 
03035     if ( minimize || ! ( step % slowFreq ) )
03036     {
03037       electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
03038 //fepb
03039       electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
03040 
03041       electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
03042       electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
03043       electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
03044       electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
03045 //fepe
03046     }
03047 
03048     if (simParameters->LJcorrection && volume) {
03049 #ifdef MEM_OPT_VERSION
03050       NAMD_bug("LJcorrection not supported in memory optimized build.");
03051 #else
03052       // Apply tail correction to energy.
03053       BigReal alchLambda = simParameters->getCurrentLambda(step);
03054       ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
03055 
03056       if (simParameters->alchOn) {
03057         if (simParameters->alchFepOn) {
03058           BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
03059           ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
03060         } else if (simParameters->alchThermIntOn) {
03061           ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
03062           ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
03063         }
03064       }
03065 #endif
03066     }
03067 
03068 //fepb BKR - Compute alchemical work if using dynamic lambda.  This is here so
03069 //           that the cumulative work can be given during a callback.
03070     if (simParameters->alchLambdaFreq > 0) {
03071       if (step <= 
03072           simParameters->firstTimestep + simParameters->alchEquilSteps) {
03073         cumAlchWork = 0.0;
03074       }
03075       alchWork = computeAlchWork(step);
03076       cumAlchWork += alchWork;
03077     }
03078 //fepe
03079 
03080     momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
03081     momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
03082     momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
03083     angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
03084     angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
03085     angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
03086 
03087     // Ported by JLai
03088     potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
03089   + improperEnergy + electEnergy + electEnergySlow + ljEnergy
03090   + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy 
03091         + groLJEnergy + groGaussEnergy);
03092     // End of port
03093     totalEnergy = potentialEnergy + kineticEnergy;
03094     flatEnergy = (totalEnergy
03095         + (1.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
03096     if ( !(step%slowFreq) ) {
03097       // only adjust based on most accurate energies
03098       BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
03099       if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
03100       if ( step != simParams->firstTimestep ) {
03101         smooth2_avg *= 0.9375;
03102         smooth2_avg += 0.0625 * s;
03103       }
03104     }
03105     smoothEnergy = (flatEnergy + smooth2_avg
03106         - (4.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
03107 
03108     // Reset values for accumulated heat and work.
03109     if (step <= simParameters->firstTimestep &&
03110         (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
03111       heat = 0.0;
03112       totalEnergy0 = totalEnergy;
03113     }
03114     if ( simParameters->outputMomenta && ! minimize &&
03115          ! ( step % simParameters->outputMomenta ) )
03116     {
03117       iout << "MOMENTUM: " << step 
03118            << " P: " << momentum
03119            << " L: " << angularMomentum
03120            << "\n" << endi;
03121     }
03122 
03123     if ( simParameters->outputPressure ) {
03124       pressure_tavg += pressure;
03125       groupPressure_tavg += groupPressure;
03126       tavg_count += 1;
03127       if ( minimize || ! ( step % simParameters->outputPressure ) ) {
03128         iout << "PRESSURE: " << step << " "
03129            << PRESSUREFACTOR * pressure << "\n"
03130            << "GPRESSURE: " << step << " "
03131            << PRESSUREFACTOR * groupPressure << "\n";
03132         if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
03133            << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
03134            << "GPRESSAVG: " << step << " "
03135            << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
03136         iout << endi;
03137         pressure_tavg = 0;
03138         groupPressure_tavg = 0;
03139         tavg_count = 0;
03140       }
03141     }
03142 
03143     // pressure profile reductions
03144     if (pressureProfileSlabs) {
03145       const int freq = simParams->pressureProfileFreq;
03146       const int arraysize = 3*pressureProfileSlabs;
03147       
03148       BigReal *total = new BigReal[arraysize];
03149       memset(total, 0, arraysize*sizeof(BigReal));
03150       const int first = simParams->firstTimestep;
03151 
03152       if (ppbonded)    ppbonded->getData(first, step, lattice, total);
03153       if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
03154       if (ppint)       ppint->getData(first, step, lattice, total);
03155       for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
03156       pressureProfileCount++;
03157 
03158       if (!(step % freq)) {
03159         // convert NAMD internal virial to pressure in units of bar 
03160         BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
03161    
03162         iout << "PRESSUREPROFILE: " << step << " ";
03163         if (step == first) {
03164           // output pressure profile for this step
03165           for (int i=0; i<arraysize; i++) {
03166             iout << total[i] * scalefac << " ";
03167           }
03168         } else {
03169           // output pressure profile averaged over the last count steps.
03170           scalefac /= pressureProfileCount;
03171           for (int i=0; i<arraysize; i++) 
03172             iout << pressureProfileAverage[i]*scalefac << " ";
03173         }
03174         iout << "\n" << endi; 
03175        
03176         // Clear the average for the next block
03177         memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
03178         pressureProfileCount = 0;
03179       }
03180       delete [] total;
03181     }
03182   
03183    if ( step != simParams->firstTimestep || stepInFullRun == 0 ) {  // skip repeated first step
03184     if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
03185      int benchPhase = stepInFullRun / simParams->firstLdbStep;
03186      if ( benchPhase > 0 && benchPhase < 10 ) {
03187 #ifdef NAMD_CUDA
03188       if ( simParams->outputEnergies < 60 ) {
03189         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
03190       }
03191 #endif
03192       iout << iINFO;
03193       if ( benchPhase < 4 ) iout << "Initial time: ";
03194       else iout << "Benchmark time: ";
03195       iout << CkNumPes() << " CPUs ";
03196       {
03197         BigReal wallPerStep =
03198     (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
03199   BigReal ns = simParams->dt / 1000000.0;
03200   BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
03201   BigReal daysPerNano = wallPerStep * days / ns;
03202   iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
03203         iout << memusage_MB() << " MB memory\n" << endi;
03204       }
03205      }
03206      startBenchTime = CmiWallTimer();
03207     }
03208     ++stepInFullRun;
03209    }
03210 
03211     printTiming(step);
03212 
03213     Vector pairVDWForce, pairElectForce;
03214     if ( simParameters->pairInteractionOn ) {
03215       GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
03216       GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
03217     }
03218 
03219     // Compute cumulative nonequilibrium work (including shadow work).
03220     if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03221       work = totalEnergy - totalEnergy0 - heat;
03222     }
03223 
03224     // callback to Tcl with whatever we can
03225 #ifdef NAMD_TCL
03226 #define CALLBACKDATA(LABEL,VALUE) \
03227                 labels << (LABEL) << " "; values << (VALUE) << " ";
03228 #define CALLBACKLIST(LABEL,VALUE) \
03229                 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
03230     if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
03231       std::ostringstream labels, values;
03232 #if CMK_BLUEGENEL
03233       // the normal version below gives a compiler error
03234       values.precision(16);
03235 #else
03236       values << std::setprecision(16);
03237 #endif
03238       CALLBACKDATA("TS",step);
03239       CALLBACKDATA("BOND",bondEnergy);
03240       CALLBACKDATA("ANGLE",angleEnergy);
03241       CALLBACKDATA("DIHED",dihedralEnergy);
03242       CALLBACKDATA("CROSS",crosstermEnergy);
03243       CALLBACKDATA("IMPRP",improperEnergy);
03244       CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
03245       CALLBACKDATA("VDW",ljEnergy);
03246       CALLBACKDATA("BOUNDARY",boundaryEnergy);
03247       CALLBACKDATA("MISC",miscEnergy);
03248       CALLBACKDATA("POTENTIAL",potentialEnergy);
03249       CALLBACKDATA("KINETIC",kineticEnergy);
03250       CALLBACKDATA("TOTAL",totalEnergy);
03251       CALLBACKDATA("TEMP",temperature);
03252       CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
03253       CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
03254       CALLBACKDATA("VOLUME",lattice.volume());
03255       CALLBACKLIST("CELL_A",lattice.a());
03256       CALLBACKLIST("CELL_B",lattice.b());
03257       CALLBACKLIST("CELL_C",lattice.c());
03258       CALLBACKLIST("CELL_O",lattice.origin());
03259       labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
03260                 << lattice.b_p() << " " << lattice.c_p() << "} ";
03261       if (simParameters->drudeOn) {
03262         CALLBACKDATA("DRUDEBOND",drudeBondTemp);
03263       }
03264       if ( simParameters->pairInteractionOn ) {
03265         CALLBACKLIST("VDW_FORCE",pairVDWForce);
03266         CALLBACKLIST("ELECT_FORCE",pairElectForce);
03267       }
03268       if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03269         CALLBACKLIST("HEAT",heat);
03270         CALLBACKLIST("WORK",work);
03271       }
03272       if (simParameters->alchOn) {
03273         if (simParameters->alchThermIntOn) {
03274           CALLBACKLIST("BOND1", bondedEnergy_ti_1);
03275           CALLBACKLIST("ELEC1", electEnergy_ti_1 + electEnergySlow_ti_1 +
03276                                 electEnergyPME_ti_1);
03277           CALLBACKLIST("VDW1", ljEnergy_ti_1);
03278           CALLBACKLIST("BOND2", bondedEnergy_ti_2);
03279           CALLBACKLIST("ELEC2", electEnergy_ti_2 + electEnergySlow_ti_2 +
03280                                 electEnergyPME_ti_2);
03281           CALLBACKLIST("VDW2", ljEnergy_ti_2);
03282           if (simParameters->alchLambdaFreq > 0) {
03283             CALLBACKLIST("CUMALCHWORK", cumAlchWork);
03284           }
03285         } else if (simParameters->alchFepOn) {
03286           CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
03287                                 improperEnergy + bondedEnergyDiff_f);
03288           CALLBACKLIST("ELEC2", electEnergy_f + electEnergySlow_f);
03289           CALLBACKLIST("VDW2", ljEnergy_f);
03290         } 
03291       }
03292 
03293       labels << '\0';  values << '\0';  // insane but makes Linux work
03294       state->callback_labelstring = labels.str();
03295       state->callback_valuestring = values.str();
03296       // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
03297     }
03298 #undef CALLBACKDATA
03299 #endif
03300 
03301     drudeBondTempAvg += drudeBondTemp;
03302 
03303     temp_avg += temperature;
03304     pressure_avg += trace(pressure)/3.;
03305     groupPressure_avg += trace(groupPressure)/3.;
03306     avg_count += 1;
03307 
03308     if ( simParams->outputPairlists && pairlistWarnings &&
03309                                 ! (step % simParams->outputPairlists) ) {
03310       iout << iINFO << pairlistWarnings <<
03311         " pairlist warnings in past " << simParams->outputPairlists <<
03312         " steps.\n" << endi;
03313       pairlistWarnings = 0;
03314     }
03315 
03316     BigReal enthalpy;
03317     if (simParameters->multigratorOn && ((step % simParameters->outputEnergies) == 0)) {
03318       enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
03319     }
03320 
03321     // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
03322     if ( ! minimize &&  step % simParameters->outputEnergies ) return;
03323     // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
03324 
03325     if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
03326       IMDEnergies energies;
03327       energies.tstep = step;
03328       energies.T = temp_avg/avg_count;
03329       energies.Etot = totalEnergy;
03330       energies.Epot = potentialEnergy;
03331       energies.Evdw = ljEnergy;
03332       energies.Eelec = electEnergy + electEnergySlow;
03333       energies.Ebond = bondEnergy;
03334       energies.Eangle = angleEnergy;
03335       energies.Edihe = dihedralEnergy + crosstermEnergy;
03336       energies.Eimpr = improperEnergy;
03337       Node::Object()->imd->gather_energies(&energies);
03338     }
03339 
03340     if ( marginViolations ) {
03341       iout << iERROR << marginViolations <<
03342         " margin violations detected since previous energy output.\n" << endi;
03343     }
03344     marginViolations = 0;
03345 
03346     if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
03347     {
03348         iout << "ETITLE:      TS";
03349         iout << FORMAT("BOND");
03350         iout << FORMAT("ANGLE");
03351         iout << FORMAT("DIHED");
03352         if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
03353         iout << FORMAT("IMPRP");
03354         iout << "     ";
03355         iout << FORMAT("ELECT");
03356         iout << FORMAT("VDW");
03357         iout << FORMAT("BOUNDARY");
03358         iout << FORMAT("MISC");
03359         iout << FORMAT("KINETIC");
03360         iout << "     ";
03361         iout << FORMAT("TOTAL");
03362         iout << FORMAT("TEMP");
03363         iout << FORMAT("POTENTIAL");
03364         // iout << FORMAT("TOTAL2");
03365         iout << FORMAT("TOTAL3");
03366         iout << FORMAT("TEMPAVG");
03367         if ( volume != 0. ) {
03368           iout << "     ";
03369           iout << FORMAT("PRESSURE");
03370           iout << FORMAT("GPRESSURE");
03371           iout << FORMAT("VOLUME");
03372           iout << FORMAT("PRESSAVG");
03373           iout << FORMAT("GPRESSAVG");
03374         }
03375         if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03376           iout << "     ";
03377           iout << FORMAT("HEAT");
03378           iout << FORMAT("WORK");
03379         }
03380         if (simParameters->drudeOn) {
03381           iout << "     ";
03382           iout << FORMAT("DRUDEBOND");
03383           iout << FORMAT("DRBONDAVG");
03384         }
03385         // Ported by JLai
03386         if (simParameters->goGroPair) {
03387           iout << "     ";
03388           iout << FORMAT("GRO_PAIR_LJ");
03389           iout << FORMAT("GRO_PAIR_GAUSS");
03390         }
03391 
03392         if (simParameters->goForcesOn) {
03393           iout << "     ";
03394           iout << FORMAT("NATIVE");
03395           iout << FORMAT("NONNATIVE");
03396           //iout << FORMAT("REL_NATIVE");
03397           //iout << FORMAT("REL_NONNATIVE");
03398           iout << FORMAT("GOTOTAL");
03399           //iout << FORMAT("GOAVG");
03400         }
03401         // End of port -- JLai
03402 
03403         if (simParameters->alchOn) {
03404           if (simParameters->alchThermIntOn) {
03405             iout << "\nTITITLE:     TS";
03406             iout << FORMAT("BOND1");
03407             iout << FORMAT("ELECT1");
03408             iout << FORMAT("VDW1");
03409             iout << FORMAT("BOND2");
03410             iout << "     ";
03411             iout << FORMAT("ELECT2");
03412             iout << FORMAT("VDW2");
03413             if (simParameters->alchLambdaFreq > 0) {
03414               iout << FORMAT("LAMBDA");
03415               iout << FORMAT("ALCHWORK");
03416               iout << FORMAT("CUMALCHWORK");
03417             }
03418           } else if (simParameters->alchFepOn) {
03419             iout << "\nFEPTITLE:    TS";
03420             iout << FORMAT("BOND2");
03421             iout << FORMAT("ELECT2");
03422             iout << FORMAT("VDW2");
03423             if (simParameters->alchLambdaFreq > 0) {
03424               iout << FORMAT("LAMBDA");
03425             }
03426           }
03427         }
03428 
03429         iout << "\n\n" << endi;
03430         
03431         if (simParameters->qmForcesOn) {
03432             iout << "QMETITLE:      TS";
03433             iout << FORMAT("QMID");
03434             iout << FORMAT("ENERGY");
03435             if (simParameters->PMEOn) iout << FORMAT("PMECORRENERGY");
03436             iout << "\n\n" << endi;
03437         }
03438         
03439     }
03440 
03441     // N.B.  HP's aCC compiler merges FORMAT calls in the same expression.
03442     //       Need separate statements because data returned in static array.
03443     iout << ETITLE(step);
03444     iout << FORMAT(bondEnergy);
03445     iout << FORMAT(angleEnergy);
03446     if ( simParameters->mergeCrossterms ) {
03447       iout << FORMAT(dihedralEnergy+crosstermEnergy);
03448     } else {
03449       iout << FORMAT(dihedralEnergy);
03450       iout << FORMAT(crosstermEnergy);
03451     }
03452     iout << FORMAT(improperEnergy);
03453     iout << "     ";
03454     iout << FORMAT(electEnergy+electEnergySlow);
03455     iout << FORMAT(ljEnergy);
03456     iout << FORMAT(boundaryEnergy);
03457     iout << FORMAT(miscEnergy);
03458     iout << FORMAT(kineticEnergy);
03459     iout << "     ";
03460     iout << FORMAT(totalEnergy);
03461     iout << FORMAT(temperature);
03462     iout << FORMAT(potentialEnergy);
03463     // iout << FORMAT(flatEnergy);
03464     iout << FORMAT(smoothEnergy);
03465     iout << FORMAT(temp_avg/avg_count);
03466     if ( volume != 0. )
03467     {
03468         iout << "     ";
03469         iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
03470         iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
03471         iout << FORMAT(volume);
03472         iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
03473         iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
03474     }
03475     if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03476           iout << "     ";
03477           iout << FORMAT(heat);
03478           iout << FORMAT(work);
03479     }
03480     if (simParameters->drudeOn) {
03481         iout << "     ";
03482         iout << FORMAT(drudeBondTemp);
03483         iout << FORMAT(drudeBondTempAvg/avg_count);
03484     }
03485     // Ported by JLai
03486     if (simParameters->goGroPair) {
03487       iout << "     ";
03488       iout << FORMAT(groLJEnergy);
03489       iout << FORMAT(groGaussEnergy);
03490     }
03491 
03492     if (simParameters->goForcesOn) {
03493       iout << "     ";
03494       iout << FORMAT(goNativeEnergy);
03495       iout << FORMAT(goNonnativeEnergy);
03496       //iout << FORMAT(relgoNativeEnergy);
03497       //iout << FORMAT(relgoNonnativeEnergy);
03498       iout << FORMAT(goTotalEnergy);
03499       //iout << FORMAT("not implemented");
03500     } // End of port -- JLai
03501 
03502     if (simParameters->alchOn) { 
03503       if (simParameters->alchThermIntOn) {
03504         iout << "\n";
03505         iout << TITITLE(step);
03506         iout << FORMAT(bondedEnergy_ti_1);
03507         iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 + 
03508                        electEnergyPME_ti_1);
03509         iout << FORMAT(ljEnergy_ti_1);
03510         iout << FORMAT(bondedEnergy_ti_2);
03511         iout << "     ";
03512         iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 +
03513                        electEnergyPME_ti_2);
03514         iout << FORMAT(ljEnergy_ti_2);
03515         if (simParameters->alchLambdaFreq > 0) {
03516           iout << FORMAT(simParameters->getCurrentLambda(step));
03517           iout << FORMAT(alchWork);
03518           iout << FORMAT(cumAlchWork);
03519         }
03520       } else if (simParameters->alchFepOn) {
03521         iout << "\n";
03522         iout << FEPTITLE2(step);
03523         iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy 
03524                        + improperEnergy + bondedEnergyDiff_f);
03525         iout << FORMAT(electEnergy_f + electEnergySlow_f);
03526         iout << FORMAT(ljEnergy_f);
03527         if (simParameters->alchLambdaFreq > 0) {
03528           iout << FORMAT(simParameters->getCurrentLambda(step));
03529         }
03530       }
03531     }
03532 
03533     iout << "\n\n" << endi;
03534 
03535 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
03536      char webout[80];
03537      sprintf(webout,"%d %d %d %d",(int)totalEnergy,
03538              (int)(potentialEnergy),
03539              (int)kineticEnergy,(int)temperature);
03540      CApplicationDepositNode0Data(webout);
03541 #endif
03542 
03543     if (simParameters->pairInteractionOn) {
03544       iout << "PAIR INTERACTION:";
03545       iout << " STEP: " << step;
03546       iout << " VDW_FORCE: ";
03547       iout << FORMAT(pairVDWForce.x);
03548       iout << FORMAT(pairVDWForce.y);
03549       iout << FORMAT(pairVDWForce.z);
03550       iout << " ELECT_FORCE: ";
03551       iout << FORMAT(pairElectForce.x);
03552       iout << FORMAT(pairElectForce.y);
03553       iout << FORMAT(pairElectForce.z);
03554       iout << "\n" << endi;
03555     }
03556     drudeBondTempAvg = 0;
03557     temp_avg = 0;
03558     pressure_avg = 0;
03559     groupPressure_avg = 0;
03560     avg_count = 0;
03561 
03562     if ( fflush_count ) {
03563       --fflush_count;
03564       fflush(stdout);
03565     }
03566 }
03567 
03568 void Controller::writeExtendedSystemLabels(ofstream_namd &file) {
03569   Lattice &lattice = state->lattice;
03570   file << "#$LABELS step";
03571   if ( lattice.a_p() ) file << " a_x a_y a_z";
03572   if ( lattice.b_p() ) file << " b_x b_y b_z";
03573   if ( lattice.c_p() ) file << " c_x c_y c_z";
03574   file << " o_x o_y o_z";
03575   if ( simParams->langevinPistonOn ) {
03576     file << " s_x s_y s_z s_u s_v s_w";
03577   }
03578   file << std::endl;
03579 }
03580 
03581 void Controller::writeExtendedSystemData(int step, ofstream_namd &file) {
03582   Lattice &lattice = state->lattice;
03583   file.precision(12);
03584   file << step;
03585     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
03586     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
03587     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
03588     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
03589   if ( simParams->langevinPistonOn ) {
03590     Vector strainRate = diagonal(langevinPiston_strainRate);
03591     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
03592     file << " " << strainRate.x;
03593     file << " " << strainRate.y;
03594     file << " " << strainRate.z;
03595     file << " " << strainRate2.x;
03596     file << " " << strainRate2.y;
03597     file << " " << strainRate2.z;
03598   }
03599   file << std::endl;
03600 }
03601 
03602 void Controller::enqueueCollections(int timestep)
03603 {
03604   if ( Output::coordinateNeeded(timestep) ){
03605     collection->enqueuePositions(timestep,state->lattice);
03606   }
03607   if ( Output::velocityNeeded(timestep) )
03608     collection->enqueueVelocities(timestep);
03609   if ( Output::forceNeeded(timestep) )
03610     collection->enqueueForces(timestep);
03611 }
03612 
03613 //Modifications for alchemical fep
03614 void Controller::outputFepEnergy(int step) {
03615   if (simParams->alchOn && simParams->alchFepOn) {
03616     const int stepInRun = step - simParams->firstTimestep;
03617     const int alchEquilSteps = simParams->alchEquilSteps;
03618     const BigReal alchLambda = simParams->alchLambda;
03619     const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03620 
03621     if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
03622       FepNo = 0;
03623       exp_dE_ByRT = 0.0;
03624       net_dE = 0.0;
03625     }
03626     dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
03627          electEnergy - electEnergySlow - ljEnergy;
03628     BigReal RT = BOLTZMANN * simParams->alchTemp;
03629 
03630     if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
03631         (step / simParams->alchIDWSFreq) % 2 == 1 )) {
03632       // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
03633       FepNo++;
03634       exp_dE_ByRT += exp(-dE/RT);
03635       net_dE += dE;
03636     }
03637 
03638     if (simParams->alchOutFreq) {
03639       if (stepInRun == 0) {
03640         if (!fepFile.is_open()) {
03641           NAMD_backup_file(simParams->alchOutFile);
03642           fepFile.open(simParams->alchOutFile);
03643           iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
03644           if(alchEnsembleAvg){
03645             fepSum = 0.0;
03646             fepFile << "#            STEP                 Elec                            "
03647                     << "vdW                    dE           dE_avg         Temp             dG\n"
03648                     << "#                           l             l+dl      "
03649                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03650           }
03651           else{
03652             fepFile << "#            STEP                 Elec                            "
03653                     << "vdW                    dE         Temp\n"
03654                     << "#                           l             l+dl      "
03655                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03656           }
03657         }
03658         if(!step){
03659           fepFile << "#NEW FEP WINDOW: "
03660                   << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
03661                   << simParams->alchLambda2;
03662           if ( simParams->alchLambdaIDWS >= 0. ) {
03663             fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
03664           }
03665           fepFile << std::endl;
03666         }
03667       }
03668       if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03669         fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03670                 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03671                 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03672       }
03673       if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
03674         writeFepEnergyData(step, fepFile);
03675         fepFile.flush();
03676       }
03677       if (alchEnsembleAvg && (step == simParams->N)) {
03678         fepSum = fepSum + dG;
03679         fepFile << "#Free energy change for lambda window [ " << alchLambda
03680                 << " " << simParams->alchLambda2 << " ] is " << dG
03681                 << " ; net change until now is " << fepSum << std::endl;
03682         fepFile.flush();
03683       }
03684     }
03685   }
03686 }
03687 
03688 void Controller::outputTiEnergy(int step) {
03689  if (simParams->alchOn && simParams->alchThermIntOn) {
03690   const int stepInRun = step - simParams->firstTimestep;
03691   const int alchEquilSteps = simParams->alchEquilSteps;
03692   const int alchLambdaFreq = simParams->alchLambdaFreq;
03693 
03694   if (stepInRun == 0 || stepInRun == alchEquilSteps) {
03695     TiNo = 0;
03696     net_dEdl_bond_1 = 0;
03697     net_dEdl_bond_2 = 0;
03698     net_dEdl_elec_1 = 0;
03699     net_dEdl_elec_2 = 0;
03700     net_dEdl_lj_1 = 0;
03701     net_dEdl_lj_2 = 0;
03702   }
03703   TiNo++;
03704   net_dEdl_bond_1 += bondedEnergy_ti_1;
03705   net_dEdl_bond_2 += bondedEnergy_ti_2;
03706   net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
03707                       + electEnergyPME_ti_1);
03708   net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
03709                       + electEnergyPME_ti_2);
03710   net_dEdl_lj_1 += ljEnergy_ti_1;
03711   net_dEdl_lj_2 += ljEnergy_ti_2;
03712 
03713   // Don't accumulate block averages (you'll get a divide by zero!) or even 
03714   // open the TI output if alchOutFreq is zero.
03715   if (simParams->alchOutFreq) {
03716     if (stepInRun == 0 || stepInRun == alchEquilSteps 
03717         || (! ((step - 1) % simParams->alchOutFreq))) {
03718       // output of instantaneous dU/dl now replaced with running average
03719       // over last alchOutFreq steps (except for step 0)
03720       recent_TiNo = 0;
03721       recent_dEdl_bond_1 = 0;
03722       recent_dEdl_bond_2 = 0;
03723       recent_dEdl_elec_1 = 0;
03724       recent_dEdl_elec_2 = 0;
03725       recent_dEdl_lj_1 = 0;
03726       recent_dEdl_lj_2 = 0;
03727       recent_alchWork = 0;
03728     }
03729     recent_TiNo++;
03730     recent_dEdl_bond_1 += bondedEnergy_ti_1;
03731     recent_dEdl_bond_2 += bondedEnergy_ti_2;
03732     recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 
03733                            + electEnergyPME_ti_1); 
03734     recent_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 
03735                            + electEnergyPME_ti_2); 
03736     recent_dEdl_lj_1 += ljEnergy_ti_1;
03737     recent_dEdl_lj_2 += ljEnergy_ti_2;
03738     recent_alchWork += alchWork;
03739 
03740     if (stepInRun == 0) {
03741       if (!tiFile.is_open()) {
03742         NAMD_backup_file(simParams->alchOutFile);
03743         tiFile.open(simParams->alchOutFile);
03744         /* BKR - This has been rather drastically updated to better match 
03745            stdout. This was necessary for several reasons:
03746            1) PME global scaling is obsolete (now removed)
03747            2) scaling of bonded terms was added
03748            3) alchemical work is now accumulated when switching is active
03749          */
03750         iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
03751         tiFile << "#TITITLE:    TS";
03752         tiFile << FORMAT("BOND1");
03753         tiFile << FORMAT("AVGBOND1");
03754         tiFile << FORMAT("ELECT1");
03755         tiFile << FORMAT("AVGELECT1");
03756         tiFile << "     ";
03757         tiFile << FORMAT("VDW1");
03758         tiFile << FORMAT("AVGVDW1");
03759         tiFile << FORMAT("BOND2");
03760         tiFile << FORMAT("AVGBOND2");
03761         tiFile << FORMAT("ELECT2");
03762         tiFile << "     ";
03763         tiFile << FORMAT("AVGELECT2");
03764         tiFile << FORMAT("VDW2");
03765         tiFile << FORMAT("AVGVDW2");
03766         if (alchLambdaFreq > 0) {
03767           tiFile << FORMAT("ALCHWORK");
03768           tiFile << FORMAT("CUMALCHWORK");
03769         }
03770         tiFile << std::endl;
03771       }
03772 
03773       if (alchLambdaFreq > 0) {
03774         tiFile << "#ALCHEMICAL SWITCHING ACTIVE " 
03775                << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
03776                << "\n#LAMBDA SCHEDULE: " 
03777                << "dL: " << simParams->getLambdaDelta() 
03778                << " Freq: " << alchLambdaFreq;
03779       }
03780       else {
03781         const BigReal alchLambda = simParams->alchLambda;    
03782         const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
03783         const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
03784         const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
03785         const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
03786         const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
03787         const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
03788         tiFile << "#NEW TI WINDOW: "
03789                << "LAMBDA " << alchLambda 
03790                << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
03791                << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
03792                << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
03793                << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
03794       }
03795       tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
03796              << std::endl;
03797     }
03798 
03799     if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03800       tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03801              << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03802              << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03803     }
03804     if ((step%simParams->alchOutFreq) == 0) {
03805       writeTiEnergyData(step, tiFile);
03806       tiFile.flush();
03807     }
03808   }
03809  }
03810 }
03811 
03812 /* 
03813  Work is accumulated whenever alchLambda changes.  In the current scheme we
03814  always increment lambda _first_, then integrate in time.  Therefore the work 
03815  is wrt the "old" lambda before the increment.
03816 */
03817 BigReal Controller::computeAlchWork(const int step) {
03818   // alchemical scaling factors for groups 1/2 at the previous lambda
03819   const BigReal oldLambda = simParams->getCurrentLambda(step-1);
03820   const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
03821   const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
03822   const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
03823   const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
03824   const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
03825   const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
03826   // alchemical scaling factors for groups 1/2 at the new lambda
03827   const BigReal alchLambda = simParams->getCurrentLambda(step);
03828   const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
03829   const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
03830   const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
03831   const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
03832   const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
03833   const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda); 
03834 
03835   return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
03836           (elec_lambda_12 - elec_lambda_1)*
03837           (electEnergy_ti_1 + electEnergySlow_ti_1 +  electEnergyPME_ti_1) +
03838           (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
03839           (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
03840           (elec_lambda_22 - elec_lambda_2)*
03841           (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) + 
03842           (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
03843   );
03844 }
03845 
03846 void Controller::writeFepEnergyData(int step, ofstream_namd &file) {
03847   BigReal eeng = electEnergy + electEnergySlow;
03848   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
03849   BigReal RT = BOLTZMANN * simParams->alchTemp;
03850 
03851   const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03852   const int stepInRun = step - simParams->firstTimestep;
03853   const BigReal alchLambda = simParams->alchLambda;
03854 
03855   if(stepInRun){
03856     if ( simParams->alchLambdaIDWS >= 0. &&
03857         (step / simParams->alchIDWSFreq) % 2 == 0 ) {
03858       // IDWS is active and we are on a "backward" timestep
03859       fepFile << FEPTITLE_BACK(step);
03860     } else {
03861       // "forward" timestep
03862       fepFile << FEPTITLE(step);
03863     }
03864     fepFile << FORMAT(eeng);
03865     fepFile << FORMAT(eeng_f);
03866     fepFile << FORMAT(ljEnergy);
03867     fepFile << FORMAT(ljEnergy_f);
03868     fepFile << FORMAT(dE);
03869     if(alchEnsembleAvg){
03870       BigReal dE_avg = net_dE / FepNo;
03871       fepFile << FORMAT(dE_avg);
03872     }
03873     fepFile << FORMAT(temperature);
03874     if(alchEnsembleAvg){
03875       dG = -(RT * log(exp_dE_ByRT / FepNo));
03876       fepFile << FORMAT(dG);
03877     }
03878     fepFile << std::endl;
03879   }
03880 }
03881 
03882 void Controller::writeTiEnergyData(int step, ofstream_namd &file) {
03883   tiFile << TITITLE(step);
03884   tiFile << FORMAT(recent_dEdl_bond_1 / recent_TiNo);
03885   tiFile << FORMAT(net_dEdl_bond_1 / TiNo);
03886   tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
03887   tiFile << FORMAT(net_dEdl_elec_1/TiNo);
03888   tiFile << "     ";
03889   tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
03890   tiFile << FORMAT(net_dEdl_lj_1/TiNo);
03891   tiFile << FORMAT(recent_dEdl_bond_2 / recent_TiNo);
03892   tiFile << FORMAT(net_dEdl_bond_2 / TiNo);
03893   tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
03894   tiFile << "     ";
03895   tiFile << FORMAT(net_dEdl_elec_2/TiNo);
03896   tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
03897   tiFile << FORMAT(net_dEdl_lj_2/TiNo);
03898   if (simParams->alchLambdaFreq > 0) {
03899     tiFile << FORMAT(recent_alchWork / recent_TiNo);
03900     tiFile << FORMAT(cumAlchWork);
03901   }
03902   tiFile << std::endl;
03903 }
03904 //fepe
03905 
03906 void Controller::outputExtendedSystem(int step)
03907 {
03908 
03909   if ( step >= 0 ) {
03910 
03911     // Write out eXtended System Trajectory (XST) file
03912     if ( simParams->xstFrequency &&
03913          ((step % simParams->xstFrequency) == 0) )
03914     {
03915       if ( ! xstFile.is_open() )
03916       {
03917         iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
03918         NAMD_backup_file(simParams->xstFilename);
03919         xstFile.open(simParams->xstFilename);
03920         while (!xstFile) {
03921           if ( errno == EINTR ) {
03922             CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
03923             xstFile.clear();
03924             xstFile.open(simParams->xstFilename);
03925             continue;
03926           }
03927           char err_msg[257];
03928           sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
03929           NAMD_err(err_msg);
03930         }
03931         xstFile << "# NAMD extended system trajectory file" << std::endl;
03932         writeExtendedSystemLabels(xstFile);
03933       }
03934       writeExtendedSystemData(step,xstFile);
03935       xstFile.flush();
03936     }
03937 
03938     // Write out eXtended System Configuration (XSC) files
03939     //  Output a restart file
03940     if ( simParams->restartFrequency &&
03941          ((step % simParams->restartFrequency) == 0) &&
03942          (step != simParams->firstTimestep) )
03943     {
03944       iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
03945                 << step << "\n" << endi;
03946       char fname[140];
03947       const char *bsuffix = ".old";
03948       strcpy(fname, simParams->restartFilename);
03949       if ( simParams->restartSave ) {
03950         char timestepstr[20];
03951         sprintf(timestepstr,".%d",step);
03952         strcat(fname, timestepstr);
03953         bsuffix = ".BAK";
03954       }
03955       strcat(fname, ".xsc");
03956       NAMD_backup_file(fname,bsuffix);
03957       ofstream_namd xscFile(fname);
03958       while (!xscFile) {
03959         if ( errno == EINTR ) {
03960           CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
03961           xscFile.clear();
03962           xscFile.open(fname);
03963           continue;
03964         }
03965         char err_msg[257];
03966         sprintf(err_msg, "Error opening XSC restart file %s",fname);
03967         NAMD_err(err_msg);
03968       } 
03969       xscFile << "# NAMD extended system configuration restart file" << std::endl;
03970       writeExtendedSystemLabels(xscFile);
03971       writeExtendedSystemData(step,xscFile);
03972       if (!xscFile) {
03973         char err_msg[257];
03974         sprintf(err_msg, "Error writing XSC restart file %s",fname);
03975         NAMD_err(err_msg);
03976       } 
03977     }
03978 
03979   }
03980 
03981   //  Output final coordinates
03982   if (step == FILE_OUTPUT || step == END_OF_RUN)
03983   {
03984     int realstep = ( step == FILE_OUTPUT ?
03985         simParams->firstTimestep : simParams->N );
03986     iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
03987                 << realstep << "\n" << endi;
03988     static char fname[140];
03989     strcpy(fname, simParams->outputFilename);
03990     strcat(fname, ".xsc");
03991     NAMD_backup_file(fname);
03992     ofstream_namd xscFile(fname);
03993     while (!xscFile) {
03994       if ( errno == EINTR ) {
03995         CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
03996         xscFile.clear();
03997         xscFile.open(fname);
03998         continue;
03999       }
04000       char err_msg[257];
04001       sprintf(err_msg, "Error opening XSC output file %s",fname);
04002       NAMD_err(err_msg);
04003     } 
04004     xscFile << "# NAMD extended system configuration output file" << std::endl;
04005     writeExtendedSystemLabels(xscFile);
04006     writeExtendedSystemData(realstep,xscFile);
04007     if (!xscFile) {
04008       char err_msg[257];
04009       sprintf(err_msg, "Error writing XSC output file %s",fname);
04010       NAMD_err(err_msg);
04011     } 
04012   }
04013 
04014   //  Close trajectory file
04015   if (step == END_OF_RUN) {
04016     if ( xstFile.is_open() ) {
04017       xstFile.close();
04018       iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
04019     }
04020   }
04021 
04022 }
04023 
04024 
04025 void Controller::recvCheckpointReq(const char *key, int task, checkpoint &cp) {  // responding replica
04026   switch ( task ) {
04027     case SCRIPT_CHECKPOINT_STORE:
04028       if ( ! checkpoints.count(key) ) {
04029         checkpoints[key] = new checkpoint;
04030       }
04031       *checkpoints[key] = cp;
04032       break;
04033     case SCRIPT_CHECKPOINT_LOAD:
04034       if ( ! checkpoints.count(key) ) {
04035         NAMD_die("Unable to load checkpoint, requested key was never stored.");
04036       }
04037       cp = *checkpoints[key];
04038       break;
04039     case SCRIPT_CHECKPOINT_SWAP:
04040       if ( ! checkpoints.count(key) ) {
04041         NAMD_die("Unable to swap checkpoint, requested key was never stored.");
04042       }
04043       std::swap(cp,*checkpoints[key]);
04044       break;
04045     case SCRIPT_CHECKPOINT_FREE:
04046       if ( ! checkpoints.count(key) ) {
04047         NAMD_die("Unable to free checkpoint, requested key was never stored.");
04048       }
04049       delete checkpoints[key];
04050       checkpoints.erase(key);
04051       break;
04052   }
04053 }
04054 
04055 void Controller::recvCheckpointAck(checkpoint &cp) {  // initiating replica
04056   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
04057     state->lattice = cp.lattice;
04058     *(ControllerState*)this = cp.state;
04059   }
04060   CkpvAccess(_qd)->process();
04061 }
04062 
04063 
04064 void Controller::rebalanceLoad(int step)
04065 {
04066   if ( ! ldbSteps ) { 
04067     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
04068   }
04069   if ( ! --ldbSteps ) {
04070     startBenchTime -= CmiWallTimer();
04071         Node::Object()->outputPatchComputeMaps("before_ldb", step);
04072     LdbCoordinator::Object()->rebalance(this);  
04073         startBenchTime += CmiWallTimer();
04074     fflush_count = 3;
04075   }
04076 }
04077 
04078 void Controller::cycleBarrier(int doBarrier, int step) {
04079 #if USE_BARRIER
04080         if (doBarrier) {
04081           broadcast->cycleBarrier.publish(step,1);
04082           CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
04083                   CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04084         }
04085 #endif
04086 }
04087 
04088 void Controller::traceBarrier(int turnOnTrace, int step) {
04089         CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);       
04090         CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04091         nd.traceBarrier(turnOnTrace, step);
04092         CthSuspend();
04093 }
04094 
04095 void Controller::resumeAfterTraceBarrier(int step){
04096         broadcast->traceBarrier.publish(step,1);
04097         CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime); 
04098         awaken();
04099 }
04100 
04101 #ifdef MEASURE_NAMD_WITH_PAPI
04102 void Controller::papiMeasureBarrier(int turnOnMeasure, int step){
04103         CkPrintf("Cycle time at PAPI measurement sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);    
04104         CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04105         nd.papiMeasureBarrier(turnOnMeasure, step);
04106         CthSuspend();
04107 }
04108 
04109 void Controller::resumeAfterPapiMeasureBarrier(int step){
04110         broadcast->papiMeasureBarrier.publish(step,1);
04111         CkPrintf("Cycle time at PAPI measurement sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);      
04112         awaken();
04113 }
04114 #endif
04115 
04116 void Controller::terminate(void) {
04117   BackEnd::awaken();
04118   CthFree(thread);
04119   CthSuspend();
04120 }
04121