Controller.C

Go to the documentation of this file.

/*****************************************************************************
 * $Source: /home/cvs/namd/cvsroot/namd2/src/Controller.C,v $
 * $Author: jim $
 * $Date: 2017/03/29 21:29:03 $
 * $Revision: 1.1326 $
 *****************************************************************************/

#include "InfoStream.h"
#include "common.h"
#include "converse.h"
#include "memusage.h"
#include "Node.h"
#include "Molecule.h"
#include "SimParameters.h"
#include "Controller.h"
#include "ReductionMgr.h"
#include "CollectionMaster.h"
#include "Output.h"
#include "strlib.h"
#include "BroadcastObject.h"
#include "NamdState.h"
#include "ScriptTcl.h"
#include "Broadcasts.h"
#include "LdbCoordinator.h"
#include "Thread.h"
#include <math.h>
#include <signal.h>
#include "NamdOneTools.h"
#include "PatchMap.h"
#include "PatchMap.inl"
#include "Random.h"
#include "imd.h"
#include "IMDOutput.h"
#include "BackEnd.h"
#include <fstream>
#include <iomanip>
#include <errno.h>
#include "qd.h"
#include "PatchData.h"

#if defined(USE_GETTIMEOFDAY)
#include <sys/time.h>
// Replace CmiWallTimer with our own based on gettimeofday
static double namdWallTimer() {
  struct timeval t;
  gettimeofday(&t, NULL);
  double sec = double(t.tv_sec); // get absolute seconds
  double usec = double (t.tv_usec); // plus microseconds
  double totalSec = sec + 1e-6 * usec;
  return totalSec;
}
#else
static double namdWallTimer() {
  return CmiWallTimer();
}
#endif

#if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
extern "C" void CApplicationDepositNode0Data(char *);
#endif

//#define DEBUG_PRESSURE
#define MIN_DEBUG_LEVEL 3
//#define DEBUGM
#include "Debug.h"

#define XXXBIGREAL 1.0e32

//Modifications for alchemical fep
static char *FEPTITLE(int X)
{
  static char tmp_string[21];
  sprintf(tmp_string, "FepEnergy: %6d ",X);
  return tmp_string;
}

static char *FEPTITLE_BACK(int X)
{
  static char tmp_string[21];
  sprintf(tmp_string, "FepE_back: %6d ",X);
  return tmp_string;
}

static char *FEPTITLE2(int X)
{
  static char tmp_string[21];
  sprintf(tmp_string, "FEP:    %7d",X);
  return tmp_string;
}

static char *TITITLE(int X)
{
  static char tmp_string[21];
  sprintf(tmp_string, "TI:     %7d",X);
  return tmp_string;
}
//fepe

class PressureProfileReduction {
private:
  RequireReduction *reduction;
  const int nslabs;
  const int freq;
  int nelements;
  int count;
  char *name;

  BigReal *average;

public:
  PressureProfileReduction(int rtag, int numslabs, int numpartitions,
      const char *myname, int outputfreq)
  : nslabs(numslabs), freq(outputfreq) {
    name = strdup(myname);
    nelements = 3*nslabs * numpartitions;
    reduction = ReductionMgr::Object()->willRequire(rtag,nelements);

    average = new BigReal[nelements];
    count = 0;
  }
  ~PressureProfileReduction() {
    delete [] average;
    delete reduction;
    free(name);
  }
  // 
  void getData(int firsttimestep, int step, const Lattice &lattice, 
      BigReal *total) {
    reduction->require();

    int i;
    double inv_volume = 1.0 / lattice.volume();
    // accumulate the pressure profile computed for this step into the average.
    int arraysize = 3*nslabs;
    for (i=0; i<nelements; i++) {
      BigReal val = reduction->item(i) * inv_volume;
      average[i] += val;
      total[i % arraysize] += val;
    }
    count++;
    if (!(step % freq)) {
      // convert NAMD internal virial to pressure in units of bar 
      BigReal scalefac = PRESSUREFACTOR * nslabs;
 
      iout << "PPROFILE" << name << ": " << step << " ";
      if (step == firsttimestep) {
        // output pressure profile for this step
        for (i=0; i<nelements; i++) 
          iout << reduction->item(i)*scalefac*inv_volume << " ";
      } else {
        // output pressure profile averaged over the last Freq steps.
        scalefac /= count;
        for (i=0; i<nelements; i++) 
          iout << average[i]*scalefac << " ";
      }
      iout << "\n" << endi; 
      // Clear the average for the next block
      memset(average, 0, nelements*sizeof(BigReal));
      count = 0;
    }
  }
};

struct Controller::CthThreadWrapper {
  CthThread thread;
};

Controller::Controller(NamdState *s) :
        computeChecksum(0), marginViolations(0), pairlistWarnings(0),
        simParams(Node::Object()->simParameters),
        state(s),
        collection(CollectionMaster::Object()),
        startCTime(0),
        firstCTime(CmiTimer()),
        startWTime(0),
        firstWTime(namdWallTimer()),
        startBenchTime(0),
        stepInFullRun(0),
        ldbSteps(0),
        fflush_count(3)
{
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    if (simParams->CUDASOAintegrateMode) {
        CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
        PatchData* patchData = cpdata.ckLocalBranch();
        broadcast = new ControllerBroadcasts(0 /* ldObjPtr */, patchData->nodeBroadcast);
    } else 
#endif  // defined(NAMD_CUDA) || defined(NAMD_HIP)
    {
        broadcast = new ControllerBroadcasts;
    }
    min_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MINIMIZER,1);
    // for accelMD
    if (simParams->accelMDOn) {
       // REDUCTIONS_BASIC wil contain all potential energies and arrive first
       amd_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
       // contents of amd_reduction be submitted to REDUCTIONS_AMD
       submit_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_AMD);
       // REDUCTIONS_AMD will contain Sequencer contributions and arrive second
       reductionBasic = ReductionMgr::Object()->willRequire(REDUCTIONS_AMD);
    } else {
       amd_reduction = NULL;
       submit_reduction = NULL;
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
       if (simParams->CUDASOAintegrateMode) {
         reductionGpuResident = ReductionMgr::Object()->willRequire(REDUCTIONS_GPURESIDENT);
       }
#endif  // defined(NAMD_CUDA) || defined(NAMD_HIP)
       reductionBasic = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
    }
    // pressure profile reductions
    pressureProfileSlabs = 0;
    pressureProfileCount = 0;
    ppbonded = ppnonbonded = ppint = NULL;
    if (simParams->pressureProfileOn || simParams->pressureProfileEwaldOn) {
      int ntypes = simParams->pressureProfileAtomTypes;
      int nslabs = pressureProfileSlabs = simParams->pressureProfileSlabs;
      // number of partitions for pairwise interactions
      int npairs = (ntypes * (ntypes+1))/2;
      pressureProfileAverage = new BigReal[3*nslabs];
      memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
      int freq = simParams->pressureProfileFreq;
      if (simParams->pressureProfileOn) {
        ppbonded = new PressureProfileReduction(REDUCTIONS_PPROF_BONDED, 
            nslabs, npairs, "BONDED", freq);
        ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
            nslabs, npairs, "NONBONDED", freq);
        // internal partitions by atom type, but not in a pairwise fashion
        ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL, 
            nslabs, ntypes, "INTERNAL", freq);
      } else {
        // just doing Ewald, so only need nonbonded
        ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
            nslabs, npairs, "NONBONDED", freq);
      }
    }
    random = new Random(simParams->randomSeed);
    random->split(0,PatchMap::Object()->numPatches()+1);

    heat = totalEnergy0 = 0.0;
    rescaleVelocities_sumTemps = 0;  
    rescaleVelocities_numTemps = 0;
    stochRescale_count = 0;
    if (simParams->stochRescaleOn) {
      stochRescaleTimefactor = exp(-simParams->stochRescaleFreq *
          simParams->dt * 0.001 / simParams->stochRescalePeriod);
    }
    berendsenPressure_avg = 0; berendsenPressure_count = 0;
    // strainRate tensor is symmetric to avoid rotation
    langevinPiston_strainRate =
        Tensor::symmetric(simParams->strainRate,simParams->strainRate2);
    if ( ! simParams->useFlexibleCell ) {
      BigReal avg = trace(langevinPiston_strainRate) / 3.;
      langevinPiston_strainRate = Tensor::identity(avg);
    } else if ( simParams->useConstantRatio ) {
#define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
                 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
      AVGXY(langevinPiston_strainRate);
#undef AVGXY
    }
    langevinPiston_origStrainRate = langevinPiston_strainRate;

    // Monte Carlo pressure control
    if (simParams->monteCarloPressureOn) {
      monteCarloMaxVolume.x = simParams->monteCarloMaxVolume.x;
      monteCarloMaxVolume.y = simParams->monteCarloMaxVolume.y;
      monteCarloMaxVolume.z = simParams->monteCarloMaxVolume.z;
      mc_totalTry = mc_totalAccept = 0;
      for(int ax = 0; ax < MC_AXIS_TOTAL; ++ax) {
        mc_trial[ax] = 0;
        mc_accept[ax] = 0;
      }
    }

    if (simParams->multigratorOn) {
      multigratorXi = 0.0;
      int n = simParams->multigratorNoseHooverChainLength;
      BigReal tau = simParams->multigratorTemperatureRelaxationTime;
      Node *node = Node::Object();
      Molecule *molecule = node->molecule;
      BigReal Nf = molecule->num_deg_freedom();
      BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
      multigratorNu.resize(n);
      multigratorNuT.resize(n);
      multigratorZeta.resize(n);
      multigratorOmega.resize(n);
      for (int i=0;i < n;i++) {
        multigratorNu[i] = 0.0;
        multigratorZeta[i] = 0.0;
        if (i == 0) {
          multigratorOmega[i] = Nf*kT0*tau*tau;
        } else {
          multigratorOmega[i] = kT0*tau*tau;
        }
      }
      multigratorReduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
    } else {
      multigratorReduction = NULL;
    }
    origLattice = state->lattice;
    smooth2_avg = XXXBIGREAL;
    temp_avg = 0;
    pressure_avg = 0;
    groupPressure_avg = 0;
    avg_count = 0;
    pressure_tavg = 0;
    groupPressure_tavg = 0;
    tavg_count = 0;
    checkpoint_stored = 0;
    drudeBondTemp = 0;
    drudeBondTempAvg = 0;
    cumAlchWork = 0;

    bondedEnergy_ti_1 = 0;
    bondedEnergy_ti_2 = 0;
    ljEnergy_ti_1 = 0;
    ljEnergy_ti_2 = 0;
    electEnergy_ti_1 = 0;
    electEnergy_ti_2 = 0;
    electEnergySlow_ti_1 = 0;
    electEnergySlow_ti_2 = 0;
    TIderivative = 0;
}

Controller::~Controller(void)
{
    delete broadcast;
    if (reductionBasic) {
      delete reductionBasic;
    }
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    if (reductionGpuResident) {
      delete reductionGpuResident;
    }
#endif
    delete min_reduction;
    delete amd_reduction;
    delete submit_reduction;
    delete ppbonded;
    delete ppnonbonded;
    delete ppint;
    delete [] pressureProfileAverage;
    delete random;
    if (multigratorReduction) delete multigratorReduction;
}

void Controller::threadRun(Controller* arg)
{
    arg->algorithm();
}

void Controller::run(void)
{
    // create a Thread and invoke it
    DebugM(4, "Starting thread in controller on this=" << this << "\n");
    threadWrapper = new CthThreadWrapper;
    threadWrapper->thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
    CthSetStrategyDefault(threadWrapper->thread);
#if CMK_BLUEGENE_CHARM
    BgAttach(thread);
#endif
    awaken();
}

void Controller::awaken(void) {
  CthAwaken(threadWrapper->thread);
}

void Controller::suspend(void) {
  CthSuspend();
}

void Controller::algorithm(void)
{
  int scriptTask;
  int scriptSeq = 0;
  BackEnd::awaken();
  while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
#ifdef NODEGROUP_FORCE_REGISTER
    CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
    PatchData *patchData = cpdata.ckLocalBranch();
#endif
    switch ( scriptTask ) {
      case SCRIPT_OUTPUT:
        enqueueCollections(FILE_OUTPUT);
        outputExtendedSystem(FILE_OUTPUT);
        break;
      case SCRIPT_FORCEOUTPUT:
        enqueueCollections(FORCE_OUTPUT);
        break;
      case SCRIPT_MEASURE:
        enqueueCollections(EVAL_MEASURE);
        break;
      case SCRIPT_REINITVELS:
        iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
          << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
        break;
      case SCRIPT_RESCALEVELS:
        iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
          << " BY " << simParams->scriptArg1 << "\n" << endi;
        break;
      case SCRIPT_RESCALESOLUTECHARGES:
        // Parameter setting already reported in ScriptTcl
        // Nothing to do!
        break;
      case SCRIPT_CHECKPOINT:
        iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
          << "\n" << endi;
        checkpoint_stored = 1;
        checkpoint_lattice = state->lattice;
        checkpoint_state = *(ControllerState*)this;
        break;
      case SCRIPT_REVERT:
        iout << "REVERTING AT STEP " << simParams->firstTimestep
          << "\n" << endi;
        if ( ! checkpoint_stored )
          NAMD_die("Unable to revert, checkpoint was never called!");
        state->lattice = checkpoint_lattice;
        *(ControllerState*)this = checkpoint_state;
        break;
      case SCRIPT_CHECKPOINT_STORE:
        iout << "STORING CHECKPOINT AT STEP " << simParams->firstTimestep
          << " TO KEY " << simParams->scriptStringArg1;
        if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
        iout << "\n" << endi;
        if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
          if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
            checkpoints[simParams->scriptStringArg1] = new checkpoint;
          }
          checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
          cp.lattice = state->lattice;
          cp.state = *(ControllerState*)this;
        } else {
          Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                            scriptTask, state->lattice, *(ControllerState*)this);
        }
        break;
      case SCRIPT_CHECKPOINT_LOAD:
        iout << "LOADING CHECKPOINT AT STEP " << simParams->firstTimestep
          << " FROM KEY " << simParams->scriptStringArg1;
        if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
        iout << "\n" << endi;
        if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
          if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
            NAMD_die("Unable to load checkpoint, requested key was never stored.");
          }
          checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
          state->lattice = cp.lattice;
          *(ControllerState*)this = cp.state;
        } else {
          Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                            scriptTask, state->lattice, *(ControllerState*)this);
        }
        break;
      case SCRIPT_CHECKPOINT_SWAP:
        iout << "SWAPPING CHECKPOINT AT STEP " << simParams->firstTimestep
          << " FROM KEY " << simParams->scriptStringArg1;
        if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
        iout << "\n" << endi;
        if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
          if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
            NAMD_die("Unable to swap checkpoint, requested key was never stored.");
          }
          checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
          std::swap(state->lattice,cp.lattice);
          std::swap(*(ControllerState*)this,cp.state);
        } else {
          Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                            scriptTask, state->lattice, *(ControllerState*)this);
        }
        break;
      case SCRIPT_CHECKPOINT_FREE:
        iout << "FREEING CHECKPOINT AT STEP " << simParams->firstTimestep
          << " FROM KEY " << simParams->scriptStringArg1;
        if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
        iout << "\n" << endi;
        if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
          if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
            NAMD_die("Unable to free checkpoint, requested key was never stored.");
          }
          delete checkpoints[simParams->scriptStringArg1];
          checkpoints.erase(simParams->scriptStringArg1);
        } else {
          Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                            scriptTask, state->lattice, *(ControllerState*)this);
        }
        break;
      case SCRIPT_ATOMSENDRECV:
      case SCRIPT_ATOMSEND:
      case SCRIPT_ATOMRECV:
        break;
      case SCRIPT_MINIMIZE:
        minimize();
        break;
      case SCRIPT_RUN:
      case SCRIPT_CONTINUE:
        integrate(scriptTask);
        break;
      default:
        NAMD_bug("Unknown task in Controller::algorithm");
    }
    BackEnd::awaken();
  }
  enqueueCollections(END_OF_RUN);
  outputExtendedSystem(END_OF_RUN);
  terminate();
}


//
// XXX static and global variables are unsafe for shared memory builds.
// The use of global and static vars should be eliminated.
//
extern int eventEndOfTimeStep;

// Handle SIGINT so that restart files get written completely.
static int gotsigint = 0;
static void my_sigint_handler(int sig) {
  if (sig == SIGINT) gotsigint = 1;
}
extern "C" {
  typedef void (*namd_sighandler_t)(int);
}

void Controller::integrate(int scriptTask) {
    char traceNote[24];
  
    int step = simParams->firstTimestep;

    const int numberOfSteps = simParams->N;
    const int stepsPerCycle = simParams->stepsPerCycle;

    const int zeroMomentum = simParams->zeroMomentum;

    nbondFreq = simParams->nonbondedFrequency;
    const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
    if (dofull)
      slowFreq = simParams->fullElectFrequency;
    else
      slowFreq = simParams->nonbondedFrequency;
    if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;

   if(!simParams->CUDASOAintegrate){
     
     if ( scriptTask == SCRIPT_RUN ) {

       reassignVelocities(step);  // only for full-step velecities
       rescaleaccelMD(step);
       adaptTempUpdate(step); // Init adaptive tempering;

       receivePressure(step);
       if ( zeroMomentum && dofull && ! (step % slowFreq) )
                                                correctMomentum(step);
       printFepMessage(step);
       printTiMessage(step);
       printDynamicsEnergies(step);
       outputFepEnergy(step);
       outputTiEnergy(step);
       if(traceIsOn()){
           traceUserEvent(eventEndOfTimeStep);
           sprintf(traceNote, "s:%d", step);
           traceUserSuppliedNote(traceNote);
       }
       outputExtendedSystem(step);
       rebalanceLoad(step);

     }

    // Handling SIGINT doesn't seem to be working on Lemieux, and it
    // sometimes causes the net-xxx versions of NAMD to segfault on exit, 
    // so disable it for now.
    // namd_sighandler_t oldhandler = signal(SIGINT, 
    //  (namd_sighandler_t)my_sigint_handler);
    for ( ++step ; step <= numberOfSteps; ++step )
    {
        adaptTempUpdate(step);
        rescaleVelocities(step);
        tcoupleVelocities(step);
        stochRescaleVelocities(step);
        berendsenPressure(step);
        langevinPiston1(step);
        rescaleaccelMD(step);
        enqueueCollections(step);  // after lattice scaling!
        receivePressure(step);
        if ( zeroMomentum && dofull && ! (step % slowFreq) )
                                                correctMomentum(step);
            langevinPiston2(step);
        reassignVelocities(step);

        multigratorTemperature(step, 1);
        multigratorPressure(step, 1);
        multigratorPressure(step, 2);
        multigratorTemperature(step, 2);

        printDynamicsEnergies(step);
        outputFepEnergy(step);
        outputTiEnergy(step);
        if(traceIsOn()){
            traceUserEvent(eventEndOfTimeStep);
            sprintf(traceNote, "s:%d", step);
            traceUserSuppliedNote(traceNote);
        }
  // if (gotsigint) {
  //   iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
  //   NAMD_quit();
  // }
        outputExtendedSystem(step);
#if CYCLE_BARRIER
        cycleBarrier(!((step+1) % stepsPerCycle),step);
#elif  PME_BARRIER
        cycleBarrier(dofull && !(step%slowFreq),step);
#elif  STEP_BARRIER
        cycleBarrier(1, step);
#endif
                
        if(Node::Object()->specialTracing || simParams->statsOn){               
                 int bstep = simParams->traceStartStep;
                 int estep = bstep + simParams->numTraceSteps;          
                 if(step == bstep){
                         traceBarrier(1, step);
                 }
                 if(step == estep){                     
                         traceBarrier(0, step);
                 }
         }

#ifdef MEASURE_NAMD_WITH_PAPI
        if(simParams->papiMeasure) {
                int bstep = simParams->papiMeasureStartStep;
                int estep = bstep + simParams->numPapiMeasureSteps;
                if(step == bstep) {
                        papiMeasureBarrier(1, step);
                }
                if(step == estep) {
                        papiMeasureBarrier(0, step);
                }
        }
#endif
         
        rebalanceLoad(step);

#if  PME_BARRIER
        cycleBarrier(dofull && !((step+1)%slowFreq),step);   // step before PME
#endif
    }
  }else CthSuspend();
    // signal(SIGINT, oldhandler);
}

void Controller::resetMovingAverage() {
  const int size = simParams->movingAverageWindowSize;
  totalEnergyAverage.reset(size);
  temperatureAverage.reset(size);
  pressureAverage.reset(size);
  groupPressureAverage.reset(size);
}

#ifdef NODEGROUP_FORCE_REGISTER

