Controller.C

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

Generated on Fri May 25 01:17:12 2018 for NAMD by  doxygen 1.4.7