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

Generated on Sat Sep 23 01:17:13 2017 for NAMD by  doxygen 1.4.7