void Controller::printStep(int step){
  char traceNote[24];
  const int numberOfSteps = simParams->N;
  const int stepsPerCycle = simParams->stepsPerCycle;
  const int zeroMomentum = simParams->zeroMomentum;
  nbondFreq = simParams->nonbondedFrequency;
  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
  if (dofull)
    slowFreq = simParams->fullElectFrequency;
  else
    slowFreq = simParams->nonbondedFrequency;
  if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;

if(step == simParams->firstTimestep /* && scriptTask == SCRIPT_RUN */){
    reassignVelocities(step);  // only for full-step velecities
    rescaleaccelMD(step);
    adaptTempUpdate(step); // Init adaptive tempering;

    receivePressure(step);
    if ( zeroMomentum && dofull && ! (step % slowFreq) )
           correctMomentum(step);
    printFepMessage(step);
    printTiMessage(step);
    printDynamicsEnergies(step);
    outputFepEnergy(step);
    outputTiEnergy(step);
    if(traceIsOn()){
        traceUserEvent(eventEndOfTimeStep);
        sprintf(traceNote, "s:%d", step);
        traceUserSuppliedNote(traceNote);
    }
    outputExtendedSystem(step);
  }else{
    adaptTempUpdate(step);
    rescaleVelocities(step);
    tcoupleVelocities(step);
    // XXX for single GPU version, call directly from Sequencer
    // stochRescaleVelocities(step);
    // XXX for single GPU version, call directly from Sequencer
    // berendsenPressure(step);
    // XXX for single GPU version, call directly from Sequencer
    // langevinPiston1(step);
    // XXX for single GPU version, call directly from Sequencer
    //  monteCarloPressure_prepare(step);
    //  monteCarloPressure_accept(step);
    rescaleaccelMD(step);
    enqueueCollections(step);  // after lattice scaling!
    receivePressure(step);
    if ( zeroMomentum && dofull && ! (step % slowFreq) )
        correctMomentum(step);
    langevinPiston2(step);
    reassignVelocities(step);

    multigratorTemperature(step, 1);
    multigratorPressure(step, 1);
    multigratorPressure(step, 2);
    multigratorTemperature(step, 2);

    printDynamicsEnergies(step);
    outputFepEnergy(step);
    outputTiEnergy(step);
    if(traceIsOn()){
      traceUserEvent(eventEndOfTimeStep);
      sprintf(traceNote, "s:%d", step);
      traceUserSuppliedNote(traceNote);
    }
    outputExtendedSystem(step);
  }
}

#endif


#define CALCULATE \
  printMinimizeEnergies(step); \
  outputExtendedSystem(step); \
  rebalanceLoad(step); \
  if ( step == numberOfSteps ) return; \
  else ++step;

#define MOVETO(X) \
  if ( step == numberOfSteps ) { \
    if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \
    if ( newDir || (mid.x-last.x) ) { \
      broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \
    } else { \
      broadcast->minimizeCoefficient.publish(minSeq++,0.); \
      broadcast->minimizeCoefficient.publish(minSeq++,0.); \
      min_reduction->require(); \
      broadcast->minimizeCoefficient.publish(minSeq++,0.); \
    } \
    enqueueCollections(step); \
    CALCULATE \
  } else if ( (X)-last.x ) { \
    broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \
    newDir = 0; \
    last.x = (X); \
    enqueueCollections(step); \
    CALCULATE \
    last.u = min_energy; \
    last.dudx = -1. * min_f_dot_v; \
    last.noGradient = min_huge_count; \
    if ( minVerbose ) { \
      iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \
      if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \
      iout << "\n" << endi; \
    } \
  }

struct minpoint {
  BigReal x, u, dudx, noGradient;
  minpoint() : x(0), u(0), dudx(0), noGradient(1) { ; }
};

void Controller::minimize() {
  // iout << "Controller::minimize() called.\n" << endi;

  const int minVerbose = simParams->minVerbose;
  const int numberOfSteps = simParams->N;
  int step = simParams->firstTimestep;
  slowFreq = nbondFreq = 1;
  BigReal linegoal = simParams->minLineGoal;  // 1.0e-3
  const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );

  CALCULATE

  int minSeq = 0;

  // just move downhill until initial bad contacts go away or 100 steps
  int old_min_huge_count = 2000000000;
  int steps_at_same_min_huge_count = 0;
  for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
    if ( min_huge_count >= old_min_huge_count ) {
      if ( ++steps_at_same_min_huge_count > 10 ) break;
    } else {
      old_min_huge_count = min_huge_count;
      steps_at_same_min_huge_count = 0;
    }
    iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
    broadcast->minimizeCoefficient.publish(minSeq++,1.);
    enqueueCollections(step);
    CALCULATE
  }
  if ( min_huge_count ) {
    iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
  }
  iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;

  int atStart = 2;
  int errorFactor = 10;
  BigReal old_f_dot_f = min_f_dot_f;
  broadcast->minimizeCoefficient.publish(minSeq++,0.);
  broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
  int newDir = 1;
  min_f_dot_v = min_f_dot_f;
  min_v_dot_v = min_f_dot_f;
  while ( 1 ) {
    // line minimization
    // bracket minimum on line
    minpoint lo,hi,mid,last;
    BigReal x = 0;
    lo.x = x;
    lo.u = min_energy;
    lo.dudx = -1. * min_f_dot_v;
    lo.noGradient = min_huge_count;
    mid = lo;
    last = mid;
    if ( minVerbose ) {
      iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
      if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
      iout << "\n" << endi;
    }
    BigReal tol = fabs( linegoal * min_f_dot_v );
    iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
            fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
    int start_with_huge = last.noGradient;
    min_reduction->require();
    BigReal maxstep = 0.1 / sqrt(min_reduction->item(0));
    x = maxstep; MOVETO(x);
    // bracket minimum on line
    while ( last.u < mid.u ) {
      if ( last.dudx < mid.dudx * (5.5 - x/maxstep)/5. ) {
        x = 6 * maxstep; break;
      }
      lo = mid; mid = last;
      x += maxstep;
      if ( x > 5.5 * maxstep ) break;
      MOVETO(x)
    }
    if ( x > 5.5 * maxstep ) {
      iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" << endi;
      broadcast->minimizeCoefficient.publish(minSeq++,0.);
      broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
      newDir = 1;
      old_f_dot_f = min_f_dot_f;
      min_f_dot_v = min_f_dot_f;
      min_v_dot_v = min_f_dot_f;
      continue;
    }
    hi = last;
#define PRINT_BRACKET \
    iout << "LINE MINIMIZER BRACKET: DX " \
         << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
        " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
        lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
    PRINT_BRACKET
    // converge on minimum on line
    int itcnt;
    for ( itcnt = 10; itcnt > 0; --itcnt ) {
      int progress = 1;
      // select new position
      if ( mid.noGradient ) {
       if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) {  // subdivide left side
        x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
        MOVETO(x)
        if ( last.u <= mid.u ) {
          hi = mid; mid = last;
        } else {
          lo = last;
        }
       } else {  // subdivide right side
        x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
        MOVETO(x)
        if ( last.u <= mid.u ) {
          lo = mid; mid = last;
        } else {
          hi = last;
        }
       }
      } else {
       if ( mid.dudx > 0. ) {  // subdivide left side
        BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
        BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
        x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
        BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
        if ( xdiv ) x /= xdiv; else x = last.x;
        if ( x > altxhi ) x = altxhi;
        if ( x < altxlo ) x = altxlo;
        if ( x-last.x == 0 ) break;
        MOVETO(x)
        if ( last.u <= mid.u ) {
          hi = mid; mid = last;
        } else {
          if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
          lo = last;
        }
       } else {  // subdivide right side
        BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
        BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
        x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
        BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
        if ( xdiv ) x /= xdiv; else x = last.x;
        if ( x < altxlo ) x = altxlo;
        if ( x > altxhi ) x = altxhi;
        if ( x-last.x == 0 ) break;
        MOVETO(x)
        if ( last.u <= mid.u ) {
          lo = mid; mid = last;
        } else {
          if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
          hi = last;
        }
       }
      }
      PRINT_BRACKET
      if ( ! progress ) {
        MOVETO(mid.x);
        break;
      }
      if ( fabs(last.dudx) < tol ) break;
      if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
      if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
    }
    // new direction
    broadcast->minimizeCoefficient.publish(minSeq++,0.);
    BigReal c = min_f_dot_f / old_f_dot_f;
    c = ( c > 1.5 ? 1.5 : c );
    if ( atStart ) { c = 0; --atStart; }
    if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
      c = 0;
      if ( errorFactor < 100 ) errorFactor += 10;
    }
    if ( c == 0 ) {
      iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
    }
    broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
    newDir = 1;
    old_f_dot_f = min_f_dot_f;
    min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
    min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
  }

}

#undef MOVETO
#undef CALCULATE

// NOTE: Only isotropic case implemented here!
void Controller::multigratorPressure(int step, int callNumber) {
  if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
    BigReal P0 = simParams->multigratorPressureTarget;
    BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
    BigReal NG = controlNumDegFreedom;
    BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
    BigReal tau = simParams->multigratorPressureRelaxationTime;
    BigReal s = 0.5*simParams->multigratorPressureFreq*simParams->dt/tau;
    {
      // Compute new scaling factors and send them to Sequencer
      BigReal V = state->lattice.volume();
      BigReal Pinst = trace(controlPressure)/3.0;
      BigReal PGsum = trace(momentumSqrSum);
      //
      multigratorXiT = multigratorXi + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
      BigReal scale = exp(s*NGfac*multigratorXiT);
      BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
      // fprintf(stderr, "%d | T %lf P %lf V %1.3lf\n", step, temperature, Pinst*PRESSUREFACTOR, V);
      Tensor scaleTensor = Tensor::identity(scale);
      Tensor volScaleTensor = Tensor::identity(scale);
      Tensor velScaleTensor = Tensor::identity(velScale);
      state->lattice.rescale(volScaleTensor);
      if (callNumber == 1) {
        broadcast->positionRescaleFactor.publish(step,scaleTensor);
        broadcast->velocityRescaleTensor.publish(step,velScaleTensor);
      } else {
        broadcast->positionRescaleFactor2.publish(step,scaleTensor);
        broadcast->velocityRescaleTensor2.publish(step,velScaleTensor);      
      }
    }

    {
      // Wait here for Sequencer to finish scaling and force computation
      RequireReduction* reduction = getCurrentReduction();
      reduction->require();
      Tensor virial_normal;
      Tensor virial_nbond;
      Tensor virial_slow;
      Tensor intVirial_normal;
      Tensor intVirial_nbond;
      Tensor intVirial_slow;
      Vector extForce_normal;
      Vector extForce_nbond;
      Vector extForce_slow;
      GET_TENSOR(momentumSqrSum, reduction, REDUCTION_MOMENTUM_SQUARED);
      GET_TENSOR(virial_normal, reduction, REDUCTION_VIRIAL_NORMAL);
      GET_TENSOR(virial_nbond, reduction, REDUCTION_VIRIAL_NBOND);
      GET_TENSOR(virial_slow, reduction, REDUCTION_VIRIAL_SLOW);
      GET_TENSOR(intVirial_normal, reduction, REDUCTION_INT_VIRIAL_NORMAL);
      GET_TENSOR(intVirial_nbond, reduction, REDUCTION_INT_VIRIAL_NBOND);
      GET_TENSOR(intVirial_slow, reduction, REDUCTION_INT_VIRIAL_SLOW);
      GET_VECTOR(extForce_normal, reduction, REDUCTION_EXT_FORCE_NORMAL);
      GET_VECTOR(extForce_nbond, reduction, REDUCTION_EXT_FORCE_NBOND);
      GET_VECTOR(extForce_slow, reduction, REDUCTION_EXT_FORCE_SLOW);
      calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
        intVirial_normal, intVirial_nbond, intVirial_slow,
        extForce_normal, extForce_nbond, extForce_slow);
      if (callNumber == 2) {
        // Update temperature for the Temperature Cycle that is coming next
        BigReal kineticEnergy = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
        temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
      }
    }

    {
      // Update pressure integrator
      BigReal V = state->lattice.volume();
      BigReal Pinst = trace(controlPressure)/3.0;
      BigReal PGsum = trace(momentumSqrSum);
      //
      multigratorXi = multigratorXiT + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
    }

  }
}

