Controller.C

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

Generated on Tue Oct 16 01:17:14 2018 for NAMD by  doxygen 1.4.7