void Controller::multigratorTemperature(int step, int callNumber) {
  if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
    BigReal tau = simParams->multigratorTemperatureRelaxationTime;
    BigReal t = 0.5*simParams->multigratorTemperatureFreq*simParams->dt;
    BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
    BigReal Nf = numDegFreedom;
    int n = simParams->multigratorNoseHooverChainLength;
    BigReal T1, T2, v;
    {
      T1 = temperature;
      BigReal kTinst = BOLTZMANN * temperature;
      for (int i=n-1;i >= 0;i--) {
        if (i == 0) {
          BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
          multigratorNuT[0] = exp(-0.5*t*NuOmega)*multigratorNu[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
        } else if (i == n-1) {
          multigratorNuT[n-1] = multigratorNu[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
        } else {
          BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
          multigratorNuT[i] = exp(-0.5*t*NuOmega)*multigratorNu[i] + 
          0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
        }
      }
      BigReal velScale = exp(-t*multigratorNuT[0]/multigratorOmega[0]);
      v = velScale;
      if (callNumber == 1)
        broadcast->velocityRescaleFactor.publish(step,velScale);
      else
        broadcast->velocityRescaleFactor2.publish(step,velScale);
    }

    {
      // Wait here for Sequencer to finish scaling and re-calculating kinetic energy
      multigratorReduction->require();
      BigReal kineticEnergy = multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY);
      temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
      T2 = temperature;
      if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
        // If this is pressure cycle, receive new momentum product
        GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
      }
    }

    // fprintf(stderr, "%d | T %lf scale %lf T' %lf\n", step, T1, v, T2);

    {
      BigReal kTinst = BOLTZMANN * temperature;
      for (int i=0;i < n;i++) {
        if (i == 0) {
          BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
          multigratorNu[0] = exp(-0.5*t*NuOmega)*multigratorNuT[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
        } else if (i == n-1) {
          multigratorNu[n-1] = multigratorNuT[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
        } else {
          BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
          multigratorNu[i] = exp(-0.5*t*NuOmega)*multigratorNuT[i] + 
          0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
        }
        multigratorZeta[i] += t*multigratorNuT[i]/multigratorOmega[i];
      }
    }

  }
}

// Calculates Enthalpy for multigrator
BigReal Controller::multigatorCalcEnthalpy(BigReal potentialEnergy, int step, int minimize) {
  Node *node = Node::Object();
  Molecule *molecule = node->molecule;
  SimParameters *simParameters = node->simParameters;

  BigReal V = state->lattice.volume();
  BigReal P0 = simParams->multigratorPressureTarget;
  BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
  BigReal NG = controlNumDegFreedom;
  BigReal Nf = numDegFreedom;
  BigReal tauV = simParams->multigratorPressureRelaxationTime;
  BigReal sumZeta = 0.0;
  for (int i=1;i < simParams->multigratorNoseHooverChainLength;i++) {
    sumZeta += multigratorZeta[i];
  }
  BigReal nuOmegaSum = 0.0;
  for (int i=0;i < simParams->multigratorNoseHooverChainLength;i++) {
    nuOmegaSum += multigratorNu[i]*multigratorNu[i]/(2.0*multigratorOmega[i]);
  }
  BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
  BigReal eta = sqrt(kT0*W)*multigratorXi;

  BigReal enthalpy = kineticEnergy + potentialEnergy + eta*eta/(2.0*W) + P0*V + nuOmegaSum + kT0*(Nf*multigratorZeta[0] + sumZeta);

//  if (!(step % 100))
    // fprintf(stderr, "enthalpy %lf %lf %lf %lf %lf %lf %lf\n", enthalpy,
    //   kineticEnergy, potentialEnergy, eta*eta/(2.0*W), P0*V, nuOmegaSum, kT0*(Nf*multigratorZeta[0] + sumZeta));

  return enthalpy;
}

void Controller::monteCarloPressure_accept(int step)
{
  if (simParams->monteCarloPressureOn) {
    if ( !(step % simParams->monteCarloPressureFreq) ) {    
      Node *node = Node::Object();
      Molecule *molecule = node->molecule;
      // need to consider the fix atom case!
      int numAtom = molecule->numMolecules;
      BigReal volumeOld = mc_oldLattice.volume();
      BigReal volumeNew = state->lattice.volume();
      BigReal areaNew = (volumeNew / state->lattice.c().z);
      BigReal areaOld = (volumeOld / mc_oldLattice.c().z);
      BigReal kbT = BOLTZMANN * simParams->monteCarloTemp;
      BigReal pressureTarget = simParams->monteCarloPressureTarget;
      BigReal surfTensionTarget = simParams->surfaceTensionTarget;

      // When using the GPU resident code path in single process mode, the reduction
      // data needs to persist so it can be consumed by the printStep function
      const bool clear_reduction_data = !simParams->GPUresidentSingleProcessMode;
      RequireReduction* reduction = getCurrentReduction();
      reduction->require(clear_reduction_data);
      
      BigReal totalEnergyOld = mc_totalEnergyOld;
      BigReal totalEnergyNew = getTotalPotentialEnergy(step);
      BigReal weight = (totalEnergyNew - totalEnergyOld);
      weight += (pressureTarget * (volumeNew - volumeOld));
      weight -= (numAtom * kbT * log(volumeNew / volumeOld));
      weight -= (surfTensionTarget * (areaNew - areaOld));
      BigReal totalWeight = exp(-weight / kbT);
      int accepted = (random->uniform() < totalWeight);
      //accepted = 0;
      #if 0
      printf("MC-data (accept): step: %d, Pe: %d, numAtom: %d\n", step, CkMyPe(), numAtom);
      printf("                  oldE: %f, newE: %f, deltaE: %f\n", totalEnergyOld, totalEnergyNew, totalEnergyNew - totalEnergyOld);
      printf("                  oldV: %f, newV: %f, deltaV: %f\n", volumeOld, volumeNew, volumeNew - volumeOld);
      printf("                  surfTen: %f, deltaA: %f \n", surfTensionTarget, areaNew - areaOld);
      printf("                  weight: %f, accepted: %d\n\n", totalWeight, accepted);
      #endif
      if(accepted) {
        // what should I do if its accepted?
        ++mc_totalAccept;
        ++mc_accept[mc_picked_axis];
      } else {
        // what should I do if its rejected?
        state->lattice = mc_oldLattice;
      }
      // update the maximum scale no matter if it's accepted or not
      {
        if (mc_trial[mc_picked_axis] == simParams->monteCarloAdjustmentFreq) {
          BigReal factor = 1.0;
          BigReal accRate = (BigReal) mc_accept[mc_picked_axis] / (BigReal) mc_trial[mc_picked_axis];
          // safe check to make sure there is any accepted move
          if (accRate > 0.0) {
            factor = accRate / simParams->monteCarloAcceptanceRate;
          } else {
            factor = 0.5; // if no move was accepted, we scale by half
          }
          // reset the trial for picked axis
          mc_trial[mc_picked_axis] = mc_accept[mc_picked_axis] = 0;
          BigReal maxVolume = state->lattice.volume() * 0.3;
          switch (mc_picked_axis) {
            case MC_X:
              monteCarloMaxVolume.x *= factor;
              if (monteCarloMaxVolume.x < 0.001) {
                iout << iWARN << "Small volume change in MC barostat!\n" << endi;
                iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
                monteCarloMaxVolume.x = 0.001;
              } else if (monteCarloMaxVolume.x > maxVolume ) {
                iout << iWARN << "Large volume change in MC barostat!\n" << endi;
                iout << iWARN << "Maximum volume change is set to " <<  maxVolume << " A^3.\n" << endi;
                monteCarloMaxVolume.x = maxVolume;
              }
              break;

            case MC_Y:
              monteCarloMaxVolume.y *= factor;
              if (monteCarloMaxVolume.y < 0.001) {
                iout << iWARN << "Small volume change in MC barostat!\n" << endi;
                iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
                monteCarloMaxVolume.y = 0.001;
              } else if (monteCarloMaxVolume.y > maxVolume ) {
                iout << iWARN << "Large volume change in MC barostat!\n" << endi;
                iout << iWARN << "Maximum volume change is set to " <<  maxVolume << " A^3.\n" << endi;
                monteCarloMaxVolume.y = maxVolume;
              }
              break;

            case MC_Z:
              monteCarloMaxVolume.z *= factor;
              if (monteCarloMaxVolume.z < 0.001) {
                iout << iWARN << "Small volume change in MC barostat!\n" << endi;
                iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
                monteCarloMaxVolume.z = 0.001;
              } else if (monteCarloMaxVolume.z > maxVolume ) {
                iout << iWARN << "Large volume change in MC barostat!\n" << endi;
                iout << iWARN << "Maximum volume change is set to " <<  maxVolume << " A^3.\n" << endi;
                monteCarloMaxVolume.z = maxVolume;
              }
              break;
          }
        }
      }

      if ( !(step % simParams->outputEnergies) ) {
        BigReal acceptPerc = 100.0 * ((BigReal) mc_totalAccept /
                            (BigReal) mc_totalTry);
        iout << "MC_BAROSTAT STEP: " << step << " TRIALS: " << mc_totalTry << 
          " ACCEPTED: " << mc_totalAccept << " %ACCEPTANCE: " << acceptPerc;
        if (simParams->useFlexibleCell) {
          if (simParams->useConstantArea) {
            iout << " MAX_VOLUME_Z: " << monteCarloMaxVolume.z << "\n\n";
          } else if (simParams->useConstantRatio) {
            iout << " MAX_VOLUME_XY: " << monteCarloMaxVolume.x << 
              " MAX_VOLUME_Z: " << monteCarloMaxVolume.z << "\n\n";
          } else {
            iout << " MAX_VOLUME_X: " << monteCarloMaxVolume.x
                 << " MAX_VOLUME_Y: " << monteCarloMaxVolume.y
                 << " MAX_VOLUME_Z: " << monteCarloMaxVolume.z
                 << "\n\n";
          }
        } else {
          iout << " MAX_VOLUME_XYZ: " << monteCarloMaxVolume.x << "\n\n";
        }
      }

      broadcast->monteCarloBarostatAcceptance.publish(step, accepted);
    }
  }
}

void Controller::monteCarloPressure_prepare(int step)
{
  if (simParams->monteCarloPressureOn) {
    if ( !(step % simParams->monteCarloPressureFreq) ) {   
      // Wait here for Sequencer to finish energy/force computation
      // When using the GPU resident code path in single process mode, the reduction
      // data needs to persist so it can be consumed by the MCPressure accept function
      const bool clear_reduction_data = !simParams->GPUresidentSingleProcessMode;
      RequireReduction* reduction = getCurrentReduction();
      reduction->require(clear_reduction_data);
      mc_oldLattice = state->lattice;
      mc_totalEnergyOld = getTotalPotentialEnergy(step);
      Tensor factor = Tensor::identity(1.0);
      // pick a random deltaV and calc scaling factor
      {
        BigReal oldVolume = mc_oldLattice.volume();
        bool flexible = simParams->useFlexibleCell;
        bool constArea = simParams->useConstantArea;
        bool constRatio = simParams->useConstantRatio;
        bool constZ, constY, constX;
        if (simParams->fixCellDims) {
          constX = simParams->fixCellDimX;
          constY = simParams->fixCellDimY;
          constZ = simParams->fixCellDimZ;
        } else {
          constZ = constY = constX = false;
        }

        if(flexible) {
          //Anisotropic fluctuation
          if(constArea) {
            // always pick z
            mc_picked_axis = MC_Z;
            BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
            factor.zz = (oldVolume + dV) / oldVolume;
          } else if (constRatio) {
            while(true) {
              // decide to scale x-y or scale z
              BigReal probAxis = 2.0 * random->uniform();
              if (probAxis < 1.0) {
                //picking x or y. We use MC_X for constant ratio
                mc_picked_axis = MC_X;
                BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
                BigReal scale = sqrt((oldVolume + dV)/oldVolume);
                factor.xx = factor.yy = scale;
                break;
              } else if (probAxis < 2.0 && !constZ) {
                //picking z
                mc_picked_axis = MC_Z;
                BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
                BigReal scale = (oldVolume + dV)/oldVolume;
                factor.zz = scale;
                break;
              }
            }
          } else {
            while(true) {
              // decide to scale x, y or z
              BigReal probAxis = 3.0 * random->uniform();
              if (probAxis < 1.0 && !constX) {
                // picking x
                mc_picked_axis = MC_X;
                BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
                factor.xx = (oldVolume + dV) / oldVolume;
                break;
              } else if (probAxis < 2.0 && !constY) {
                // picking y
                mc_picked_axis = MC_Y;
                BigReal dV = monteCarloMaxVolume.y * (2.0 * random->uniform() - 1.0);
                factor.yy = (oldVolume + dV) / oldVolume;
                break;
              } else if (probAxis < 3.0 && !constZ) {
                // picking z
                mc_picked_axis = MC_Z;
                BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
                factor.zz = (oldVolume + dV) / oldVolume;
                break;
              }
            }
          }
        } else {
          //Isotropic fluctuation.  We use MC_X for isotropic fluctuation
          mc_picked_axis = MC_X;
          BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
          BigReal scale = pow((oldVolume + dV)/oldVolume, 1.0/3.0);
          factor.xx = factor.yy = factor.zz = scale;
        } 
        // increament the MC trial
        ++mc_totalTry;
        ++mc_trial[mc_picked_axis];
      }

      broadcast->positionRescaleFactor.publish(step,factor);
      state->lattice.rescale(factor);
    }
  }
}


#define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
  if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
  if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
  
void Controller::berendsenPressure(int step)
{
  if ( simParams->berendsenPressureOn ) {
   berendsenPressure_count += 1;
   berendsenPressure_avg += controlPressure;
   const int freq = simParams->berendsenPressureFreq;
   if ( ! (berendsenPressure_count % freq) ) {
    Tensor factor = berendsenPressure_avg / berendsenPressure_count;
    berendsenPressure_avg = 0;
    berendsenPressure_count = 0;
    // We only use on-diagonal terms (for now)
    factor = Tensor::diagonal(diagonal(factor));
    factor -= Tensor::identity(simParams->berendsenPressureTarget);
    factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
       simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
    factor += Tensor::identity(1.0);
    int limited = 0;
    LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
    LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
    LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
    if ( limited ) {
      iout << iERROR << "Step " << step <<
        " cell rescaling factor limited.\n" << endi;
    }
    broadcast->positionRescaleFactor.publish(step,factor);
    state->lattice.rescale(factor);
   }
  } else {
    berendsenPressure_avg = 0;
    berendsenPressure_count = 0;
  }
}

#undef LIMIT_SCALING

void Controller::langevinPiston1(int step)
{
  if ( simParams->langevinPistonOn )
  {
    Tensor &strainRate = langevinPiston_strainRate;
    int cellDims = simParams->useFlexibleCell ? 1 : 3;
    BigReal dt = simParams->dt;
    BigReal dt_long = slowFreq * dt;
    BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
    BigReal tau = simParams->langevinPistonPeriod;
    BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;

#ifdef DEBUG_PRESSURE
    iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
#endif

    if ( ! ( (step-1) % slowFreq ) )
    {
      BigReal gamma = 1 / simParams->langevinPistonDecay;
      BigReal f1 = exp( -0.5 * dt_long * gamma );
      BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
      strainRate *= f1;
      if ( simParams->useFlexibleCell ) {
        // We only use on-diagonal terms (for now)
        if ( simParams->useConstantRatio ) {
          BigReal r = f2 * random->gaussian();
          strainRate.xx += r;
          strainRate.yy += r;
          strainRate.zz += f2 * random->gaussian();
        } else {
          strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
        }
      } else {
        strainRate += f2 * Tensor::identity(random->gaussian());
      }

#ifdef DEBUG_PRESSURE
      iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
#endif
    }

    // Apply surface tension.  If surfaceTensionTarget is zero, we get
    // the default (isotropic pressure) case.
    
    Tensor ptarget;
    ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
    if ( simParams->surfaceTensionTarget != 0. ) {
      ptarget.xx = ptarget.yy = simParams->langevinPistonTarget - 
        simParams->surfaceTensionTarget / state->lattice.c().z;
    }

    strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
      ( controlPressure - ptarget );

    if (simParams->fixCellDims) {
      if (simParams->fixCellDimX) strainRate.xx = 0;
      if (simParams->fixCellDimY) strainRate.yy = 0;
      if (simParams->fixCellDimZ) strainRate.zz = 0;
    }


#ifdef DEBUG_PRESSURE
    iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
#endif

    if ( simParams->langevinPistonBarrier ) {
    if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
    {
      // We only use on-diagonal terms (for now)
      Tensor factor;
      if ( !simParams->useConstantArea ) {
        factor.xx = exp( dt_long * strainRate.xx );
        factor.yy = exp( dt_long * strainRate.yy );
      } else {
        factor.xx = factor.yy = 1;
      }
      factor.zz = exp( dt_long * strainRate.zz );
      broadcast->positionRescaleFactor.publish(step,factor);
      state->lattice.rescale(factor);
#ifdef DEBUG_PRESSURE
      iout << iINFO << "rescaling by: " << factor << "\n";
#endif
    }
    } else { // langevinPistonBarrier
    if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
    {
      if ( ! positionRescaleFactor.xx ) {  // first pass without MTS
      // We only use on-diagonal terms (for now)
      Tensor factor;
      if ( !simParams->useConstantArea ) {
        factor.xx = exp( dt_long * strainRate.xx );
        factor.yy = exp( dt_long * strainRate.yy );
      } else {
        factor.xx = factor.yy = 1;
      }
      factor.zz = exp( dt_long * strainRate.zz );
      broadcast->positionRescaleFactor.publish(step,factor);
      positionRescaleFactor = factor;
      strainRate_old = strainRate;
      }
      state->lattice.rescale(positionRescaleFactor);
#ifdef DEBUG_PRESSURE
      iout << iINFO << "rescaling by: " << positionRescaleFactor << "\n";
#endif
    }
    if ( ! ( (step-slowFreq/2) % slowFreq ) )
    {
      // We only use on-diagonal terms (for now)
      Tensor factor;
      if ( !simParams->useConstantArea ) {
        factor.xx = exp( (dt+dt_long) * strainRate.xx - dt * strainRate_old.xx );
        factor.yy = exp( (dt+dt_long) * strainRate.yy - dt * strainRate_old.yy );
      } else {
        factor.xx = factor.yy = 1;
      }
      factor.zz = exp( (dt+dt_long) * strainRate.zz - dt * strainRate_old.zz );
      broadcast->positionRescaleFactor.publish(step+1,factor);      
      positionRescaleFactor = factor;
      strainRate_old = strainRate;
    }
    }

    // corrections to integrator
    if ( ! ( step % nbondFreq ) )
    {
#ifdef DEBUG_PRESSURE
      iout << iINFO << "correcting strain rate for nbond, ";
#endif
      strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                ( (nbondFreq - 1) * controlPressure_nbond );
#ifdef DEBUG_PRESSURE
      iout << "strain rate: " << strainRate << "\n";
#endif
    }
    if ( ! ( step % slowFreq ) )
    {
#ifdef DEBUG_PRESSURE
      iout << iINFO << "correcting strain rate for slow, ";
#endif
      strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                ( (slowFreq - 1) * controlPressure_slow );
#ifdef DEBUG_PRESSURE
      iout << "strain rate: " << strainRate << "\n";
#endif
    }
    if (simParams->fixCellDims) {
      if (simParams->fixCellDimX) strainRate.xx = 0;
      if (simParams->fixCellDimY) strainRate.yy = 0;
      if (simParams->fixCellDimZ) strainRate.zz = 0;
    }

  }

}

void Controller::langevinPiston2(int step)
{
  if ( simParams->langevinPistonOn )
  {
    Tensor &strainRate = langevinPiston_strainRate;
    int cellDims = simParams->useFlexibleCell ? 1 : 3;
    BigReal dt = simParams->dt;
    BigReal dt_long = slowFreq * dt;
    BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
    BigReal tau = simParams->langevinPistonPeriod;
    BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;

    // corrections to integrator
    if ( ! ( step % nbondFreq ) )
    {
#ifdef DEBUG_PRESSURE
      iout << iINFO << "correcting strain rate for nbond, ";
#endif
      strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                ( (nbondFreq - 1) * controlPressure_nbond );
#ifdef DEBUG_PRESSURE
      iout << "strain rate: " << strainRate << "\n";
#endif
    }
    if ( ! ( step % slowFreq ) )
    {
#ifdef DEBUG_PRESSURE
      iout << iINFO << "correcting strain rate for slow, ";
#endif
      strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                ( (slowFreq - 1) * controlPressure_slow );
#ifdef DEBUG_PRESSURE
      iout << "strain rate: " << strainRate << "\n";
#endif
    }

    // Apply surface tension.  If surfaceTensionTarget is zero, we get
    // the default (isotropic pressure) case.
   
    Tensor ptarget;
    ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
    if ( simParams->surfaceTensionTarget != 0. ) {
      ptarget.xx = ptarget.yy = simParams->langevinPistonTarget - 
        simParams->surfaceTensionTarget / state->lattice.c().z;
    }

    strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
      ( controlPressure - ptarget );

 
#ifdef DEBUG_PRESSURE
    iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
#endif

    if ( ! ( step % slowFreq ) )
    {
      BigReal gamma = 1 / simParams->langevinPistonDecay;
      BigReal f1 = exp( -0.5 * dt_long * gamma );
      BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
      strainRate *= f1;
      if ( simParams->useFlexibleCell ) {
        // We only use on-diagonal terms (for now)
        if ( simParams->useConstantRatio ) {
          BigReal r = f2 * random->gaussian();
          strainRate.xx += r;
          strainRate.yy += r;
          strainRate.zz += f2 * random->gaussian();
        } else {
          strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
        }
      } else {
        strainRate += f2 * Tensor::identity(random->gaussian());
      }
#ifdef DEBUG_PRESSURE
      iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
#endif
    }

#ifdef DEBUG_PRESSURE
    iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
#endif
    if (simParams->fixCellDims) {
      if (simParams->fixCellDimX) strainRate.xx = 0;
      if (simParams->fixCellDimY) strainRate.yy = 0;
      if (simParams->fixCellDimZ) strainRate.zz = 0;
    }
  }


}

void Controller::rescaleVelocities(int step)
{
  const int rescaleFreq = simParams->rescaleFreq;
  if ( rescaleFreq > 0 ) {
    rescaleVelocities_sumTemps += temperature;  ++rescaleVelocities_numTemps;
    if ( rescaleVelocities_numTemps == rescaleFreq ) {
      BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
      BigReal rescaleTemp = simParams->rescaleTemp;
      if ( simParams->adaptTempOn && simParams->adaptTempRescale && (step > simParams->adaptTempFirstStep ) && 
                (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
        rescaleTemp = adaptTempT;
      }
      BigReal factor = sqrt(rescaleTemp/avgTemp);
      broadcast->velocityRescaleFactor.publish(step,factor);
      //iout << "RESCALING VELOCITIES AT STEP " << step
      //     << " FROM AVERAGE TEMPERATURE OF " << avgTemp
      //     << " TO " << rescaleTemp << " KELVIN.\n" << endi;
      rescaleVelocities_sumTemps = 0;  rescaleVelocities_numTemps = 0;
    }
  }
}

void Controller::correctMomentum(int step) {
    RequireReduction* reduction = getCurrentReduction();

    Vector momentum;
    BigReal mass;
    momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
    momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
    momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
    mass = reduction->item(REDUCTION_MOMENTUM_MASS);

    if ( momentum.length2() == 0. )
      iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
    if ( mass == 0. )
      NAMD_die("Total mass is zero in Controller::correctMomentum");

    momentum *= (-1./mass);
    // What does this mean???
    broadcast->momentumCorrection.publish(step+slowFreq,momentum);
}

//Modifications for alchemical fep
void Controller::printFepMessage(int step)
{
  if (simParams->alchOn && simParams->alchFepOn 
      && !simParams->alchLambdaFreq) {
    const BigReal alchLambda = simParams->alchLambda;
    const BigReal alchLambda2 = simParams->alchLambda2;
    const BigReal alchLambdaIDWS = simParams->alchLambdaIDWS;
    const BigReal alchTemp = simParams->alchTemp;
    const int alchEquilSteps = simParams->alchEquilSteps;
    iout << "FEP: RESETTING FOR NEW FEP WINDOW "
         << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2;
    if ( alchLambdaIDWS >= 0. ) {
      iout << " LAMBDA_IDWS " << alchLambdaIDWS;
    }
    iout << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
         << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
         << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp 
         << " K FOR FEP CALCULATION\n" << endi;
  }
} 
void Controller::printTiMessage(int step)
{
  if (simParams->alchOn && simParams->alchThermIntOn 
      && !simParams->alchLambdaFreq) {
    const BigReal alchLambda = simParams->alchLambda;
    const int alchEquilSteps = simParams->alchEquilSteps;
    iout << "TI: RESETTING FOR NEW WINDOW "
         << "LAMBDA SET TO " << alchLambda 
         << "\nTI: WINDOW TO HAVE " << alchEquilSteps 
         << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
  }
} 
//fepe

void Controller::reassignVelocities(int step)
{
  const int reassignFreq = simParams->reassignFreq;
  if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
    BigReal newTemp = simParams->reassignTemp;
    newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
    if ( simParams->reassignIncr > 0.0 ) {
      if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
        newTemp = simParams->reassignHold;
    } else {
      if ( newTemp < simParams->reassignHold )
        newTemp = simParams->reassignHold;
    }
    iout << "REASSIGNING VELOCITIES AT STEP " << step
         << " TO " << newTemp << " KELVIN.\n" << endi;
  }
}

void Controller::tcoupleVelocities(int step)
{
  if ( simParams->tCoupleOn )
  {
    const BigReal tCoupleTemp = simParams->tCoupleTemp;
    BigReal coefficient = 1.;
    if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
    broadcast->tcoupleCoefficient.publish(step,coefficient);
  }
}

void Controller::stochRescaleVelocities(int step)
{
  if ( simParams->stochRescaleOn ) {
    ++stochRescale_count;
    if ( stochRescale_count == simParams->stochRescaleFreq ) {
      double coefficient = stochRescaleCoefficient();
      broadcast->stochRescaleCoefficient.publish(step,coefficient);
      stochRescale_count = 0;
    }
  }
}

double Controller::stochRescaleCoefficient() {
  const double stochRescaleTemp = simParams->stochRescaleTemp;
  double coefficient = 1;
  if ( temperature > 0 ) {
    double R1 = random->gaussian();
    // double gammaShape = 0.5*(numDegFreedom - 1);
    // R2sum is the sum of (numDegFreedom - 1) squared normal variables,
    // which is chi-squared distributed.
    // This is in turn a special case of the Gamma distribution,
    // which converges to a normal distribution in the limit of a
    // large shape parameter.
    // double R2sum = 2*(gammaShape+sqrt(gammaShape)*random->gaussian())+R1*R1;
    double R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
    double tempfactor = stochRescaleTemp/(temperature*numDegFreedom);

    coefficient = sqrt(stochRescaleTimefactor +
        (1 - stochRescaleTimefactor)*tempfactor*R2sum +
        2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*
          stochRescaleTimefactor)); 
  }
  heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient-1);
  return coefficient;
}

static char *FORMAT(BigReal X, int decimal = 4)
{
  static char tmp_string[50];
  static char format_string[50];
  const double maxnum = 99999999999.9999;
  if ( X > maxnum ) X = maxnum;
  if ( X < -maxnum ) X = -maxnum;

  int whole = (decimal <= 4 ? 14 : 10 + decimal);
  sprintf(format_string, " %%%d.%df", whole, decimal);
  sprintf(tmp_string, format_string, X); 
  return tmp_string;
}

static char *FORMAT(const char *X, int decimal = 4)
{
  static char tmp_string[50];
  static char format_string[50];

  int width = (decimal <= 4 ? 14 : 10 + decimal);
  sprintf(format_string, " %%%ds", width);
  sprintf(tmp_string, format_string, X); 
  return tmp_string;
}

static char *ETITLE(int X)
{
  static char tmp_string[21];
  sprintf(tmp_string,"ENERGY: %7d",X); 
  return  tmp_string;
}

void Controller::receivePressure(int step, int minimize)
{
    Node *node = Node::Object();
    Molecule *molecule = node->molecule;
    Lattice &lattice = state->lattice;
    bool cudaIntegrator = simParams->CUDASOAintegrate;

    RequireReduction* reduction = getCurrentReduction();
    reduction->require();

    Tensor virial_normal;
    Tensor virial_nbond;
    Tensor virial_slow;
#ifdef ALTVIRIAL
    Tensor altVirial_normal;
    Tensor altVirial_nbond;
    Tensor altVirial_slow;
#endif
    Tensor intVirial_normal;
    Tensor intVirial_nbond;
    Tensor intVirial_slow;
    Vector extForce_normal;
    Vector extForce_nbond;
    Vector extForce_slow;

#if 1
    numDegFreedom = molecule->num_deg_freedom();
    int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
    int numFixedGroups = molecule->num_fixed_groups();
    int numFixedAtoms = molecule->num_fixed_atoms();
#endif
#if 0
    int numAtoms = molecule->numAtoms;
    numDegFreedom = 3 * numAtoms;
    int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
    int numFixedAtoms =
        ( simParams->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
    int numLonepairs = molecule->numLonepairs;
    int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
    if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
    if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
    if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
    if ( ! ( numFixedAtoms || molecule->numConstraints
        || simParams->comMove || simParams->langevinOn ) ) {
      numDegFreedom -= 3;
      numGroupDegFreedom -= 3;
    }
    if (simParams->pairInteractionOn) {
      // this doesn't attempt to deal with fixed atoms or constraints
      numDegFreedom = 3 * molecule->numFepInitial;
    }
    int numRigidBonds = molecule->numRigidBonds;
    int numFixedRigidBonds =
        ( simParams->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );

    // numLonepairs is subtracted here because all lonepairs have a rigid bond
    // to oxygen, but all of the LP degrees of freedom are dealt with above
    numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
#endif
    BigReal groupKineticEnergyHalfstep, groupKineticEnergyCentered;
    kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
    kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
    groupKineticEnergyHalfstep = kineticEnergyHalfstep -
      reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
    groupKineticEnergyCentered = kineticEnergyCentered -
      reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);

    BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
                                        / ( numDegFreedom * BOLTZMANN );
    BigReal atomTempCentered = 2.0 * kineticEnergyCentered
                                        / ( numDegFreedom * BOLTZMANN );
    BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
                                        / ( numGroupDegFreedom * BOLTZMANN );
    BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
                                        / ( numGroupDegFreedom * BOLTZMANN );

    /*  test code for comparing different temperatures  
    iout << "TEMPTEST: " << step << " " << 
        atomTempHalfstep << " " <<
        atomTempCentered << " " <<
        groupTempHalfstep << " " <<
        groupTempCentered << "\n" << endi;
  iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
    numGroupDegFreedom << "\n" << endi;
     */
    GET_TENSOR(momentumSqrSum, reduction,REDUCTION_MOMENTUM_SQUARED);
    
    GET_TENSOR(virial_normal, reduction,REDUCTION_VIRIAL_NORMAL);
    GET_TENSOR(virial_nbond,  reduction,REDUCTION_VIRIAL_NBOND);
    GET_TENSOR(virial_slow,   reduction,REDUCTION_VIRIAL_SLOW);

#ifdef ALTVIRIAL
    GET_TENSOR(altVirial_normal, reduction,REDUCTION_ALT_VIRIAL_NORMAL);
    GET_TENSOR(altVirial_nbond,  reduction,REDUCTION_ALT_VIRIAL_NBOND);
    GET_TENSOR(altVirial_slow,   reduction,REDUCTION_ALT_VIRIAL_SLOW);
#endif

    GET_TENSOR(intVirial_normal, reduction,REDUCTION_INT_VIRIAL_NORMAL);
    GET_TENSOR(intVirial_nbond,  reduction,REDUCTION_INT_VIRIAL_NBOND);
    GET_TENSOR(intVirial_slow,   reduction,REDUCTION_INT_VIRIAL_SLOW);

    GET_VECTOR(extForce_normal, reduction,REDUCTION_EXT_FORCE_NORMAL);
    GET_VECTOR(extForce_nbond, reduction,REDUCTION_EXT_FORCE_NBOND);
    GET_VECTOR(extForce_slow, reduction,REDUCTION_EXT_FORCE_SLOW);
    
    // APH NOTE: These four lines are now done in calcPressure()
    // Vector extPosition = lattice.origin();
    // virial_normal -= outer(extForce_normal,extPosition);
    // virial_nbond -= outer(extForce_nbond,extPosition);
    // virial_slow -= outer(extForce_slow,extPosition);

    kineticEnergy = kineticEnergyCentered;
    temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );

    if (simParams->drudeOn) {
      BigReal drudeComKE, drudeBondKE;
      drudeComKE = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
      drudeBondKE = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
      int g_bond = 3 * molecule->numDrudeAtoms;
      int g_com = numDegFreedom - g_bond;
      temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
      drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
    }

    // Calculate number of degrees of freedom (controlNumDegFreedom)
    if ( simParams->useGroupPressure )
    {
      controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
      if ( ! ( numFixedAtoms || molecule->numConstraints
  || simParams->comMove || simParams->langevinOn ) ) {
        controlNumDegFreedom -= 1;
      }
    }
    else
    {
      controlNumDegFreedom = numDegFreedom / 3;
    }
    if (simParams->fixCellDims) {
      if (simParams->fixCellDimX) controlNumDegFreedom -= 1;
      if (simParams->fixCellDimY) controlNumDegFreedom -= 1;
      if (simParams->fixCellDimZ) controlNumDegFreedom -= 1;
    }

    // Calculate pressure tensors using virials
    calcPressure(step, minimize,
      virial_normal, virial_nbond, virial_slow,
      intVirial_normal, intVirial_nbond, intVirial_slow,
      extForce_normal, extForce_nbond, extForce_slow);

#ifdef DEBUG_PRESSURE
    iout << iINFO << "Control pressure = " << controlPressure <<
      " = " << controlPressure_normal << " + " <<
      controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
#endif
    if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParams->accelMDOutFreq) ) {
        iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
             << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
             << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
             << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
             << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
             << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
             << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
             << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
             << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
             << endi;
   }
}

//
// Calculates all pressure tensors using virials
//
// Sets variables:
// pressure, pressure_normal, pressure_nbond, pressure_slow
// groupPressure, groupPressure_normal, groupPressure_nbond, groupPressure_slow
// controlPressure, controlPressure_normal, controlPressure_nbond, controlPressure_slow
// pressure_amd
// 
void Controller::calcPressure(int step, int minimize,
  const Tensor& virial_normal_in, const Tensor& virial_nbond_in, const Tensor& virial_slow_in,
  const Tensor& intVirial_normal, const Tensor& intVirial_nbond, const Tensor& intVirial_slow,
  const Vector& extForce_normal, const Vector& extForce_nbond, const Vector& extForce_slow) {

  Tensor virial_normal = virial_normal_in;
  Tensor virial_nbond = virial_nbond_in;
  Tensor virial_slow = virial_slow_in;

#if 0
  Tensor tmp = virial_normal;
  fprintf(stderr, "VIRIAL_NORMAL  %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
     tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
  
  tmp = virial_nbond;
  fprintf(stderr, "VIRIAL_NBOND  %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
     tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
  
  tmp = virial_slow;
  fprintf(stderr, "VIRIAL_SLOW  %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
     tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
     
  tmp = intVirial_normal;
  fprintf(stderr, "INT_VIRIAL_NORMAL  %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
     tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
  
  tmp = intVirial_nbond;
  fprintf(stderr, "INT_VIRIAL_NBOND  %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
     tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
  
  tmp = intVirial_slow;
  fprintf(stderr, "INT_VIRIAL_SLOW  %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
     tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
  
  Vector ext = extForce_normal;
  fprintf(stderr, "EXT_FORCE_NORMAL  %1.2lf %1.2lf %1.2lf\n",
     ext.x, ext.y, ext.z);
  
  ext = extForce_nbond;
  fprintf(stderr, "EXT_FORCE_NBOND  %1.2lf %1.2lf %1.2lf\n",
     ext.x, ext.y, ext.z);
  
  ext = extForce_slow;
  fprintf(stderr, "EXT_FORCE_SLOW  %1.2lf %1.2lf %1.2lf\n",
     ext.x, ext.y, ext.z);

  //fprintf(stderr, "State Lattice for timestep %d\n", step);
  fprintf(stderr, "LATTICE  %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", 
    state->lattice.a().x, state->lattice.a().y, state->lattice.a().z,
    state->lattice.b().x, state->lattice.b().y, state->lattice.b().z,
    state->lattice.c().x, state->lattice.c().y, state->lattice.c().z);
    
  fprintf(stderr, "VOLUME  %3.10lf\n", state->lattice.volume());
  
#endif


  Node *node = Node::Object();
  Molecule *molecule = node->molecule;
  Lattice &lattice = state->lattice;
  BigReal volume;

  Vector extPosition = lattice.origin();
  virial_normal -= outer(extForce_normal,extPosition);
  virial_nbond -= outer(extForce_nbond,extPosition);
  virial_slow -= outer(extForce_slow,extPosition);

  if ( (volume=lattice.volume()) != 0. )
  {

    if ((simParams->LJcorrection || simParams->LJcorrectionAlt) && volume) {
#ifdef MEM_OPT_VERSION
      NAMD_bug("LJcorrection not supported in memory optimized build.");
#else
      // Apply tail correction to pressure
      BigReal alchLambda = simParams->getCurrentLambda(step);
      virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
#endif
    }


    // kinetic energy component included in virials
    if(simParams->CUDASOAintegrate && step != 0 && !simParams->isMultiTimeStepping()){
      pressure_normal = virial_normal / volume;
      // groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
      if (simParams->accelMDOn) {
        pressure_amd = virial_amd / volume;
        pressure_normal += pressure_amd;
        groupPressure_normal +=  pressure_amd;
      }

      if ( minimize || ! ( step % nbondFreq ) )
      {
        pressure_nbond = virial_nbond / volume;
        // groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
      }

      if ( minimize || ! ( step % slowFreq ) )
      {
        pressure_slow = virial_slow / volume;
        // groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
      }

      pressure = pressure_normal + pressure_nbond + pressure_slow; 
      groupPressure = pressure - intVirial_normal / volume;
      // groupPressure = groupPressure_normal + groupPressure_nbond +
      //       groupPressure_slow;
    }else{
      pressure_normal = virial_normal / volume;

      groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
      if (simParams->accelMDOn) {
        pressure_amd = virial_amd / volume;
        pressure_normal += pressure_amd;
        groupPressure_normal +=  pressure_amd;
      }

      if ( minimize || ! ( step % nbondFreq ) )
      {
        pressure_nbond = virial_nbond / volume;
        groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
      }

      if ( minimize || ! ( step % slowFreq ) )
      {
        pressure_slow = virial_slow / volume;
        groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
      }

      pressure = pressure_normal + pressure_nbond + pressure_slow; 
      groupPressure = groupPressure_normal + groupPressure_nbond +
            groupPressure_slow;
    }
  }
  else
  {
    pressure = Tensor();
    groupPressure = Tensor();
  }

  if ( simParams->useGroupPressure )
  {
    controlPressure_normal = groupPressure_normal;
    controlPressure_nbond = groupPressure_nbond;
    controlPressure_slow = groupPressure_slow;
    controlPressure = groupPressure;
  }
  else
  {
    controlPressure_normal = pressure_normal;
    controlPressure_nbond = pressure_nbond;
    controlPressure_slow = pressure_slow;
    controlPressure = pressure;
  }

  if ( simParams->useFlexibleCell ) {
    // use symmetric pressure to control rotation
    // controlPressure_normal = symmetric(controlPressure_normal);
    // controlPressure_nbond = symmetric(controlPressure_nbond);
    // controlPressure_slow = symmetric(controlPressure_slow);
    // controlPressure = symmetric(controlPressure);
    // only use on-diagonal components for now
    controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
    controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
    controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
    controlPressure = Tensor::diagonal(diagonal(controlPressure));
    if ( simParams->useConstantRatio ) {
#define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
   T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
      AVGXY(controlPressure_normal);
      AVGXY(controlPressure_nbond);
      AVGXY(controlPressure_slow);
      AVGXY(controlPressure);
#undef AVGXY
    }
  } else {
    controlPressure_normal =
  Tensor::identity(trace(controlPressure_normal)/3.);
    controlPressure_nbond =
  Tensor::identity(trace(controlPressure_nbond)/3.);
    controlPressure_slow =
  Tensor::identity(trace(controlPressure_slow)/3.);
    controlPressure =
  Tensor::identity(trace(controlPressure)/3.);
  }
}

inline void Controller :: calc_accelMDG_mean_std
(BigReal testV, int n, 
 BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV){
   BigReal delta; 

   if(testV > *Vmax){
      *Vmax = testV;
   }
   else if(testV < *Vmin){
      *Vmin = testV;
   }

   //mean and std calculated by Online Algorithm
   delta = testV - *Vavg;
   *Vavg += delta / (BigReal)n;
   *M2 += delta * (testV - (*Vavg));

   *sigmaV = sqrt(*M2/n);
}

inline void Controller :: calc_accelMDG_E_k
(int iE, int step, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, 
 BigReal* k0, BigReal* k, BigReal* E, int *iEused, char *warn){
   switch(iE){
      case 2:
         *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
         // if k0 <=0 OR k0 > 1, switch to iE=1 mode
         if(*k0 > 1.)
            *warn |= 1;
         else if(*k0 <= 0.)
            *warn |= 2;
         //else stay in iE=2 mode
         else{
            *E = Vmin + (Vmax-Vmin)/(*k0);
            *iEused = 2;
            break;
         }
      case 1:
         *E = Vmax;
         *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
         if(*k0 > 1.)
            *k0 = 1.;
         *iEused = 1;
         break;
   }

   *k = *k0 / (Vmax - Vmin);
}

inline void Controller :: calc_accelMDG_force_factor
(BigReal k, BigReal E, BigReal testV, Tensor vir_orig, 
 BigReal *dV, BigReal *factor, Tensor *vir){
   BigReal VE = testV - E;
   //if V < E, apply boost
   if(VE < 0){
      *factor = k * VE;
      *vir += vir_orig * (*factor);
      *dV += (*factor) * VE/2.;
      *factor += 1.;
   }
   else{
      *factor = 1.;
   }
}

void Controller :: write_accelMDG_rest_file
(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){
   FILE *rest;
   char timestepstr[20];
   char *filename;
   int baselen;
   char *restart_name;
   const char *bsuffix;

   if(lasttime || simParams->restartFrequency == 0){
           filename = simParams->outputFilename;
           bsuffix = ".BAK";
   }
   else{
           filename = simParams->restartFilename;
           bsuffix = ".old";
   }

   baselen = strlen(filename);
   restart_name = new char[baselen+26];

   strcpy(restart_name, filename);
   if ( simParams->restartSave && !lasttime) {
      sprintf(timestepstr,".%d",step_n);
      strcat(restart_name, timestepstr);
      bsuffix = ".BAK";
   }
   strcat(restart_name, ".gamd");

   if(write_topic){
      NAMD_backup_file(restart_name,bsuffix);

      rest = fopen(restart_name, "w");
      if(rest){
         fprintf(rest, "# NAMD accelMDG restart file\n"
               "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
               "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
               type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
         fclose(rest);
         iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
      }
      else
         iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
   }
   else{
      rest = fopen(restart_name, "a");
      if(rest){
         fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
                          type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
         fclose(rest);
      }
      else
         iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
   }

   delete[] restart_name;
}


void Controller::rescaleaccelMD(int step, int minimize)
{
    if ( !simParams->accelMDOn ) return;

    amd_reduction->require();

    // copy all to submit_reduction
    for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
      submit_reduction->item(i) += amd_reduction->item(i);
    }
    submit_reduction->submit();

    if (step == simParams->firstTimestep) {
        accelMDdVAverage = 0;
        accelMDdV = 0;
    }
//    if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
    if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;

    Node *node = Node::Object();
    Molecule *molecule = node->molecule;
    Lattice &lattice = state->lattice;

    const BigReal accelMDE = simParams->accelMDE;
    const BigReal accelMDalpha = simParams->accelMDalpha;
    const BigReal accelMDTE = simParams->accelMDTE;
    const BigReal accelMDTalpha = simParams->accelMDTalpha;
    const int accelMDOutFreq = simParams->accelMDOutFreq;

    //GaMD
    static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
    static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
    static BigReal k0P, kP, EP;
    static BigReal k0D, kD, ED;
    static int V_n = 1, iEusedD, iEusedP;
    static char warnD, warnP;
    static int restfreq;

    const int firststep = (simParams->accelMDFirstStep > simParams->firstTimestep) ? simParams->accelMDFirstStep : simParams->firstTimestep;
    const int ntcmdprep = (simParams->accelMDGcMDPrepSteps > 0) ? simParams->accelMDFirstStep + simParams->accelMDGcMDPrepSteps : 0;
    const int ntcmd     = simParams->accelMDFirstStep + simParams->accelMDGcMDSteps;
    const int ntebprep  = ntcmd + simParams->accelMDGEquiPrepSteps;
    const int nteb      = ntcmd + simParams->accelMDGEquiSteps;
    const int ntave     = simParams->accelMDGStatWindow;

    BigReal volume;
    BigReal bondEnergy;
    BigReal angleEnergy;
    BigReal dihedralEnergy;
    BigReal improperEnergy;
    BigReal crosstermEnergy;
    BigReal boundaryEnergy;
    BigReal miscEnergy;
    BigReal amd_electEnergy;
    BigReal amd_ljEnergy;
    BigReal amd_ljEnergy_Corr = 0.;
    BigReal amd_electEnergySlow;
    BigReal amd_groLJEnergy;
    BigReal amd_groGaussEnergy;
    BigReal amd_goTotalEnergy;
    BigReal potentialEnergy;
    BigReal factor_dihe = 1;
    BigReal factor_tot = 1;
    BigReal testV;
    BigReal dV = 0.;
    Vector  accelMDfactor;
    Tensor vir; //auto initialized to 0
    Tensor vir_dihe;
    Tensor vir_normal;
    Tensor vir_nbond;
    Tensor vir_slow;

    volume = lattice.volume();

    bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
    angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
    dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
    improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
    crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
    boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
    miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);

    GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
    GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
    GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
    GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);

    if ( !( step % nbondFreq ) ) {
      amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
      amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
      amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
      amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
      BigReal goNativeEnergy = amd_reduction->item(REDUCTION_GO_NATIVE_ENERGY);
      BigReal goNonnativeEnergy = amd_reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
      amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
    } else {
      amd_electEnergy = electEnergy;
      amd_ljEnergy = ljEnergy;
      amd_groLJEnergy = groLJEnergy;
      amd_groGaussEnergy = groGaussEnergy;
      amd_goTotalEnergy = goTotalEnergy;
    }

    if( (simParams->LJcorrection || simParams->LJcorrectionAlt) && volume ) {
        // Obtain tail correction (copied from printEnergies())
        // This value is only printed for sanity check
        // accelMD currently does not 'boost' tail correction
        amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
    }

    if ( !( step % slowFreq ) ) {
      amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
    } else {
      amd_electEnergySlow = electEnergySlow;
    }

    potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
        improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
        crosstermEnergy + boundaryEnergy + miscEnergy +
        amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;

    //GaMD
    if(simParams->accelMDG){
       // if first time running accelMDG module
       if(step == firststep) {
          iEusedD = iEusedP = simParams->accelMDGiE;
          warnD   = warnP   = '\0';

          //restart file reading
          if(simParams->accelMDGRestart){
             FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
             char line[256];
             int dihe_n=0, tot_n=0;
             if(!rest){
                sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
                NAMD_die(line);
             }

             while(fgets(line, 256, rest)){
                if(line[0] == 'D'){
                   dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
                                   &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
                }
                else if(line[0] == 'T'){
                   tot_n  = sscanf(line+1, " %d %la %la %la %la %la %la %la",
                                   &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
                }
             }

             fclose(rest);

             V_n++;

             //check dihe settings
             if(simParams->accelMDdihe || simParams->accelMDdual){
                if(dihe_n !=8)
                   NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
                k0D = kD * (VmaxD - VminD);
                iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
                   << " Vmax " << VmaxD << " Vmin " << VminD 
                   << " Vavg " << VavgD << " sigmaV " << sigmaVD
                   << " E " << ED << " k " << kD
                   << "\n" << endi;
             }

             //check tot settings
             if(!simParams->accelMDdihe || simParams->accelMDdual){
                if(tot_n !=8)
                   NAMD_die("Invalid total potential energy format in accelMDG restart file");
                k0P = kP * (VmaxP - VminP);
                iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
                   << " Vmax " << VmaxP << " Vmin " << VminP 
                   << " Vavg " << VavgP << " sigmaV " << sigmaVP
                   << " E " << EP << " k " << kP
                   << "\n" << endi;
            }

            iEusedD = (ED == VmaxD) ? 1 : 2;
            iEusedP = (EP == VmaxP) ? 1 : 2;
          }
          //local restfreq follows NAMD restartfreq (default: 500)
          restfreq = simParams->restartFrequency ? simParams->restartFrequency : 500;
       }

       //for dihedral potential
       if(simParams->accelMDdihe || simParams->accelMDdual){
          testV = dihedralEnergy + crosstermEnergy;

          //write restart file every restartfreq steps
          if(step > firststep && step % restfreq == 0)
             write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
                   true, false);
          //write restart file at the end of the simulation
          if( simParams->accelMDLastStep > 0 ){
              if( step == simParams->accelMDLastStep )
                  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
                          true, true);
          }
          else if(step == simParams->N)
              write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
                      true, true);

          //conventional MD
          if(step < ntcmd){
             //very first step
             if(step == firststep && !simParams->accelMDGRestart){
                //initialize stat
                VmaxD = VminD = VavgD = testV;
                M2D = sigmaVD = 0.;
             }
             //first step after cmdprep
             else if(step == ntcmdprep && ntcmdprep != 0){
                //reset stat
                VmaxD = VminD = VavgD = testV;
                M2D = sigmaVD = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
             }
             //every ntave steps
             else if(ntave > 0 && step % ntave == 0){
                //update Vmax Vmin
                if(testV > VmaxD) VmaxD = testV;
                if(testV < VminD) VminD = testV;
                //reset avg and std
                VavgD = testV;
                M2D = sigmaVD = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
             }
             //normal steps
             else
                calc_accelMDG_mean_std(testV, V_n,
                      &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;

             //last cmd step
             if(step == ntcmd - 1){
                calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D, 
                      VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
             }
          }
          //equilibration
          else if(step < nteb){
             calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
                   &dV, &factor_dihe, &vir);

             //first step after cmd
             if(step == ntcmd && simParams->accelMDGresetVaftercmd){
                //reset stat
                VmaxD = VminD = VavgD = testV;
                M2D = sigmaVD = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
             }
             //every ntave steps
             else if(ntave > 0 && step % ntave == 0){
                //update Vmax Vmin
                if(testV > VmaxD) VmaxD = testV;
                if(testV < VminD) VminD = testV;
                //reset avg and std
                VavgD = testV;
                M2D = sigmaVD = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
             }
             else
                calc_accelMDG_mean_std(testV, V_n, 
                      &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;

             //steps after ebprep
             if(step >= ntebprep)
                if((ntave > 0 && step % ntave == 0) || ntave <= 0)
                    calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D,
                          VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
          }
          //production
          else{
             calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
                   &dV, &factor_dihe, &vir);
          }

       }
       //for total potential
       if(!simParams->accelMDdihe || simParams->accelMDdual){
          Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
          testV = potentialEnergy;
          if(simParams->accelMDdual){
             testV -= dihedralEnergy + crosstermEnergy;
             vir_tot -= vir_dihe;
          }

          //write restart file every restartfreq steps
          if(step > firststep && step % restfreq == 0)
             write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
                   !simParams->accelMDdual, false);
          //write restart file at the end of simulation
          if( simParams->accelMDLastStep > 0 ){
              if( step == simParams->accelMDLastStep )
                  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
                          !simParams->accelMDdual, true);
          }
          else if(step == simParams->N)
             write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
                   !simParams->accelMDdual, true);

          //conventional MD
          if(step < ntcmd){
             //very first step
             if(step == firststep && !simParams->accelMDGRestart){
                //initialize stat
                VmaxP = VminP = VavgP = testV;
                M2P = sigmaVP = 0.;
             }
             //first step after cmdprep
             else if(step == ntcmdprep && ntcmdprep != 0){
                //reset stat
                VmaxP = VminP = VavgP = testV;
                M2P = sigmaVP = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
             }
             //every ntave steps
             else if(ntave > 0 && step % ntave == 0){
                //update Vmax Vmin
                if(testV > VmaxP) VmaxP = testV;
                if(testV < VminP) VminP = testV;
                //reset avg and std
                VavgP = testV;
                M2P = sigmaVP = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
             }
             //normal steps
             else
                calc_accelMDG_mean_std(testV, V_n,
                      &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
             //last cmd step
             if(step == ntcmd - 1){
                calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P, 
                      VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
             }
          }
          //equilibration
          else if(step < nteb){
             calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
                   &dV, &factor_tot, &vir);

             //first step after cmd
             if(step == ntcmd && simParams->accelMDGresetVaftercmd){
                //reset stat
                VmaxP = VminP = VavgP = testV;
                M2P = sigmaVP = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
             }
             //every ntave steps
             else if(ntave > 0 && step % ntave == 0){
                //update Vmax Vmin
                if(testV > VmaxP) VmaxP = testV;
                if(testV < VminP) VminP = testV;
                //reset avg and std
                VavgP = testV;
                M2P = sigmaVP = 0.;
                iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
             }
             else
                calc_accelMDG_mean_std(testV, V_n,
                      &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);

             //steps after ebprep
             if(step >= ntebprep)
                if((ntave > 0 && step % ntave == 0) || ntave <= 0)
                    calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P,
                          VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
          }
          //production
          else{
             calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
                   &dV, &factor_tot, &vir);
          }

       }
       accelMDdVAverage += dV;

       //first step after ntcmdprep
       if((ntcmdprep > 0 && step == ntcmdprep) || 
          (step < nteb && ntave > 0 && step % ntave == 0) ||
          (simParams->accelMDGresetVaftercmd && step == ntcmd)){
          V_n = 1;
       }

       if(step < nteb)
           V_n++;

    }
    //aMD
    else{
        if (simParams->accelMDdihe) {

            testV = dihedralEnergy + crosstermEnergy;
            if ( testV < accelMDE ) {
                factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
                factor_dihe *= factor_dihe;
                vir = vir_dihe * (factor_dihe - 1.0);
                dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
                accelMDdVAverage += dV;
            }  

        } else if (simParams->accelMDdual) {

            testV = dihedralEnergy + crosstermEnergy;
            if ( testV < accelMDE ) {
                factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
                factor_dihe *= factor_dihe;
                vir = vir_dihe * (factor_dihe - 1.0) ;
                dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
            }

            testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
            if ( testV < accelMDTE ) {
                factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
                factor_tot *= factor_tot;
                vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
                dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
            }
            accelMDdVAverage += dV;

        } else {

            testV = potentialEnergy;
            if ( testV < accelMDE ) {
                factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
                factor_tot *= factor_tot;
                vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
                dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
                accelMDdVAverage += dV;
            }
        } 
    }

    accelMDfactor[0]=factor_dihe;
    accelMDfactor[1]=factor_tot;
    accelMDfactor[2]=1;
    broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
    virial_amd = vir;
    accelMDdV = dV;

    if ( factor_tot < 0.001 ) {
       iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
            << "\n" << endi;
    }

    if ( ! (step % accelMDOutFreq) ) {
        if ( !(step == simParams->firstTimestep) ) {
             accelMDdVAverage = accelMDdVAverage/accelMDOutFreq; 
        }
        iout << "ACCELERATED MD: STEP " << step
             << " dV "   << dV
             << " dVAVG " << accelMDdVAverage 
             << " BOND " << bondEnergy
             << " ANGLE " << angleEnergy
             << " DIHED " << dihedralEnergy+crosstermEnergy
             << " IMPRP " << improperEnergy
             << " ELECT " << amd_electEnergy+amd_electEnergySlow
             << " VDW "  << amd_ljEnergy
             << " POTENTIAL "  << potentialEnergy;
        if(amd_ljEnergy_Corr)
            iout << " LJcorr " << amd_ljEnergy_Corr;
        iout << "\n" << endi;
        if(simParams->accelMDG){
            if(simParams->accelMDdihe || simParams->accelMDdual)
                iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD 
                    << " Vmax " << VmaxD << " Vmin " << VminD 
                    << " Vavg " << VavgD << " sigmaV " << sigmaVD
                    << " E " << ED << " k0 " << k0D << " k " << kD
                    << "\n" << endi;
            if(warnD & 1)
                iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
                    << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
            if(warnD & 2)
                iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
                    << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
            if(!simParams->accelMDdihe || simParams->accelMDdual)
                iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP 
                    << " Vmax " << VmaxP << " Vmin " << VminP 
                    << " Vavg " << VavgP << " sigmaV " << sigmaVP
                    << " E " << EP << " k0 " << k0P << " k " << kP
                    << "\n" << endi;
            if(warnP & 1)
                iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
                    << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
            if(warnP & 2)
                iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
                    << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
            warnD = warnP = '\0';
        }

        accelMDdVAverage = 0;

        if (simParams->accelMDDebugOn) {
           Tensor p_normal;
           Tensor p_nbond;
           Tensor p_slow;
           Tensor p;
           if ( volume != 0. )  {
                 p_normal = vir_normal/volume;
                 p_nbond  = vir_nbond/volume;
                 p_slow = vir_slow/volume;
                 p = vir/volume;
           }    
           iout << " accelMD Scaling Factor: " << accelMDfactor << "\n" 
                << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
                << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
                << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
                << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n" 
                << endi;
        }
   }
}

void Controller::adaptTempInit(int step) {
    if (!simParams->adaptTempOn) return;
    iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
    adaptTempDtMin = 0;
    adaptTempDtMax = 0;
    adaptTempAutoDt = false;
    if (simParams->adaptTempBins == 0) {
      iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
      std::ifstream adaptTempRead(simParams->adaptTempInFile);
      if (adaptTempRead) {
        int readInt;
        BigReal readReal;
        bool readBool;
        // step
        adaptTempRead >> readInt;
        // Start with min and max temperatures
        adaptTempRead >> adaptTempT;     // KELVIN
        adaptTempRead >> adaptTempBetaMin;  // KELVIN
        adaptTempRead >> adaptTempBetaMax;  // KELVIN
        adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
        adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
        // In case file is manually edited
        if (adaptTempBetaMin > adaptTempBetaMax){
            readReal = adaptTempBetaMax;
            adaptTempBetaMax = adaptTempBetaMin;
            adaptTempBetaMin = adaptTempBetaMax;
        }
        adaptTempRead >> adaptTempBins;     
        adaptTempRead >> adaptTempCg;
        adaptTempRead >> adaptTempDt;
        adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
        adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
        adaptTempPotEnergySamples = new int[adaptTempBins];
        adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
        adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
        adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
        adaptTempBetaN           = new BigReal[adaptTempBins];
        adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
        for(int j = 0; j < adaptTempBins; ++j) {
          adaptTempRead >> adaptTempPotEnergyAve[j];
          adaptTempRead >> adaptTempPotEnergyVar[j];
          adaptTempRead >> adaptTempPotEnergySamples[j];
          adaptTempRead >> adaptTempPotEnergyAveNum[j];
          adaptTempRead >> adaptTempPotEnergyVarNum[j];
          adaptTempRead >> adaptTempPotEnergyAveDen[j];
          adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
        } 
        adaptTempRead.close();
      }
      else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
    } 
    else {
      adaptTempBins = simParams->adaptTempBins;
      adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
      adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
      adaptTempPotEnergySamples = new int[adaptTempBins];
      adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
      adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
      adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
      adaptTempBetaN           = new BigReal[adaptTempBins];
      adaptTempBetaMax = 1./simParams->adaptTempTmin;
      adaptTempBetaMin = 1./simParams->adaptTempTmax;
      adaptTempCg = simParams->adaptTempCgamma;   
      adaptTempDt  = simParams->adaptTempDt;
      adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
      adaptTempT = simParams->initialTemp; 
      if (simParams->langevinOn)
        adaptTempT = simParams->langevinTemp;
      else if (simParams->rescaleFreq > 0)
        adaptTempT = simParams->rescaleTemp;
      for(int j = 0; j < adaptTempBins; ++j){
          adaptTempPotEnergyAveNum[j] = 0.;
          adaptTempPotEnergyAveDen[j] = 0.;
          adaptTempPotEnergySamples[j] = 0;
          adaptTempPotEnergyVarNum[j] = 0.;
          adaptTempPotEnergyVar[j] = 0.;
          adaptTempPotEnergyAve[j] = 0.;
          adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
      }
    }
    if (simParams->adaptTempAutoDt > 0.0) {
       adaptTempAutoDt = true;
       adaptTempDtMin =  simParams->adaptTempAutoDt - 0.01;
       adaptTempDtMax =  simParams->adaptTempAutoDt + 0.01;
    }
    adaptTempDTave = 0;
    adaptTempDTavenum = 0;
    iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
    iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
    iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
    iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
    iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
    if (simParams->adaptTempRestartFreq > 0) {
      NAMD_backup_file(simParams->adaptTempRestartFile);
      adaptTempRestartFile.open(simParams->adaptTempRestartFile);
       if(!adaptTempRestartFile)
        NAMD_die("Error opening restart file");
    }
}

void Controller::adaptTempWriteRestart(int step) {
    if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
        adaptTempRestartFile.seekp(std::ios::beg);        
        iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
        adaptTempRestartFile << step << " ";
        // Start with min and max temperatures
        adaptTempRestartFile << adaptTempT<< " ";     // KELVIN
        adaptTempRestartFile << 1./adaptTempBetaMin << " ";  // KELVIN
        adaptTempRestartFile << 1./adaptTempBetaMax << " ";  // KELVIN
        adaptTempRestartFile << adaptTempBins << " ";     
        adaptTempRestartFile << adaptTempCg << " ";
        adaptTempRestartFile << adaptTempDt ;
        adaptTempRestartFile << "\n" ;
        for(int j = 0; j < adaptTempBins; ++j) {
          adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
          adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
          adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
          adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
          adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
          adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
          adaptTempRestartFile << "\n";          
        }
        adaptTempRestartFile.flush(); 
    }
}    

void Controller::adaptTempUpdate(int step, int minimize)
{
    //Beta = 1./T
    if ( !simParams->adaptTempOn ) return;
    int j = 0;
    if (step == simParams->firstTimestep) {
        adaptTempInit(step);
        return;
    }
    if ( minimize || (step < simParams->adaptTempFirstStep ) || 
        ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
    const int adaptTempOutFreq  = simParams->adaptTempOutFreq;
    const bool adaptTempDebug  = simParams->adaptTempDebug;
    //Calculate Current inverse temperature and bin 
    BigReal adaptTempBeta = 1./adaptTempT;
    adaptTempBin   = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);

    if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
        iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin  
                               << " adaptTempBeta: " << adaptTempBeta 
                              << " adaptTempDBeta: " << adaptTempDBeta 
                               << " betaMin:" << adaptTempBetaMin 
                               << " betaMax: " << adaptTempBetaMax << "\n";
    adaptTempPotEnergySamples[adaptTempBin] += 1;
    BigReal gammaAve = 1.-adaptTempCg/adaptTempPotEnergySamples[adaptTempBin];

    BigReal potentialEnergy;
    BigReal potEnergyAverage;
    BigReal potEnergyVariance;
    potentialEnergy = totalEnergy-kineticEnergy;

    //calculate new bin average and variance using adaptive averaging
    adaptTempPotEnergyAveNum[adaptTempBin] = adaptTempPotEnergyAveNum[adaptTempBin]*gammaAve + potentialEnergy;
    adaptTempPotEnergyAveDen[adaptTempBin] = adaptTempPotEnergyAveDen[adaptTempBin]*gammaAve + 1;
    adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
    
    potEnergyAverage = adaptTempPotEnergyAveNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
    potEnergyVariance = adaptTempPotEnergyVarNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
    potEnergyVariance -= potEnergyAverage*potEnergyAverage;

    adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
    adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
    
    // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
    // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
    // is constant for each bin. This is to estimate <E(beta)> where beta \in
    // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
    if ( ! ( step % simParams->adaptTempFreq ) ) {
      // If adaptTempBin not at the edge, calculate improved average:
      if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
          // Get Averaging Limits:
          BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
          BigReal betaPlus;
          BigReal betaMinus;
          int     nMinus =0;
          int     nPlus = 0;
          if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
          if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
          betaMinus = adaptTempBeta - deltaBeta;
          betaPlus =  adaptTempBeta + deltaBeta;
          BigReal betaMinus2 = betaMinus*betaMinus;
          BigReal betaPlus2  = betaPlus*betaPlus;
          nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
          nPlus  = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
          // Variables for <E(beta)> estimate:
          BigReal potEnergyAve0 = 0.0;
          BigReal potEnergyAve1 = 0.0;
          // Integral terms
          BigReal A0 = 0;
          BigReal A1 = 0;
          BigReal A2 = 0;
          //A0 phi_s integral for beta_minus < beta < beta_{i+1}
          BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1]; 
          BigReal DeltaE2Ave;
          BigReal DeltaE2AveSum = 0;
          BigReal betaj;
          for (j = nMinus+1; j <= adaptTempBin; ++j) {
              potEnergyAve0 += adaptTempPotEnergyAve[j]; 
              DeltaE2Ave = adaptTempPotEnergyVar[j];
              DeltaE2AveSum += DeltaE2Ave;
              A0 += j*DeltaE2Ave;
          }
          A0 *= adaptTempDBeta;
          A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
          A0 *= adaptTempDBeta;
          betaj = adaptTempBetaN[nMinus+1]-betaMinus; 
          betaj *= betaj;
          A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
          A0 /= (betaNp1-betaMinus);

          //A1 phi_s integral for beta_{i+1} < beta < beta_plus
          DeltaE2AveSum = 0;
          for (j = adaptTempBin+1; j < nPlus; j++) {
              potEnergyAve1 += adaptTempPotEnergyAve[j];
              DeltaE2Ave = adaptTempPotEnergyVar[j];
              DeltaE2AveSum += DeltaE2Ave;
              A1 += j*DeltaE2Ave;
          }
          A1 *= adaptTempDBeta;   
          A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
          A1 *= adaptTempDBeta;
          if ( nPlus < adaptTempBins ) {
            betaj = betaPlus-adaptTempBetaN[nPlus];
            betaj *= betaj;
            A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
          }
          A1 /= (betaPlus-betaNp1);

          //A2 phi_t integral for beta_i
          A2 = 0.5*adaptTempDBeta*potEnergyVariance;

          // Now calculate a+ and a-
          DeltaE2Ave = A0-A1;
          BigReal aplus = 0;
          BigReal aminus = 0;
          if (DeltaE2Ave != 0) {
            aplus = (A0-A2)/(A0-A1);
            if (aplus < 0) {
                    aplus = 0;
            }
            if (aplus > 1)  {
                    aplus = 1;
            }
            aminus = 1-aplus;
            potEnergyAve0 *= adaptTempDBeta;
            potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
            potEnergyAve0 /= (betaNp1-betaMinus);
            potEnergyAve1 *= adaptTempDBeta;
            if ( nPlus < adaptTempBins ) {
                potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
            }
            potEnergyAve1 /= (betaPlus-betaNp1);
            potEnergyAverage = aminus*potEnergyAve0;
            potEnergyAverage += aplus*potEnergyAve1;
          }
          if (simParams->adaptTempDebug) {
       iout << "ADAPTEMP DEBUG:"  << "\n"
            << "     adaptTempBin:    " << adaptTempBin << "\n"
            << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
            << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
            << "     adaptTemp:   " << adaptTempT<< "\n"
            << "     betaMin:   " << adaptTempBetaMin << "\n"
            << "     betaMax:   " << adaptTempBetaMax << "\n"
            << "     gammaAve:  " << gammaAve << "\n"
            << "     deltaBeta: " << deltaBeta << "\n"
            << "     betaMinus: " << betaMinus << "\n"
            << "     betaPlus:  " << betaPlus << "\n"
            << "     nMinus:    " << nMinus << "\n"
            << "     nPlus:     " << nPlus << "\n"
            << "     A0:        " << A0 << "\n"
            << "     A1:        " << A1 << "\n"
            << "     A2:        " << A2 << "\n"
            << "     a+:        " << aplus << "\n"
            << "     a-:        " << aminus << "\n"
            << endi;
          }
      }
      else {
          if (simParams->adaptTempDebug) {
       iout << "ADAPTEMP DEBUG:"  << "\n"
            << "     adaptTempBin:    " << adaptTempBin << "\n"
            << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
            << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
            << "     adaptTemp:   " << adaptTempT<< "\n"
            << "     betaMin:   " << adaptTempBetaMin << "\n"
            << "     betaMax:   " << adaptTempBetaMax << "\n"
            << "     gammaAve:  " << gammaAve << "\n"
            << "     deltaBeta:  N/A\n"
            << "     betaMinus:  N/A\n"
            << "     betaPlus:   N/A\n"
            << "     nMinus:     N/A\n"
            << "     nPlus:      N/A\n"
            << "     A0:         N/A\n"
            << "     A1:         N/A\n"
            << "     A2:         N/A\n"
            << "     a+:         N/A\n"
            << "     a-:         N/A\n"
            << endi;
          }
      }
      
      //dT is new temperature
      BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
      dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
      dT += adaptTempT;
      
     // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
     // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
      if ( dT > 1./adaptTempBetaMin || dT  < 1./adaptTempBetaMax ) {
        dT = ((potentialEnergy-adaptTempPotEnergyAve[adaptTempBin])/BOLTZMANN+adaptTempT)*adaptTempDt;
        dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
        dT += adaptTempT;
        // Check again, if not then keep original adaptTempTor assign random.
        if ( dT > 1./adaptTempBetaMin ) {
          if (!simParams->adaptTempRandom) {             
             //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
             //     << " K higher than adaptTempTmax."
             //     << " Keeping temperature at " 
             //     << adaptTempT<< "\n"<< endi;             
             dT = adaptTempT;
          }
          else {             
             //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
             //     << " K higher than adaptTempTmax."
             //     << " Assigning random temperature in range\n" << endi;
             dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);             
             dT = 1./dT;
          }
        } 
        else if ( dT  < 1./adaptTempBetaMax ) {
          if (!simParams->adaptTempRandom) {            
            //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
            //     << " K lower than adaptTempTmin."
            //     << " Keeping temperature at " << adaptTempT<< "\n" << endi; 
            dT = adaptTempT;
          }
          else {
            //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
            //     << " K lower than adaptTempTmin."
            //     << " Assigning random temperature in range\n" << endi;
            dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
            dT = 1./dT;
          }
        }
        else if (adaptTempAutoDt) {
          //update temperature step size counter
          //FOR "TRUE" ADAPTIVE TEMPERING 
          BigReal adaptTempTdiff = fabs(dT-adaptTempT);
          if (adaptTempTdiff > 0) {
            adaptTempDTave += adaptTempTdiff; 
            adaptTempDTavenum++;
//            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
          }
          if(adaptTempDTavenum == 100){
                BigReal Frac;
                adaptTempDTave /= adaptTempDTavenum;
                Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
                Frac = adaptTempDTave/Frac;
                //if average temperature jump is > 3% of temperature range,
                //modify jump size to match 3%
                iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
                if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
                    Frac = adaptTempDtMax/Frac;
                    iout << "ADAPTEMP: Updating adaptTempDt to ";
                    adaptTempDt *= Frac;
                    iout << adaptTempDt << "\n" << endi;
                }
                adaptTempDTave = 0;
                adaptTempDTavenum = 0;
          }
        }
      }
      else if (adaptTempAutoDt) {
          //update temperature step size counter
          // FOR "TRUE" ADAPTIVE TEMPERING
          BigReal adaptTempTdiff = fabs(dT-adaptTempT);
          if (adaptTempTdiff > 0) {
            adaptTempDTave += adaptTempTdiff; 
            adaptTempDTavenum++;
//            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
          }
          if(adaptTempDTavenum == 100){
                BigReal Frac;
                adaptTempDTave /= adaptTempDTavenum;
                Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
                Frac = adaptTempDTave/Frac;
                //if average temperature jump is > 3% of temperature range,
                //modify jump size to match 3%
                iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
                if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
                    Frac = adaptTempDtMax/Frac;
                    iout << "ADAPTEMP: Updating adaptTempDt to ";
                    adaptTempDt *= Frac;
                    iout << adaptTempDt << "\n" << endi;
                }
                adaptTempDTave = 0;
                adaptTempDTavenum = 0;

          }
          
      }
      
      adaptTempT = dT; 
      broadcast->adaptTemperature.publish(step,adaptTempT);
    }
    adaptTempWriteRestart(step);
    if ( ! (step % adaptTempOutFreq) ) {
        iout << "ADAPTEMP: STEP " << step
             << " TEMP "   << adaptTempT
             << " ENERGY " << std::setprecision(10) << potentialEnergy   
             << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
             << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
        iout << "\n" << endi;
   }
   
}


void Controller::compareChecksums(int step, int forgiving) {
    Node *node = Node::Object();
    Molecule *molecule = node->molecule;
    bool cudaIntegrator = simParams->CUDASOAintegrate;
    RequireReduction* reduction = getCurrentReduction();

    // Some basic correctness checking
    BigReal checksum, checksum_b;
    char errmsg[256];

    checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
    if ( ((int)checksum) != molecule->numAtoms && step != 0) {
      sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
              (int)checksum, molecule->numAtoms);
      NAMD_bug(errmsg);
    }
    // do we ever need this if we're on a GPU-based code?
    checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
    
    if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
      if ( ! computeChecksum ) {
        computesPartitioned = 0;
      } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
        computesPartitioned = 1;
      } else {
        NAMD_bug("Bad global compute count!\n");
      }
      computeChecksum = ((int)checksum);
    }
    checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
    if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
      sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcBonds);
      if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }

    checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);

    if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
      sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcAngles);
      if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }
    checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
    if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
      sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcDihedrals);
      if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }

    checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
    if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
      sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcImpropers);
      if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }

    checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
    if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
      sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcTholes);
      if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }

    checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
    if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
      sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcAnisos);
      if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }

    checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
    if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
      sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcCrossterms);
      if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }

    checksum = reduction->item(REDUCTION_ONEFOURNBTHOLE_CHECKSUM);
    if ( checksum_b && (((int)checksum) != molecule->numCalcOneFourNbTholes) ) {
      sprintf(errmsg, "Bad global 1-4 NbThole count! (%d vs %d)\n",
              (int)checksum, molecule->numCalcOneFourNbTholes);
      if ( forgiving && (((int)checksum) < molecule->numCalcOneFourNbTholes) )
        iout << iWARN << errmsg << endi;
      else NAMD_bug(errmsg);
    }

    checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
    if ( ((int64)checksum) > molecule->numCalcExclusions &&
         ( ! simParams->mollyOn || step % slowFreq ) ) {
      if ( forgiving )
        iout << iWARN << "High global exclusion count ("
                      << ((int64)checksum) << " vs "
                      << molecule->numCalcExclusions << "), possible error!\n"
          << iWARN << "This warning is not unusual during minimization.\n"
          << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
      else {
        char errmsg[256];
        sprintf(errmsg, "High global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
                (long long)checksum, (long long)molecule->numCalcExclusions);
        NAMD_bug(errmsg);
      }
    }

    if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
        ! simParams->goGroPair && ! simParams->qmForcesOn) {
      if ( forgiving || ! simParams->fullElectFrequency ) {
        iout << iWARN << "Low global exclusion count!  ("
          << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
        if ( forgiving ) {
          iout << "\n"
            << iWARN << "This warning is not unusual during minimization.\n"
            << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
        } else {
          iout << "  System unstable or pairlistdist or cutoff too small.\n";
        }
        iout << endi;
      } else {
        char errmsg[256];
        sprintf(errmsg, "Low global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
                (long long)checksum, (long long)molecule->numCalcExclusions);
        NAMD_bug(errmsg);
      }
    }

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
    if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
         ( ! simParams->mollyOn || step % slowFreq ) ) {
      if ( forgiving )
        iout << iWARN << "High global CUDA exclusion count ("
                      << ((int64)checksum) << " vs "
                      << molecule->numCalcFullExclusions << "), possible error!\n"
          << iWARN << "This warning is not unusual during minimization.\n"
          << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
      else {
        char errmsg[256];
        sprintf(errmsg, "High global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
                (long long)checksum, (long long)molecule->numCalcFullExclusions);
        NAMD_bug(errmsg);
      }
    }

    if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
        ! simParams->goGroPair && ! simParams->qmForcesOn) {
      if ( forgiving || ! simParams->fullElectFrequency ) {
        iout << iWARN << "Low global CUDA exclusion count!  ("
          << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ") on step " << step;
        if ( forgiving ) {
          iout << "\n"
            << iWARN << "This warning is not unusual during minimization.\n"
            << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
        } else {
          iout << "  System unstable or pairlistdist or cutoff too small.\n";
        }
        iout << endi;
      } else {
        char errmsg[256];
        sprintf(errmsg, "Low global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
                (long long)checksum, (long long)molecule->numCalcFullExclusions);
        NAMD_bug(errmsg);
      }
    }
#endif

    checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
    if ( ((int)checksum) && ! marginViolations ) {
      iout << iERROR << "Margin is too small for " << ((int)checksum) <<
        " atoms during timestep " << step << ".\n" << iERROR <<
      "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
    }
    marginViolations += (int)checksum;

    checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
    if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
      iout << iINFO <<
        "Pairlistdist is too small for " << ((int)checksum) <<
        " computes during timestep " << step << ".\n" << endi;
    }
    if ( simParams->outputPairlists )  pairlistWarnings += (int)checksum;

    checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
    if ( checksum ) {
      if ( forgiving )
        iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
      else NAMD_bug("Stray PME grid charges detected!\n");
    }
}

void Controller::printTiming(int step) {

    if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
    {
      const double endWTime = namdWallTimer() - firstWTime;
      const double endCTime = CmiTimer() - firstCTime;

      // fflush about once per minute
      if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;

      const double elapsedW = 
        (endWTime - startWTime) / simParams->outputTiming;
      const double elapsedC = 
        (endCTime - startCTime) / simParams->outputTiming;

      const double remainingW = elapsedW * (simParams->N - step);
      const double remainingW_hours = remainingW / 3600;

      startWTime = endWTime;
      startCTime = endCTime;

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
      if ( simParams->computeEnergies < 60 &&
           step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
        iout << iWARN << "Energy evaluation is expensive, increase computeEnergies to improve performance.\n" << endi;
        iout << iWARN << "If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n" << endi;
      }
#endif
      if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
        if (simParams->nsPerDayOn) {
          BigReal ns = simParams->dt / 1000000.0;
          BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
          BigReal nsPerDay = ns / (elapsedW * days);
          CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
                    ", %g ns/days"
                    ", %g hours remaining, %f MB of memory in use.\n",
                    step, endCTime, elapsedC, endWTime, elapsedW,
                    nsPerDay,
                    remainingW_hours, memusage_MB());
        }
        else {
          CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
                    ", %g hours remaining, %f MB of memory in use.\n",
                    step, endCTime, elapsedC, endWTime, elapsedW,
                    remainingW_hours, memusage_MB());
        }
        if (simParams->outputPerformance) {
          // calculate running average and sample variance of
          // the measured time per step (elapsedW)
          perfstats.add_sample(elapsedW);
          double t_avg = perfstats.average();
          double t_std = perfstats.standard_deviation();
          // convert t_avg to average number of nanoseconds simulated per day
          double ns_per_day = (simParams->dt / t_avg) * (60 * 60 * 24 * 1e-6);
          CmiPrintf("PERFORMANCE: %d  averaging %g ns/day, %g sec/step with"
              " standard deviation %g\n", step, ns_per_day, t_avg, t_std);
        }
        if ( fflush_count ) { --fflush_count; fflush(stdout); }
        double benchmarkTime = (double) simParams->benchmarkTime;
        if (benchmarkTime > 0 && benchmarkTime <= endWTime) {
          // terminate NAMD benchmark
          char s[100];
          sprintf(s, "Exceeded benchmark time limit %g seconds",
              benchmarkTime);
          NAMD_quit(s);
        }
      }
    }
}

void Controller::printMinimizeEnergies(int step) {

    rescaleaccelMD(step,1);
    receivePressure(step,1);
    RequireReduction* reduction = getCurrentReduction();

    Node *node = Node::Object();
    Molecule *molecule = node->molecule;
    compareChecksums(step,1);

    printEnergies(step,1);

    min_energy = totalEnergy;
    min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
    min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
    min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
    min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));

#ifdef DEBUG_MINIMIZE
    printf("%s, line %d\n", __FILE__, __LINE__);
    printf("Step %d:\n", step);
    printf("   min_energy = %f\n", min_energy);
    printf("   min_f_dot_f = %f\n", min_f_dot_f);
    printf("   min_f_dot_v = %f\n", min_f_dot_v);
    printf("   min_v_dot_v = %f\n", min_v_dot_v);
    printf("   min_huge_count = %d\n", min_huge_count);
    printf("\n");
#endif
}

void Controller::printDynamicsEnergies(int step) {
    compareChecksums(step);
    printEnergies(step,0);
}

BigReal Controller::getTotalPotentialEnergy(int step)
{
    compareChecksums(step, 0);
    
    Node *node = Node::Object();
    Molecule *molecule = node->molecule;
    Lattice &lattice = state->lattice;
    BigReal volume = lattice.volume();
    // store total energy
    BigReal totalPotentialEnergy = 0.0;
    BigReal bondEnergy, angleEnergy, dihedralEnergy;
    BigReal improperEnergy, crosstermEnergy;
    BigReal electEnergy, ljEnergy, electEnergySlow, ljEnergySlow;
    BigReal miscEnergy, bcEnergy;
    bool cudaIntegrator = simParams->CUDASOAintegrate;
    RequireReduction* reduction = getCurrentReduction();

    bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
    angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
    dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
    improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
    crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
    bcEnergy = reduction->item(REDUCTION_BC_ENERGY);
    miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);

    totalPotentialEnergy += bondEnergy + angleEnergy + dihedralEnergy +
                  improperEnergy + crosstermEnergy + miscEnergy
                  + bcEnergy;

    if (! ( step % nbondFreq ))
    {
      electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
      ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
      totalPotentialEnergy += electEnergy + ljEnergy;
    }

    if (! ( step % slowFreq ))
    {
      electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
      ljEnergySlow = reduction->item(REDUCTION_LJ_ENERGY_SLOW);
      totalPotentialEnergy += electEnergySlow;
      totalPotentialEnergy += ljEnergySlow;
    }

    if (simParams->LJcorrection && volume) {
#ifdef MEM_OPT_VERSION
      NAMD_bug("LJcorrection not supported in memory optimized build.");
#else
      // Apply tail correction to energy.
      BigReal alchLambda = simParams->getCurrentLambda(step);
      totalPotentialEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
#endif
    }

  #if 0
    printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, Bond: %f, Angle: %f, Dih: %f, IMP: %f, Cross: %f\n", 
        step, CkMyPe(), bondEnergy, angleEnergy, dihedralEnergy, improperEnergy, crosstermEnergy);
    printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, LJ: %f, Real: %f, Recip: %f\n", step, CkMyPe(), 
        ljEnergy, electEnergy, electEnergySlow);
    printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, MISC: %f, BC: %f, TOTAL: %f\n", step, CkMyPe(), 
        miscEnergy, bcEnergy, totalPotentialEnergy);
  #endif

    return totalPotentialEnergy;
}


void Controller::printEnergies(int step, int minimize)
{
    Node *node = Node::Object();
    Molecule *molecule = node->molecule;
    Lattice &lattice = state->lattice;
    bool cudaIntegrator = simParams->CUDASOAintegrate;
    RequireReduction* reduction = getCurrentReduction();

    // Drude model ANISO energy is added into BOND energy
    // and THOLE energy is added into ELECT energy

    BigReal bondEnergy;
    BigReal angleEnergy;
    BigReal dihedralEnergy;
    BigReal improperEnergy;
    BigReal crosstermEnergy;
    BigReal boundaryEnergy;
    BigReal miscEnergy;
    BigReal potentialEnergy;
    BigReal flatEnergy;
    BigReal smoothEnergy;
    BigReal work;

    Vector momentum;
    Vector angularMomentum;
    BigReal volume = lattice.volume();
    bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
    angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
    dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
    improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
    crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
    boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
    miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);

    if ( minimize || ! ( step % nbondFreq ) )
    {
      electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
      ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);

      // JLai
      groLJEnergy = reduction->item(REDUCTION_GRO_LJ_ENERGY);
      groGaussEnergy = reduction->item(REDUCTION_GRO_GAUSS_ENERGY);

      goNativeEnergy = reduction->item(REDUCTION_GO_NATIVE_ENERGY);
      goNonnativeEnergy = reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
      goTotalEnergy = goNativeEnergy + goNonnativeEnergy;

      //fepb
      bondedEnergyDiff_f = reduction->item(REDUCTION_BONDED_ENERGY_F);
      electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
      ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);

      bondedEnergy_ti_1 = reduction->item(REDUCTION_BONDED_ENERGY_TI_1);
      electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
      ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
      bondedEnergy_ti_2 = reduction->item(REDUCTION_BONDED_ENERGY_TI_2);
      electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
      ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
//fepe
    }

    if ( minimize || ! ( step % slowFreq ) )
    {
      electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
      ljEnergySlow = reduction->item(REDUCTION_LJ_ENERGY_SLOW);
//fepb
      electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);

      electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
      electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
      electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
      electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
//fepe
    }

    if ((simParams->LJcorrection || simParams->LJcorrectionAlt) && volume) {
#ifdef MEM_OPT_VERSION
      NAMD_bug("LJcorrection not supported in memory optimized build.");
#else
      // Apply tail correction to energy.
      BigReal alchLambda = simParams->getCurrentLambda(step);
      ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;

      if (simParams->alchOn) {
        if (simParams->alchFepOn) {
          BigReal alchLambda2 = simParams->getCurrentLambda2(step);
          ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
        } else if (simParams->alchThermIntOn) {
          ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
          ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
        }
      }
#endif
    }

//fepb BKR - Compute alchemical work if using dynamic lambda.  This is here so
//           that the cumulative work can be given during a callback.
    if (simParams->alchLambdaFreq > 0) {
      if (step <= 
          simParams->firstTimestep + simParams->alchEquilSteps) {
        cumAlchWork = 0.0;
      }
      alchWork = computeAlchWork(step);
      cumAlchWork += alchWork;
    }
//fepe
    momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
    momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
    momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
    angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
    angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
    angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);

    // Ported by JLai
    potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
  + improperEnergy + electEnergy + electEnergySlow 
  + ljEnergy + ljEnergySlow
  + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy 
        + groLJEnergy + groGaussEnergy);
    // End of port
    totalEnergy = potentialEnergy + kineticEnergy;
    if ( ! cudaIntegrator) {
      flatEnergy = (totalEnergy
          + (1.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
      if ( !(step%slowFreq) ) {
        // only adjust based on most accurate energies
        BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
        if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
        if ( step != simParams->firstTimestep ) {
          smooth2_avg *= 0.9375;
          smooth2_avg += 0.0625 * s;
        }
      }
      smoothEnergy = (flatEnergy + smooth2_avg
          - (4.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
    }

    // Reset values for accumulated heat and work.
    if (step <= simParams->firstTimestep &&
        (simParams->stochRescaleOn && simParams->stochRescaleHeat)) {
      heat = 0.0;
      totalEnergy0 = totalEnergy;
    }
    if ( simParams->outputMomenta && ! minimize &&
         ! ( step % simParams->outputMomenta ) )
    {
      iout << "MOMENTUM: " << step 
           << " P: " << momentum
           << " L: " << angularMomentum
           << "\n" << endi;
    }

    if ( simParams->outputPressure ) {
      if ( ! cudaIntegrator ) {
        pressure_tavg += pressure;
        groupPressure_tavg += groupPressure;
      }
      tavg_count += 1;
      if ( minimize || ! ( step % simParams->outputPressure ) ) {
        if (cudaIntegrator) {
          // tensors are valid for outputPressure step
          // update averages, then copy into tensor for output
          pressureAverage_xx.addSample(pressure.xx);
          pressureAverage_yx.addSample(pressure.yx);
          pressureAverage_yy.addSample(pressure.yy);
          pressureAverage_zx.addSample(pressure.zx);
          pressureAverage_zy.addSample(pressure.zy);
          pressureAverage_zz.addSample(pressure.zz);
          groupPressureAverage_xx.addSample(groupPressure.xx);
          groupPressureAverage_xy.addSample(groupPressure.xy);
          groupPressureAverage_xz.addSample(groupPressure.xz);
          groupPressureAverage_yx.addSample(groupPressure.yx);
          groupPressureAverage_yy.addSample(groupPressure.yy);
          groupPressureAverage_yz.addSample(groupPressure.yz);
          groupPressureAverage_zx.addSample(groupPressure.zx);
          groupPressureAverage_zy.addSample(groupPressure.zy);
          groupPressureAverage_zz.addSample(groupPressure.zz);
          pressure_tavg.xx = pressureAverage_xx.average();
          pressure_tavg.xy = pressureAverage_yx.average();
          pressure_tavg.xz = pressureAverage_zx.average();
          pressure_tavg.yx = pressureAverage_yx.average();
          pressure_tavg.yy = pressureAverage_yy.average();
          pressure_tavg.yz = pressureAverage_zy.average();
          pressure_tavg.zx = pressureAverage_zx.average();
          pressure_tavg.zy = pressureAverage_zy.average();
          pressure_tavg.zz = pressureAverage_zz.average();
          groupPressure_tavg.xx = groupPressureAverage_xx.average();
          groupPressure_tavg.xy = groupPressureAverage_xy.average();
          groupPressure_tavg.xz = groupPressureAverage_xz.average();
          groupPressure_tavg.yx = groupPressureAverage_yx.average();
          groupPressure_tavg.yy = groupPressureAverage_yy.average();
          groupPressure_tavg.yz = groupPressureAverage_yz.average();
          groupPressure_tavg.zx = groupPressureAverage_zx.average();
          groupPressure_tavg.zy = groupPressureAverage_zy.average();
          groupPressure_tavg.zz = groupPressureAverage_zz.average();
          iout << "PRESSURE: " << step << " "
             << PRESSUREFACTOR * pressure << "\n"
             << "GPRESSURE: " << step << " "
             << PRESSUREFACTOR * groupPressure << "\n";
          // tavg_count makes output behaves the same as for standard case
          // but is no longer part of the average
          if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
             << (PRESSUREFACTOR) * pressure_tavg << "\n"
             << "GPRESSAVG: " << step << " "
             << (PRESSUREFACTOR) * groupPressure_tavg << "\n";
          iout << endi;
          tavg_count = 0;
        }
        else {
          iout << "PRESSURE: " << step << " "
             << PRESSUREFACTOR * pressure << "\n"
             << "GPRESSURE: " << step << " "
             << PRESSUREFACTOR * groupPressure << "\n";
          if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
             << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
             << "GPRESSAVG: " << step << " "
             << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
          iout << endi;
          pressure_tavg = 0;
          groupPressure_tavg = 0;
          tavg_count = 0;
        }
      }
    }

    // pressure profile reductions
    if (pressureProfileSlabs) {
      const int freq = simParams->pressureProfileFreq;
      const int arraysize = 3*pressureProfileSlabs;
      
      BigReal *total = new BigReal[arraysize];
      memset(total, 0, arraysize*sizeof(BigReal));
      const int first = simParams->firstTimestep;

      if (ppbonded)    ppbonded->getData(first, step, lattice, total);
      if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
      if (ppint)       ppint->getData(first, step, lattice, total);
      for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
      pressureProfileCount++;

      if (!(step % freq)) {
        // convert NAMD internal virial to pressure in units of bar 
        BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
   
        iout << "PRESSUREPROFILE: " << step << " ";
        if (step == first) {
          // output pressure profile for this step
          for (int i=0; i<arraysize; i++) {
            iout << total[i] * scalefac << " ";
          }
        } else {
          // output pressure profile averaged over the last count steps.
          scalefac /= pressureProfileCount;
          for (int i=0; i<arraysize; i++) 
            iout << pressureProfileAverage[i]*scalefac << " ";
        }
        iout << "\n" << endi; 
       
        // Clear the average for the next block
        memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
        pressureProfileCount = 0;
      }
      delete [] total;
    }
  
   if ( step != simParams->firstTimestep || stepInFullRun == 0 ) {  // skip repeated first step
    if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
     int benchPhase = stepInFullRun / simParams->firstLdbStep;
     if ( benchPhase > 0 && benchPhase < 10 ) {
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
      if ( simParams->computeEnergies < 60 ) {
        iout << iWARN << "Energy evaluation is expensive, increase computeEnergies to improve performance.\n";
        iout << iWARN << "If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n";
      }
#endif
      iout << iINFO;
      if ( benchPhase < 4 ) iout << "Initial time: ";
      else iout << "Benchmark time: ";
      iout << CkNumPes() << " CPUs ";

      {
        BigReal wallPerStep =
          (namdWallTimer() - startBenchTime) / simParams->firstLdbStep;
        BigReal ns = simParams->dt / 1000000.0;
        BigReal days = 1.0 / (24.0 * 60.0 * 60.0);

        if (simParams->nsPerDayOn) {
          BigReal nsPerDay = ns / (wallPerStep * days);
          iout << wallPerStep << " s/step " << nsPerDay << " ns/day ";
        }
        else {
          BigReal daysPerNano = wallPerStep * days / ns;
          iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
        }
        iout << memusage_MB() << " MB memory\n" << endi;
      }

     }
     startBenchTime = namdWallTimer();
    }
    ++stepInFullRun;
   }

    printTiming(step);

    Vector pairVDWForce, pairElectForce;
    if ( simParams->pairInteractionOn ){
      GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
      GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
    }

    // Compute cumulative nonequilibrium work (including shadow work).
    if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
      work = totalEnergy - totalEnergy0 - heat;
    }

    // callback to Tcl with whatever we can
#ifdef NAMD_TCL
#define CALLBACKDATA(LABEL,VALUE) \
                labels << (LABEL) << " "; values << (VALUE) << " ";
#define CALLBACKLIST(LABEL,VALUE) \
                labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
    if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {

      std::ostringstream labels, values;
#if CMK_BLUEGENEL
      // the normal version below gives a compiler error
      values.precision(16);
#else
      values << std::setprecision(16);
#endif
      CALLBACKDATA("TS",step);
      CALLBACKDATA("BOND",bondEnergy);
      CALLBACKDATA("ANGLE",angleEnergy);
      CALLBACKDATA("DIHED",dihedralEnergy);
      CALLBACKDATA("CROSS",crosstermEnergy);
      CALLBACKDATA("IMPRP",improperEnergy);
      CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
      CALLBACKDATA("VDW",ljEnergy+ljEnergySlow);
      CALLBACKDATA("BOUNDARY",boundaryEnergy);
      CALLBACKDATA("MISC",miscEnergy);
      CALLBACKDATA("POTENTIAL",potentialEnergy);
      CALLBACKDATA("KINETIC",kineticEnergy);
      CALLBACKDATA("TOTAL",totalEnergy);
      CALLBACKDATA("TEMP",temperature);
      CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
      CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
      CALLBACKDATA("VOLUME",lattice.volume());
      CALLBACKLIST("CELL_A",lattice.a());
      CALLBACKLIST("CELL_B",lattice.b());
      CALLBACKLIST("CELL_C",lattice.c());
      CALLBACKLIST("CELL_O",lattice.origin());
      labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
                << lattice.b_p() << " " << lattice.c_p() << "} ";
      if (simParams->drudeOn) {
        CALLBACKDATA("DRUDEBOND",drudeBondTemp);
      }
      if ( simParams->pairInteractionOn ) {
        CALLBACKLIST("VDW_FORCE",pairVDWForce);
        CALLBACKLIST("ELECT_FORCE",pairElectForce);
      }
      if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
        CALLBACKLIST("HEAT",heat);
        CALLBACKLIST("WORK",work);
      }
      if (simParams->alchOn) {
        if (simParams->alchThermIntOn) {
          CALLBACKLIST("BOND1", bondedEnergy_ti_1);
          CALLBACKLIST("ELEC1", electEnergy_ti_1 + electEnergySlow_ti_1 +
                                electEnergyPME_ti_1);
          CALLBACKLIST("VDW1", ljEnergy_ti_1);
          CALLBACKLIST("BOND2", bondedEnergy_ti_2);
          CALLBACKLIST("ELEC2", electEnergy_ti_2 + electEnergySlow_ti_2 +
                                electEnergyPME_ti_2);
          CALLBACKLIST("VDW2", ljEnergy_ti_2);
          if (simParams->alchLambdaFreq > 0) {
            CALLBACKLIST("CUMALCHWORK", cumAlchWork);
          }
        } else if (simParams->alchFepOn) {
          CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
                                improperEnergy + bondedEnergyDiff_f);
          CALLBACKLIST("ELEC2", electEnergy_f + electEnergySlow_f);
          CALLBACKLIST("VDW2", ljEnergy_f);
        } 
      }

      labels << '\0';  values << '\0';  // insane but makes Linux work
      state->callback_labelstring = labels.str();
      state->callback_valuestring = values.str();
      // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
    }
#undef CALLBACKDATA
#endif

    if ( ! cudaIntegrator) {
      drudeBondTempAvg += drudeBondTemp;
      temp_avg += temperature;
      pressure_avg += trace(pressure)/3.;
      groupPressure_avg += trace(groupPressure)/3.;
      avg_count += 1;
    }

    if ( simParams->outputPairlists && pairlistWarnings &&
                                ! (step % simParams->outputPairlists) ) {
      iout << iINFO << pairlistWarnings <<
        " pairlist warnings in past " << simParams->outputPairlists <<
        " steps.\n" << endi;
      pairlistWarnings = 0;
    }

    BigReal enthalpy;
    if (simParams->multigratorOn && ((step % simParams->computeEnergies) == 0)) {
      enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
    }

    if (simParams->IMDon) {
      IMDOutput *imd = NULL;
#ifdef NODEGROUP_FORCE_REGISTER
      if (simParams->CUDASOAintegrate) {
        CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
        imd = cpdata.ckLocalBranch()->imd;
      } else
#endif
      {
        imd = Node::Object()->imd;
      }
      if ( !(step % simParams->IMDfreq) &&
          ((simParams->IMDversion == IMDversion_t::IMDv3)
          && (simParams->IMDsendsettings.time_switch == 1)
          && (step != simParams->firstTimestep)) ) {

        IMDTime time = {simParams->dt * 0.001, simParams->dt * 0.001 * step, step}; // convert to ps

        if (imd != NULL) imd->gather_time(&time);
      }

      if (!(step % simParams->IMDfreq) && ((simParams->IMDversion == IMDversion_t::IMDv2) ||
          ((simParams->IMDversion == IMDversion_t::IMDv3) && (simParams->IMDsendsettings.energies_switch == 1) && (step != simParams->firstTimestep)))) {
        IMDEnergies energies;
        energies.tstep = step;
        energies.T = temp_avg/avg_count;
        energies.Etot = totalEnergy;
        energies.Epot = potentialEnergy;
        energies.Evdw = ljEnergy + ljEnergySlow;
        energies.Eelec = electEnergy + electEnergySlow;
        energies.Ebond = bondEnergy;
        energies.Eangle = angleEnergy;
        energies.Edihe = dihedralEnergy + crosstermEnergy;
        energies.Eimpr = improperEnergy;
        if (imd != NULL) imd->gather_energies(&energies);
      }

      if (!simParams->CUDASOAintegrate &&
          (simParams->IMDversion == IMDversion_t::IMDv3) &&
          !(step % simParams->IMDfreq) &&
          ((simParams->IMDsendsettings.time_switch == 1) || (simParams->IMDsendsettings.energies_switch == 1)) &&
          (step != simParams->firstTimestep)) {
        broadcast->IMDTimeEnergyBarrier.publish(step, 1);
      }
    }
    // XXX
    // Important note: in CPU-only/GPU-offload modes
    // printEnergies is called EVERY STEP.
    // Reduction averages are summed above and then later are output 
    // divided by count. Sums and counter are reset after printing.
    // XXX
    // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
    if ( ! minimize &&  step % simParams->outputEnergies ) return;

    if (cudaIntegrator) {
      totalEnergyAverage.addSample(totalEnergy);
      temperatureAverage.addSample(temperature);
      pressureAverage.addSample((1./3)*trace(pressure));
      groupPressureAverage.addSample((1./3)*trace(groupPressure));
    }

    // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!

    if ( marginViolations ) {
      iout << iERROR << marginViolations <<
        " margin violations detected since previous energy output.\n" << endi;
    }
    marginViolations = 0;

    int precision = simParams->outputEnergiesPrecision;
    if ( (step % (10 * (minimize?1:simParams->outputEnergies) ) ) == 0 )
    {
        iout << "ETITLE:      TS";
        iout << FORMAT("BOND", precision);
        iout << FORMAT("ANGLE", precision);
        iout << FORMAT("DIHED", precision);
        if ( ! simParams->mergeCrossterms ) iout << FORMAT("CROSS", precision);
        iout << FORMAT("IMPRP", precision);
        iout << "     ";
        iout << FORMAT("ELECT", precision);
        iout << FORMAT("VDW", precision);
        iout << FORMAT("BOUNDARY", precision);
        iout << FORMAT("MISC", precision);
        iout << FORMAT("KINETIC", precision);
        iout << "     ";
        iout << FORMAT("TOTAL", precision);
        iout << FORMAT("TEMP", precision);
        iout << FORMAT("POTENTIAL", precision);
        if (cudaIntegrator) {
          iout << FORMAT("TOTALAVG", precision);
        }
        else {
          // iout << FORMAT("TOTAL2", precision);
          iout << FORMAT("TOTAL3", precision);
        }
        iout << FORMAT("TEMPAVG", precision);
        if ( volume != 0. ) {
          iout << "     ";
          iout << FORMAT("PRESSURE", precision);
          iout << FORMAT("GPRESSURE", precision);
          iout << FORMAT("VOLUME", precision);
          iout << FORMAT("PRESSAVG", precision);
          iout << FORMAT("GPRESSAVG", precision);
        }
        if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
          iout << "     ";
          iout << FORMAT("HEAT", precision);
          iout << FORMAT("WORK", precision);
        }
        if (simParams->drudeOn) {
          iout << "     ";
          iout << FORMAT("DRUDEBOND", precision);
          iout << FORMAT("DRBONDAVG", precision);
        }
        // Ported by JLai
        if (simParams->goGroPair) {
          iout << "     ";
          iout << FORMAT("GRO_PAIR_LJ", precision);
          iout << FORMAT("GRO_PAIR_GAUSS", precision);
        }

        if (simParams->goForcesOn) {
          iout << "     ";
          iout << FORMAT("NATIVE", precision);
          iout << FORMAT("NONNATIVE", precision);
          //iout << FORMAT("REL_NATIVE", precision);
          //iout << FORMAT("REL_NONNATIVE", precision);
          iout << FORMAT("GOTOTAL", precision);
          //iout << FORMAT("GOAVG", precision);
        }
        // End of port -- JLai

        if (simParams->alchOn) {
          if (simParams->alchThermIntOn) {
            iout << "\nTITITLE:     TS";
            iout << FORMAT("BOND1", precision);
            iout << FORMAT("ELECT1", precision);
            iout << FORMAT("VDW1", precision);
            iout << FORMAT("BOND2", precision);
            iout << "     ";
            iout << FORMAT("ELECT2", precision);
            iout << FORMAT("VDW2", precision);
            if (simParams->alchLambdaFreq > 0) {
              iout << FORMAT("LAMBDA", precision);
              iout << FORMAT("ALCHWORK", precision);
              iout << FORMAT("CUMALCHWORK", precision);
            }
          } else if (simParams->alchFepOn) {
            iout << "\nFEPTITLE:    TS";
            iout << FORMAT("BOND2", precision);
            iout << FORMAT("ELECT2", precision);
            iout << FORMAT("VDW2", precision);
            if (simParams->alchLambdaFreq > 0) {
              iout << FORMAT("LAMBDA", precision);
            }
          }
        }

        iout << "\n\n" << endi;
        
        if (simParams->qmForcesOn) {
            iout << "QMETITLE:      TS";
            iout << FORMAT("QMID", precision);
            iout << FORMAT("ENERGY", precision);
            if (simParams->PMEOn) iout << FORMAT("PMECORRENERGY", precision);
            iout << "\n\n" << endi;
        }
        
    }

    // N.B.  HP's aCC compiler merges FORMAT calls in the same expression.
    //       Need separate statements because data returned in static array.
    iout << ETITLE(step);
    iout << FORMAT(bondEnergy, precision);
    iout << FORMAT(angleEnergy, precision);
    if ( simParams->mergeCrossterms ) {
      iout << FORMAT(dihedralEnergy+crosstermEnergy, precision);
    } else {
      iout << FORMAT(dihedralEnergy, precision);
      iout << FORMAT(crosstermEnergy, precision);
    }
    iout << FORMAT(improperEnergy, precision);
    iout << "     ";
    iout << FORMAT(electEnergy+electEnergySlow, precision);
    iout << FORMAT(ljEnergy+ljEnergySlow, precision);
    iout << FORMAT(boundaryEnergy, precision);
    iout << FORMAT(miscEnergy, precision);
    iout << FORMAT(kineticEnergy, precision);
    iout << "     ";
    iout << FORMAT(totalEnergy, precision);
    iout << FORMAT(temperature, precision);
    iout << FORMAT(potentialEnergy, precision);
    if (cudaIntegrator) {
      iout << FORMAT(totalEnergyAverage.average(), precision);
      iout << FORMAT(temperatureAverage.average(), precision);
    }
    else {
      // iout << FORMAT(flatEnergy, precision);
      iout << FORMAT(smoothEnergy, precision);
      iout << FORMAT(temp_avg/avg_count, precision);
    }
    if ( volume != 0. )
    {
        iout << "     ";
        iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3., precision);
        iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3., precision);
        iout << FORMAT(volume, precision);
        if (cudaIntegrator) {
          iout << FORMAT(pressureAverage.average()*PRESSUREFACTOR, precision);
          iout << FORMAT(groupPressureAverage.average()*PRESSUREFACTOR, precision);
        }
        else {
          iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count, precision);
          iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count, precision);
        }
    }
    if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
          iout << "     ";
          iout << FORMAT(heat, precision);
          iout << FORMAT(work, precision);
    }
    if (simParams->drudeOn) {
        iout << "     ";
        iout << FORMAT(drudeBondTemp, precision);
        iout << FORMAT(drudeBondTempAvg/avg_count, precision);
    }
    // Ported by JLai
    if (simParams->goGroPair) {
      iout << "     ";
      iout << FORMAT(groLJEnergy, precision);
      iout << FORMAT(groGaussEnergy, precision);
    }

    if (simParams->goForcesOn) {
      iout << "     ";
      iout << FORMAT(goNativeEnergy, precision);
      iout << FORMAT(goNonnativeEnergy, precision);
      //iout << FORMAT(relgoNativeEnergy, precision);
      //iout << FORMAT(relgoNonnativeEnergy, precision);
      iout << FORMAT(goTotalEnergy, precision);
      //iout << FORMAT("not implemented", precision);
    } // End of port -- JLai

    if (simParams->alchOn) { 
      if (simParams->alchThermIntOn) {
        iout << "\n";
        iout << TITITLE(step);
        iout << FORMAT(bondedEnergy_ti_1, precision);
        iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 + 
                       electEnergyPME_ti_1, precision);
        iout << FORMAT(ljEnergy_ti_1, precision);
        iout << FORMAT(bondedEnergy_ti_2, precision);
        iout << "     ";
        iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 +
                       electEnergyPME_ti_2, precision);
        iout << FORMAT(ljEnergy_ti_2, precision);
        if (simParams->alchLambdaFreq > 0) {
          iout << FORMAT(simParams->getCurrentLambda(step), precision);
          iout << FORMAT(alchWork, precision);
          iout << FORMAT(cumAlchWork, precision);
        }
      } else if (simParams->alchFepOn) {
        iout << "\n";
        iout << FEPTITLE2(step);
        iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy 
                       + improperEnergy + bondedEnergyDiff_f, precision);
        iout << FORMAT(electEnergy_f + electEnergySlow_f, precision);
        iout << FORMAT(ljEnergy_f, precision);
        if (simParams->alchLambdaFreq > 0) {
          iout << FORMAT(simParams->getCurrentLambda(step), precision);
        }
      }
    }

    iout << "\n\n" << endi;

#if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
     char webout[80];
     sprintf(webout,"%d %d %d %d",(int)totalEnergy,
             (int)(potentialEnergy),
             (int)kineticEnergy,(int)temperature);
     CApplicationDepositNode0Data(webout);
#endif

    if (simParams->pairInteractionOn) {
      iout << "PAIR INTERACTION:";
      iout << " STEP: " << step;
      iout << " VDW_FORCE: ";
      iout << FORMAT(pairVDWForce.x, precision);
      iout << FORMAT(pairVDWForce.y, precision);
      iout << FORMAT(pairVDWForce.z, precision);
      iout << " ELECT_FORCE: ";
      iout << FORMAT(pairElectForce.x, precision);
      iout << FORMAT(pairElectForce.y, precision);
      iout << FORMAT(pairElectForce.z, precision);
      iout << "\n" << endi;
    }
    drudeBondTempAvg = 0;
    temp_avg = 0;
    pressure_avg = 0;
    groupPressure_avg = 0;
    avg_count = 0;

    if ( fflush_count ) {
      --fflush_count;
      fflush(stdout);
    }
}

void Controller::writeExtendedSystemLabels(ofstream_namd &file) {
  Lattice &lattice = state->lattice;
  file << "#$LABELS step";
  if ( lattice.a_p() ) file << " a_x a_y a_z";
  if ( lattice.b_p() ) file << " b_x b_y b_z";
  if ( lattice.c_p() ) file << " c_x c_y c_z";
  file << " o_x o_y o_z";
  if ( simParams->langevinPistonOn ) {
    file << " s_x s_y s_z s_u s_v s_w";
  }
  if (simParams->monteCarloPressureOn) {
    file << " v_x v_y v_z";
  }
  file << std::endl;
}

void Controller::writeExtendedSystemData(int step, ofstream_namd &file) {
  Lattice &lattice = state->lattice;
  file.precision(12);
  file << step;
    if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
    if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
    if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
    file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
  if ( simParams->langevinPistonOn ) {
    Vector strainRate = diagonal(langevinPiston_strainRate);
    Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
    file << " " << strainRate.x;
    file << " " << strainRate.y;
    file << " " << strainRate.z;
    file << " " << strainRate2.x;
    file << " " << strainRate2.y;
    file << " " << strainRate2.z;
  }
  if (simParams->monteCarloPressureOn) {
    // We could print the elements as they are, but it would be meaningless,
    // because for constant ratio and isotropic fluctuations, we only operate
    // on x element. 
    if (!simParams->useFlexibleCell) {
      // isotropic fluctuations: use x element for v_x, v_y, and v_z
      file << " " << monteCarloMaxVolume.x;
      file << " " << monteCarloMaxVolume.x;
      file << " " << monteCarloMaxVolume.x;
    } else {
      // we have flexible cell
      if (simParams->useConstantRatio) {
        // use x element for v_x and v_y
        file << " " << monteCarloMaxVolume.x;
        file << " " << monteCarloMaxVolume.x;
        file << " " << monteCarloMaxVolume.z;
      } else {
        // constant area or just flexible in each axis dimension
        file << " " << monteCarloMaxVolume.x;
        file << " " << monteCarloMaxVolume.y;
        file << " " << monteCarloMaxVolume.z;
      }
    }
  }
  file << std::endl;
}

void Controller::enqueueCollections(int timestep)
{
  int prec;
  int dcdIndex;

  std::tie(prec, dcdIndex) = Output::coordinateNeeded(timestep);
#ifndef MEM_OPT_VERSION
  if ( prec == 4 ){
    collection->enqueuePositionsDcdSelection(timestep,state->lattice);
  }
  else
#endif
    if ( prec ){
    collection->enqueuePositions(timestep,state->lattice);
  }
  if ( Output::velocityNeeded(timestep) )
    collection->enqueueVelocities(timestep);
  if ( Output::forceNeeded(timestep) )
    collection->enqueueForces(timestep);
}

//Modifications for alchemical fep
void Controller::outputFepEnergy(int step) {
  if (simParams->alchOn && simParams->alchFepOn) {
    const int stepInRun = step - simParams->firstTimestep;
    const int alchEquilSteps = simParams->alchEquilSteps;
    const BigReal alchLambda = simParams->alchLambda;
    const bool alchEnsembleAvg = simParams->alchEnsembleAvg;

    if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
      FepNo = 0;
      exp_dE_ByRT = 0.0;
      net_dE = 0.0;
    }
    dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
         electEnergy - electEnergySlow - ljEnergy;
    BigReal RT = BOLTZMANN * simParams->alchTemp;

    if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
        (step / simParams->alchIDWSFreq) % 2 == 1 ) 
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
        &&
        // some people rely on the ensemble averages (dE_avg and dG) in fep output for calculating the free energy changes.
        // in the GPU version, dE is only available every outputEnergies steps.
        // so I check if the step is multiples of outputEnergies here.
        // I am still not sure whether we should check the stepInRun or step.
        // JH: dE is also calculated at alchOutFreq
        (step % simParams->computeEnergies == 0 || step % simParams->alchOutFreq == 0)
#endif
    ) {
      // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
      FepNo++;
      exp_dE_ByRT += exp(-dE/RT);
      net_dE += dE;
    }

    if (simParams->alchOutFreq) {
      if (stepInRun == 0) {
        if (!fepFile.is_open()) {
          NAMD_backup_file(simParams->alchOutFile);
          fepFile.open(simParams->alchOutFile);
          iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
          if(alchEnsembleAvg){
            fepSum = 0.0;
            fepFile << "#            STEP                 Elec                            "
                    << "vdW                    dE           dE_avg         Temp             dG\n"
                    << "#                           l             l+dl      "
                    << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
          }
          else{
            fepFile << "#            STEP                 Elec                            "
                    << "vdW                    dE         Temp\n"
                    << "#                           l             l+dl      "
                    << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
          }
        }
        if(!step){
          fepFile << "#NEW FEP WINDOW: "
                  << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
                  << simParams->alchLambda2;
          if ( simParams->alchLambdaIDWS >= 0. ) {
            fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
          }
          fepFile << std::endl;
        }
      }
      if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
        fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
                << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
                << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
      }
      if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
        writeFepEnergyData(step, fepFile);
        fepFile.flush();
      }
      if (alchEnsembleAvg && (step == simParams->N)) {
        fepSum = fepSum + dG;
        fepFile << "#Free energy change for lambda window [ " << alchLambda
                << " " << simParams->alchLambda2 << " ] is " << dG
                << " ; net change until now is " << fepSum << std::endl;
        fepFile.flush();
      }
    }
  }
}

void Controller::outputTiEnergy(int step) {
 if (simParams->alchOn && simParams->alchThermIntOn) {
  const int stepInRun = step - simParams->firstTimestep;
  const int alchEquilSteps = simParams->alchEquilSteps;
  const int alchLambdaFreq = simParams->alchLambdaFreq;

  if (stepInRun == 0 || stepInRun == alchEquilSteps) {
    TiNo = 0;
    net_dEdl_bond_1 = 0;
    net_dEdl_bond_2 = 0;
    net_dEdl_elec_1 = 0;
    net_dEdl_elec_2 = 0;
    net_dEdl_lj_1 = 0;
    net_dEdl_lj_2 = 0;
  }
  // The energies are only computed if they are required to output in CUDA build.
  // In the case when step % simParams->computeEnergies != 0 the energies are 0.
  // All output terms in writeTiEnergyData use TiNo and recent_TiNo as denominators,
  // so the TiNo and recent_TiNo should NOT be increased.
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
  if (step % simParams->computeEnergies == 0) {
#endif
    TiNo++;
    net_dEdl_bond_1 += bondedEnergy_ti_1;
    net_dEdl_bond_2 += bondedEnergy_ti_2;
    net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
                        + electEnergyPME_ti_1);
    net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
                        + electEnergyPME_ti_2);
    net_dEdl_lj_1 += ljEnergy_ti_1;
    net_dEdl_lj_2 += ljEnergy_ti_2;
    // Compute global dE / dLambda for lambda-dynamics
    computeTIderivative(step);
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
  }
#endif

  // Don't accumulate block averages (you'll get a divide by zero!) or even 
  // open the TI output if alchOutFreq is zero.
  if (simParams->alchOutFreq) {
    if (stepInRun == 0 || stepInRun == alchEquilSteps 
        || (! ((step - 1) % simParams->alchOutFreq))) {
      // output of instantaneous dU/dl now replaced with running average
      // over last alchOutFreq steps (except for step 0)
      recent_TiNo = 0;
      recent_dEdl_bond_1 = 0;
      recent_dEdl_bond_2 = 0;
      recent_dEdl_elec_1 = 0;
      recent_dEdl_elec_2 = 0;
      recent_dEdl_lj_1 = 0;
      recent_dEdl_lj_2 = 0;
      recent_alchWork = 0;
    }
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    if (step % simParams->computeEnergies == 0) {
#endif
      recent_TiNo++;
      recent_dEdl_bond_1 += bondedEnergy_ti_1;
      recent_dEdl_bond_2 += bondedEnergy_ti_2;
      recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 
                          + electEnergyPME_ti_1); 
      recent_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 
                          + electEnergyPME_ti_2); 
      recent_dEdl_lj_1 += ljEnergy_ti_1;
      recent_dEdl_lj_2 += ljEnergy_ti_2;
      recent_alchWork += alchWork;
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    }
#endif

    if (stepInRun == 0) {
      if (!tiFile.is_open()) {
        NAMD_backup_file(simParams->alchOutFile);
        tiFile.open(simParams->alchOutFile);
        /* BKR - This has been rather drastically updated to better match 
           stdout. This was necessary for several reasons:
           1) PME global scaling is obsolete (now removed)
           2) scaling of bonded terms was added
           3) alchemical work is now accumulated when switching is active
         */
        iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
        tiFile << "#TITITLE:    TS";
        tiFile << FORMAT("BOND1");
        tiFile << FORMAT("AVGBOND1");
        tiFile << FORMAT("ELECT1");
        tiFile << FORMAT("AVGELECT1");
        tiFile << "     ";
        tiFile << FORMAT("VDW1");
        tiFile << FORMAT("AVGVDW1");
        tiFile << FORMAT("BOND2");
        tiFile << FORMAT("AVGBOND2");
        tiFile << FORMAT("ELECT2");
        tiFile << "     ";
        tiFile << FORMAT("AVGELECT2");
        tiFile << FORMAT("VDW2");
        tiFile << FORMAT("AVGVDW2");
        if (alchLambdaFreq > 0) {
          tiFile << FORMAT("ALCHWORK");
          tiFile << FORMAT("CUMALCHWORK");
        }
        tiFile << std::endl;
      }

      if (alchLambdaFreq > 0) {
        tiFile << "#ALCHEMICAL SWITCHING ACTIVE " 
               << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
               << "\n#LAMBDA SCHEDULE: " 
               << "dL: " << simParams->getLambdaDelta() 
               << " Freq: " << alchLambdaFreq;
      }
      else {
        const BigReal alchLambda = simParams->alchLambda;    
        const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
        const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
        const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
        const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
        const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
        const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
        tiFile << "#NEW TI WINDOW: "
               << "LAMBDA " << alchLambda 
               << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
               << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
               << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
               << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
      }
      tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
             << std::endl;
    }

    if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
      tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
             << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
             << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
    }
    if ((step%simParams->alchOutFreq) == 0) {
      writeTiEnergyData(step, tiFile);
      tiFile.flush();
    }
  }
 }
}

/* 
 Work is accumulated whenever alchLambda changes.  In the current scheme we
 always increment lambda _first_, then integrate in time.  Therefore the work 
 is wrt the "old" lambda before the increment.
*/
BigReal Controller::computeAlchWork(const int step) {
  // alchemical scaling factors for groups 1/2 at the previous lambda
  const BigReal oldLambda = simParams->getCurrentLambda(step-1);
  const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
  const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
  const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
  const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
  const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
  // alchemical scaling factors for groups 1/2 at the new lambda
  const BigReal alchLambda = simParams->getCurrentLambda(step);
  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
  const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
  const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
  const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
  const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda);

  return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
          (elec_lambda_12 - elec_lambda_1)*
          (electEnergy_ti_1 + electEnergySlow_ti_1 +  electEnergyPME_ti_1) +
          (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
          (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
          (elec_lambda_22 - elec_lambda_2)*
          (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) + 
          (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
  );
}

void Controller::computeTIderivative(const int step) {
    // Note: WCA is incompatible with TI
    BigReal l1 = simParams->getCurrentLambda(step);
    BigReal l2 = 1. - l1;
  
    BigReal dl_bond_dl1 = (l1 >= simParams->alchBondLambdaEnd ?
                          0. : 1. / simParams->alchBondLambdaEnd);
    BigReal dl_bond_dl2 = (l2 >= simParams->alchBondLambdaEnd ?
                          0. : -1. / simParams->alchBondLambdaEnd);

    BigReal dl_lj_dl1   = (l1 >= simParams->alchVdwLambdaEnd ?
                          0. : 1. / simParams->alchVdwLambdaEnd);
    BigReal dl_lj_dl2   = (l2 >= simParams->alchVdwLambdaEnd ?
                          0. : -1. / simParams->alchVdwLambdaEnd);

    BigReal dl_elec_dl1 = (l1 <= simParams->alchElecLambdaStart ?
                          0. : 1. / (1. - simParams->alchElecLambdaStart));
    BigReal dl_elec_dl2 = (l2 <= simParams->alchElecLambdaStart ?
                          0. : -1. / (1. - simParams->alchElecLambdaStart));

    TIderivative =
            dl_bond_dl1 * bondedEnergy_ti_1
          + dl_bond_dl2 * bondedEnergy_ti_2
          + dl_lj_dl1   * ljEnergy_ti_1
          + dl_lj_dl2   * ljEnergy_ti_2
          + dl_elec_dl1 * (electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1)
          + dl_elec_dl2 * (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2);
}

void Controller::writeFepEnergyData(int step, ofstream_namd &file) {
  BigReal eeng = electEnergy + electEnergySlow;
  BigReal eeng_f = electEnergy_f + electEnergySlow_f;
  BigReal RT = BOLTZMANN * simParams->alchTemp;

  const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
  const int stepInRun = step - simParams->firstTimestep;
  const BigReal alchLambda = simParams->alchLambda;

  if(stepInRun){
    if ( simParams->alchLambdaIDWS >= 0. &&
        (step / simParams->alchIDWSFreq) % 2 == 0 ) {
      // IDWS is active and we are on a "backward" timestep
      fepFile << FEPTITLE_BACK(step);
    } else {
      // "forward" timestep
      fepFile << FEPTITLE(step);
    }
    fepFile << FORMAT(eeng);
    fepFile << FORMAT(eeng_f);
    fepFile << FORMAT(ljEnergy);
    fepFile << FORMAT(ljEnergy_f);
    fepFile << FORMAT(dE);
    if(alchEnsembleAvg){
      BigReal dE_avg = net_dE / FepNo;
      fepFile << FORMAT(dE_avg);
    }
    fepFile << FORMAT(temperature);
    if(alchEnsembleAvg){
      dG = -(RT * log(exp_dE_ByRT / FepNo));
      fepFile << FORMAT(dG);
    }
    fepFile << std::endl;
  }
}

void Controller::writeTiEnergyData(int step, ofstream_namd &file) {
  tiFile << TITITLE(step);
  tiFile << FORMAT(recent_dEdl_bond_1 / recent_TiNo);
  tiFile << FORMAT(net_dEdl_bond_1 / TiNo);
  tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
  tiFile << FORMAT(net_dEdl_elec_1/TiNo);
  tiFile << "     ";
  tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
  tiFile << FORMAT(net_dEdl_lj_1/TiNo);
  tiFile << FORMAT(recent_dEdl_bond_2 / recent_TiNo);
  tiFile << FORMAT(net_dEdl_bond_2 / TiNo);
  tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
  tiFile << "     ";
  tiFile << FORMAT(net_dEdl_elec_2/TiNo);
  tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
  tiFile << FORMAT(net_dEdl_lj_2/TiNo);
  if (simParams->alchLambdaFreq > 0) {
    tiFile << FORMAT(recent_alchWork / recent_TiNo);
    tiFile << FORMAT(cumAlchWork);
  }
  tiFile << std::endl;
}
//fepe

void Controller::outputExtendedSystem(int step)
{

  if ( step >= 0 ) {

    // Write out eXtended System Trajectory (XST) file
    if ( simParams->xstFrequency &&
         ((step % simParams->xstFrequency) == 0) )
    {
      if ( ! xstFile.is_open() )
      {
        iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
        NAMD_backup_file(simParams->xstFilename);
        xstFile.open(simParams->xstFilename);
        while (!xstFile) {
          if ( errno == EINTR ) {
            CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
            xstFile.clear();
            xstFile.open(simParams->xstFilename);
            continue;
          }
          char err_msg[257];
          sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
          NAMD_err(err_msg);
        }
        xstFile << "# NAMD extended system trajectory file" << std::endl;
        writeExtendedSystemLabels(xstFile);
      }
      writeExtendedSystemData(step,xstFile);
      xstFile.flush();
    }

    // Write out eXtended System Configuration (XSC) files
    //  Output a restart file
    if ( simParams->restartFrequency &&
         ((step % simParams->restartFrequency) == 0) &&
         (step != simParams->firstTimestep) )
    {
      iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
                << step << "\n" << endi;
      char fname[NAMD_FILENAME_BUFFER_SIZE];
      const char *bsuffix = ".old";
      strcpy(fname, simParams->restartFilename);
      if ( simParams->restartSave ) {
        char timestepstr[20];
        sprintf(timestepstr,".%d",step);
        strcat(fname, timestepstr);
        bsuffix = ".BAK";
      }
      strcat(fname, ".xsc");
      NAMD_backup_file(fname,bsuffix);
      ofstream_namd xscFile(fname);
      while (!xscFile) {
        if ( errno == EINTR ) {
          CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
          xscFile.clear();
          xscFile.open(fname);
          continue;
        }
        char err_msg[257];
        sprintf(err_msg, "Error opening XSC restart file %s",fname);
        NAMD_err(err_msg);
      } 
      xscFile << "# NAMD extended system configuration restart file" << std::endl;
      writeExtendedSystemLabels(xscFile);
      writeExtendedSystemData(step,xscFile);
      if (!xscFile) {
        char err_msg[257];
        sprintf(err_msg, "Error writing XSC restart file %s",fname);
        NAMD_err(err_msg);
      } 
    }

  }

  //  Output final coordinates
  if (step == FILE_OUTPUT || step == END_OF_RUN)
  {
    int realstep = ( step == FILE_OUTPUT ?
        simParams->firstTimestep : simParams->N );
    iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
                << realstep << "\n" << endi;
    static char fname[NAMD_FILENAME_BUFFER_SIZE];
    strcpy(fname, simParams->outputFilename);
    strcat(fname, ".xsc");
    NAMD_backup_file(fname);
    ofstream_namd xscFile(fname);
    while (!xscFile) {
      if ( errno == EINTR ) {
        CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
        xscFile.clear();
        xscFile.open(fname);
        continue;
      }
      char err_msg[257];
      sprintf(err_msg, "Error opening XSC output file %s",fname);
      NAMD_err(err_msg);
    } 
    xscFile << "# NAMD extended system configuration output file" << std::endl;
    writeExtendedSystemLabels(xscFile);
    writeExtendedSystemData(realstep,xscFile);
    if (!xscFile) {
      char err_msg[257];
      sprintf(err_msg, "Error writing XSC output file %s",fname);
      NAMD_err(err_msg);
    } 
  }

  //  Close trajectory file
  if (step == END_OF_RUN) {
    if ( xstFile.is_open() ) {
      xstFile.close();
      iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
    }
  }

}


void Controller::recvCheckpointReq(const char *key, int task, checkpoint &cp) {  // responding replica
  switch ( task ) {
    case SCRIPT_CHECKPOINT_STORE:
      if ( ! checkpoints.count(key) ) {
        checkpoints[key] = new checkpoint;
      }
      *checkpoints[key] = cp;
      break;
    case SCRIPT_CHECKPOINT_LOAD:
      if ( ! checkpoints.count(key) ) {
        NAMD_die("Unable to load checkpoint, requested key was never stored.");
      }
      cp = *checkpoints[key];
      break;
    case SCRIPT_CHECKPOINT_SWAP:
      if ( ! checkpoints.count(key) ) {
        NAMD_die("Unable to swap checkpoint, requested key was never stored.");
      }
      std::swap(cp,*checkpoints[key]);
      break;
    case SCRIPT_CHECKPOINT_FREE:
      if ( ! checkpoints.count(key) ) {
        NAMD_die("Unable to free checkpoint, requested key was never stored.");
      }
      delete checkpoints[key];
      checkpoints.erase(key);
      break;
  }
}

void Controller::recvCheckpointAck(checkpoint &cp) {  // initiating replica
  if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
    state->lattice = cp.lattice;
    *(ControllerState*)this = cp.state;
  }
  CkpvAccess(_qd)->process();
}


void Controller::rebalanceLoad(int step)
{
  if ( ! ldbSteps ) { 
    ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
  }
  if ( ! --ldbSteps ) {
    startBenchTime -= namdWallTimer();
        Node::Object()->outputPatchComputeMaps("before_ldb", step);
    LdbCoordinator::Object()->rebalance(this);  
        startBenchTime += namdWallTimer();
    fflush_count = 3;
  }
}

void Controller::cycleBarrier(int doBarrier, int step) {
#if USE_BARRIER
        if (doBarrier) {
          broadcast->cycleBarrier.publish(step,1);
          CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
                  namdWallTimer()-firstWTime,CmiTimer()-firstCTime);
        }
#endif
}

void Controller::traceBarrier(int turnOnTrace, int step) {
        CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, namdWallTimer()-firstWTime,CmiTimer()-firstCTime);
        CProxy_Node nd(CkpvAccess(BOCclass_group).node);
        nd.traceBarrier(turnOnTrace, step);
        CthSuspend();
}

void Controller::resumeAfterTraceBarrier(int step){
        broadcast->traceBarrier.publish(step,1);
        CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, namdWallTimer()-firstWTime,CmiTimer()-firstCTime);
        awaken();
}

#ifdef MEASURE_NAMD_WITH_PAPI
void Controller::papiMeasureBarrier(int turnOnMeasure, int step){
        CkPrintf("Cycle time at PAPI measurement sync (begin) Wall at step %d: %f CPU %f\n", step, namdWallTimer()-firstWTime,CmiTimer()-firstCTime);
        CProxy_Node nd(CkpvAccess(BOCclass_group).node);
        nd.papiMeasureBarrier(turnOnMeasure, step);
        CthSuspend();
}

void Controller::resumeAfterPapiMeasureBarrier(int step){
        broadcast->papiMeasureBarrier.publish(step,1);
        CkPrintf("Cycle time at PAPI measurement sync (end) Wall at step %d: %f CPU %f\n", step, namdWallTimer()-firstWTime,CmiTimer()-firstCTime);
        awaken();
}
#endif

void Controller::terminate(void) {
  BackEnd::awaken();
  CthFree(threadWrapper->thread);
  delete threadWrapper;
  CthSuspend();
}

RequireReduction* Controller::getCurrentReduction() {
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
  return (simParams->CUDASOAintegrate) ? reductionGpuResident : 
                                         reductionBasic;
#else
  return reductionBasic;
#endif
}