NAMD: Controller.C Source File

00001 
00007 /*****************************************************************************
00008  * $Source: /home/cvs/namd/cvsroot/namd2/src/Controller.C,v $
00009  * $Author: jim $
00010  * $Date: 2009/11/12 01:54:06 $
00011  * $Revision: 1.1241 $
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 <iomanip>
00040 
00041 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
00042 extern "C" void CApplicationDepositNode0Data(char *);
00043 #endif
00044 
00045 #ifndef cbrt
00046   // cbrt() not in math.h on goneril
00047   #define cbrt(x)  pow(x,(double)(1.0/3.0))
00048 #endif
00049 
00050 //#define DEBUG_PRESSURE
00051 #define MIN_DEBUG_LEVEL 3
00052 //#define DEBUGM
00053 #include "Debug.h"
00054 
00055 #define XXXBIGREAL 1.0e32
00056 
00057 class PressureProfileReduction {
00058 private:
00059   RequireReduction *reduction;
00060   const int nslabs;
00061   const int freq;
00062   int nelements;
00063   int count;
00064   char *name;
00065 
00066   BigReal *average;
00067 
00068 public:
00069   PressureProfileReduction(int rtag, int numslabs, int numpartitions,
00070       const char *myname, int outputfreq)
00071   : nslabs(numslabs), freq(outputfreq) {
00072     name = strdup(myname);
00073     nelements = 3*nslabs * numpartitions;
00074     reduction = ReductionMgr::Object()->willRequire(rtag,nelements);
00075 
00076     average = new BigReal[nelements];
00077     count = 0;
00078   }
00079   ~PressureProfileReduction() {
00080     delete [] average;
00081     delete reduction;
00082     free(name);
00083   }
00084   // 
00085   void getData(int firsttimestep, int step, const Lattice &lattice, 
00086       BigReal *total) {
00087     reduction->require();
00088 
00089     int i;
00090     double inv_volume = 1.0 / lattice.volume();
00091     // accumulate the pressure profile computed for this step into the average.
00092     int arraysize = 3*nslabs;
00093     for (i=0; i<nelements; i++) {
00094       BigReal val = reduction->item(i) * inv_volume;
00095       average[i] += val;
00096       total[i % arraysize] += val;
00097     }
00098     count++;
00099     if (!(step % freq)) {
00100       // convert NAMD internal virial to pressure in units of bar 
00101       BigReal scalefac = PRESSUREFACTOR * nslabs;
00102  
00103       iout << "PPROFILE" << name << ": " << step << " ";
00104       if (step == firsttimestep) {
00105         // output pressure profile for this step
00106         for (i=0; i<nelements; i++) 
00107           iout << reduction->item(i)*scalefac*inv_volume << " ";
00108       } else {
00109         // output pressure profile averaged over the last Freq steps.
00110         scalefac /= count;
00111         for (i=0; i<nelements; i++) 
00112           iout << average[i]*scalefac << " ";
00113       }
00114       iout << "\n" << endi; 
00115       // Clear the average for the next block
00116       memset(average, 0, nelements*sizeof(BigReal));
00117       count = 0;
00118     }
00119   }
00120 };
00121 
00122 
00123 Controller::Controller(NamdState *s) :
00124         computeChecksum(0), marginViolations(0), pairlistWarnings(0),
00125         simParams(Node::Object()->simParameters),
00126         state(s),
00127 #ifdef MEM_OPT_VERSION
00128         collection(CollectionMasterHandler::Object()),  
00129 #else
00130         collection(CollectionMaster::Object()),
00131 #endif
00132         startCTime(0),
00133         firstCTime(CmiTimer()),
00134         startWTime(0),
00135         firstWTime(CmiWallTimer()),
00136         startBenchTime(0),
00137         ldbSteps(0),
00138         fflush_count(3)
00139 {
00140     broadcast = new ControllerBroadcasts;
00141     reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
00142     // pressure profile reductions
00143     pressureProfileSlabs = 0;
00144     pressureProfileCount = 0;
00145     ppbonded = ppnonbonded = ppint = NULL;
00146     if (simParams->pressureProfileOn || simParams->pressureProfileEwaldOn) {
00147       int ntypes = simParams->pressureProfileAtomTypes;
00148       int nslabs = pressureProfileSlabs = simParams->pressureProfileSlabs;
00149       // number of partitions for pairwise interactions
00150       int npairs = (ntypes * (ntypes+1))/2;
00151       pressureProfileAverage = new BigReal[3*nslabs];
00152       memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
00153       int freq = simParams->pressureProfileFreq;
00154       if (simParams->pressureProfileOn) {
00155         ppbonded = new PressureProfileReduction(REDUCTIONS_PPROF_BONDED, 
00156             nslabs, npairs, "BONDED", freq);
00157         ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
00158             nslabs, npairs, "NONBONDED", freq);
00159         // internal partitions by atom type, but not in a pairwise fashion
00160         ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL, 
00161             nslabs, ntypes, "INTERNAL", freq);
00162       } else {
00163         // just doing Ewald, so only need nonbonded
00164         ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
00165             nslabs, npairs, "NONBONDED", freq);
00166       }
00167     }
00168     random = new Random(simParams->randomSeed);
00169     random->split(0,PatchMap::Object()->numPatches()+1);
00170 
00171     rescaleVelocities_sumTemps = 0;  rescaleVelocities_numTemps = 0;
00172     berendsenPressure_avg = 0; berendsenPressure_count = 0;
00173     // strainRate tensor is symmetric to avoid rotation
00174     langevinPiston_strainRate =
00175         Tensor::symmetric(simParams->strainRate,simParams->strainRate2);
00176     if ( ! simParams->useFlexibleCell ) {
00177       BigReal avg = trace(langevinPiston_strainRate) / 3.;
00178       langevinPiston_strainRate = Tensor::identity(avg);
00179     } else if ( simParams->useConstantRatio ) {
00180 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
00181                  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
00182       AVGXY(langevinPiston_strainRate);
00183 #undef AVGXY
00184     }
00185     smooth2_avg = XXXBIGREAL;
00186     temp_avg = 0;
00187     pressure_avg = 0;
00188     groupPressure_avg = 0;
00189     avg_count = 0;
00190     pressure_tavg = 0;
00191     groupPressure_tavg = 0;
00192     tavg_count = 0;
00193     checkpoint_stored = 0;
00194     drudeComTemp = 0;
00195     drudeBondTemp = 0;
00196     drudeComTempAvg = 0;
00197     drudeBondTempAvg = 0;
00198 }
00199 
00200 Controller::~Controller(void)
00201 {
00202     delete broadcast;
00203     delete reduction;
00204     delete ppbonded;
00205     delete ppnonbonded;
00206     delete ppint;
00207     delete [] pressureProfileAverage;
00208     delete random;
00209 }
00210 
00211 void Controller::threadRun(Controller* arg)
00212 {
00213     arg->algorithm();
00214 }
00215 
00216 void Controller::run(void)
00217 {
00218     // create a Thread and invoke it
00219     DebugM(4, "Starting thread in controller on this=" << this << "\n");
00220     thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
00221     CthSetStrategyDefault(thread);
00222 #if CMK_BLUEGENE_CHARM
00223     BgAttach(thread);
00224 #endif
00225     awaken();
00226 }
00227 
00228 
00229 void Controller::algorithm(void)
00230 {
00231   int scriptTask;
00232   int scriptSeq = 0;
00233   BackEnd::awaken();
00234   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
00235     switch ( scriptTask ) {
00236       case SCRIPT_OUTPUT:
00237         enqueueCollections(FILE_OUTPUT);
00238         outputExtendedSystem(FILE_OUTPUT);
00239         break;
00240       case SCRIPT_MEASURE:
00241         enqueueCollections(EVAL_MEASURE);
00242         break;
00243       case SCRIPT_REINITVELS:
00244         iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
00245           << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
00246         break;
00247       case SCRIPT_RESCALEVELS:
00248         iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
00249           << " BY " << simParams->scriptArg1 << "\n" << endi;
00250         break;
00251       case SCRIPT_CHECKPOINT:
00252         iout << "CHECKPOINTING POSITIONS AT STEP " << simParams->firstTimestep
00253           << "\n" << endi;
00254         checkpoint_stored = 1;
00255         checkpoint_lattice = state->lattice;
00256         checkpoint_langevinPiston_strainRate = langevinPiston_strainRate;
00257         checkpoint_berendsenPressure_avg = berendsenPressure_avg;
00258         checkpoint_berendsenPressure_count = berendsenPressure_count;
00259         break;
00260       case SCRIPT_REVERT:
00261         iout << "REVERTING POSITIONS AT STEP " << simParams->firstTimestep
00262           << "\n" << endi;
00263         if ( ! checkpoint_stored )
00264           NAMD_die("Unable to revert, checkpoint was never called!");
00265         state->lattice = checkpoint_lattice;
00266         langevinPiston_strainRate = checkpoint_langevinPiston_strainRate;
00267         berendsenPressure_avg = checkpoint_berendsenPressure_avg;
00268         berendsenPressure_count = checkpoint_berendsenPressure_count;
00269         break;
00270       case SCRIPT_MINIMIZE:
00271         minimize();
00272         break;
00273       case SCRIPT_RUN:
00274         integrate();
00275         break;
00276     }
00277     BackEnd::awaken();
00278   }
00279   enqueueCollections(END_OF_RUN);
00280   outputExtendedSystem(END_OF_RUN);
00281   // note: this is a Converse interface call that gets translated to a
00282   // null method call if CMK_OPTIMIZE is set.
00283   if(traceAvailable())
00284     traceClose();
00285   terminate();
00286 }
00287 
00288 
00289 //
00290 // XXX static and global variables are unsafe for shared memory builds.
00291 // The use of global and static vars should be eliminated.
00292 //
00293 extern int eventEndOfTimeStep;
00294 
00295 // Handle SIGINT so that restart files get written completely.
00296 static int gotsigint = 0;
00297 static void my_sigint_handler(int sig) {
00298   if (sig == SIGINT) gotsigint = 1;
00299 }
00300 extern "C" {
00301   typedef void (*namd_sighandler_t)(int);
00302 }
00303 
00304 void Controller::integrate() {
00305 
00306     int step = simParams->firstTimestep;
00307 
00308     const int numberOfSteps = simParams->N;
00309     const int stepsPerCycle = simParams->stepsPerCycle;
00310 
00311     const int zeroMomentum = simParams->zeroMomentum;
00312 
00313     nbondFreq = simParams->nonbondedFrequency;
00314     const int dofull = ( simParams->fullDirectOn ||
00315                         simParams->FMAOn || simParams->PMEOn );
00316     if (dofull)
00317       slowFreq = simParams->fullElectFrequency;
00318     else
00319       slowFreq = simParams->nonbondedFrequency;
00320 
00321     reassignVelocities(step);  // only for full-step velecities
00322     receivePressure(step);
00323     if ( zeroMomentum && dofull && ! (step % slowFreq) )
00324                                                 correctMomentum(step);
00325     printFepMessage(step);
00326     printTiMessage(step);
00327     printDynamicsEnergies(step);
00328     outputFepEnergy(step);
00329     outputTiEnergy(step);
00330     traceUserEvent(eventEndOfTimeStep);
00331     outputExtendedSystem(step);
00332     rebalanceLoad(step);
00333 
00334     // Handling SIGINT doesn't seem to be working on Lemieux, and it
00335     // sometimes causes the net-xxx versions of NAMD to segfault on exit, 
00336     // so disable it for now.
00337     // namd_sighandler_t oldhandler = signal(SIGINT, 
00338     //  (namd_sighandler_t)my_sigint_handler);
00339     for ( ++step ; step <= numberOfSteps; ++step )
00340     {
00341 
00342         rescaleVelocities(step);
00343         tcoupleVelocities(step);
00344         berendsenPressure(step);
00345         langevinPiston1(step);
00346         enqueueCollections(step);  // after lattice scaling!
00347         receivePressure(step);
00348         if ( zeroMomentum && dofull && ! (step % slowFreq) )
00349                                                 correctMomentum(step);
00350         langevinPiston2(step);
00351         reassignVelocities(step);
00352         printDynamicsEnergies(step);
00353         outputFepEnergy(step);
00354         outputTiEnergy(step);
00355         traceUserEvent(eventEndOfTimeStep);
00356   // if (gotsigint) {
00357   //   iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
00358   //   NAMD_quit();
00359   // }
00360         outputExtendedSystem(step);
00361 #if CYCLE_BARRIER
00362         cycleBarrier(!((step+1) % stepsPerCycle),step);
00363 #elif  PME_BARRIER
00364         cycleBarrier(dofull && !(step%slowFreq),step);
00365 #elif  STEP_BARRIER
00366         cycleBarrier(1, step);
00367 #endif
00368 
00369         rebalanceLoad(step);
00370 
00371 #if  PME_BARRIER
00372         cycleBarrier(dofull && !((step+1)%slowFreq),step);   // step before PME
00373 #endif
00374     }
00375     // signal(SIGINT, oldhandler);
00376 }
00377 
00378 
00379 #define CALCULATE \
00380   printMinimizeEnergies(step); \
00381   outputExtendedSystem(step); \
00382   rebalanceLoad(step); \
00383   if ( step == numberOfSteps ) return; \
00384   else ++step;
00385 
00386 #define MOVETO(X) \
00387   if ( step == numberOfSteps ) { \
00388     if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \
00389     if ( newDir || (mid.x-last.x) ) { \
00390       broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \
00391     } else { \
00392       broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00393       broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00394       broadcast->minimizeCoefficient.publish(minSeq++,0.); \
00395     } \
00396     enqueueCollections(step); \
00397     CALCULATE \
00398   } else if ( (X)-last.x ) { \
00399     broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \
00400     newDir = 0; \
00401     last.x = (X); \
00402     enqueueCollections(step); \
00403     CALCULATE \
00404     last.u = min_energy; \
00405     last.dudx = -1. * min_f_dot_v; \
00406     last.noGradient = min_huge_count; \
00407     if ( minVerbose ) { \
00408       iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \
00409       if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \
00410       iout << "\n" << endi; \
00411     } \
00412   }
00413 
00414 struct minpoint {
00415   BigReal x, u, dudx, noGradient;
00416   minpoint() : x(0), u(0), dudx(0), noGradient(1) { ; }
00417 };
00418 
00419 void Controller::minimize() {
00420   // iout << "Controller::minimize() called.\n" << endi;
00421 
00422   const int minVerbose = simParams->minVerbose;
00423   const int numberOfSteps = simParams->N;
00424   int step = simParams->firstTimestep;
00425   slowFreq = nbondFreq = 1;
00426   BigReal tinystep = simParams->minTinyStep;  // 1.0e-6
00427   BigReal babystep = simParams->minBabyStep;  // 1.0e-2
00428   BigReal linegoal = simParams->minLineGoal;  // 1.0e-3
00429   BigReal initstep = tinystep;
00430   const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
00431 
00432   CALCULATE
00433 
00434   int minSeq = 0;
00435 
00436   // just move downhill until initial bad contacts go away or 100 steps
00437   int old_min_huge_count = -1;
00438   int steps_at_same_min_huge_count = 0;
00439   for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
00440     if ( min_huge_count == old_min_huge_count ) {
00441       if ( ++steps_at_same_min_huge_count > 10 ) break;
00442     } else {
00443       old_min_huge_count = min_huge_count;
00444       steps_at_same_min_huge_count = 0;
00445     }
00446     iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
00447     broadcast->minimizeCoefficient.publish(minSeq++,1.);
00448     enqueueCollections(step);
00449     CALCULATE
00450   }
00451   if ( min_huge_count ) {
00452     iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
00453   }
00454   iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00455 
00456   int atStart = 2;
00457   int errorFactor = 10;
00458   BigReal old_f_dot_f = min_f_dot_f;
00459   broadcast->minimizeCoefficient.publish(minSeq++,0.);
00460   broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
00461   int newDir = 1;
00462   min_f_dot_v = min_f_dot_f;
00463   min_v_dot_v = min_f_dot_f;
00464   while ( 1 ) {
00465     // line minimization
00466     // bracket minimum on line
00467     minpoint lo,hi,mid,last;
00468     BigReal x = 0;
00469     lo.x = x;
00470     lo.u = min_energy;
00471     lo.dudx = -1. * min_f_dot_v;
00472     lo.noGradient = min_huge_count;
00473     mid = lo;
00474     last = mid;
00475     if ( minVerbose ) {
00476       iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
00477       if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
00478       iout << "\n" << endi;
00479     }
00480     BigReal tol = fabs( linegoal * min_f_dot_v );
00481     if ( initstep > babystep ) initstep = babystep;
00482     if ( initstep < 1.0e-300 ) initstep = 1.0e-300;
00483     iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
00484             fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
00485     x = initstep;
00486     x *= sqrt( min_f_dot_f / min_v_dot_v ); MOVETO(x)
00487     while ( min_huge_count ) {
00488       x *= 0.25;  MOVETO(x);
00489       initstep *= 0.25;
00490     }
00491     // bracket minimum on line
00492     initstep *= 0.25;
00493     BigReal maxinitstep = initstep * 16.0;
00494     while ( last.u < mid.u ) {
00495       initstep *= 2.0;
00496       lo = mid; mid = last;
00497       x *= 2.0; MOVETO(x)
00498     }
00499     if ( initstep > maxinitstep ) initstep = maxinitstep;
00500     hi = last;
00501 #define PRINT_BRACKET \
00502     iout << "LINE MINIMIZER BRACKET: DX " \
00503          << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
00504         " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
00505         lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
00506     PRINT_BRACKET
00507     // converge on minimum on line
00508     int itcnt;
00509     for ( itcnt = 10; itcnt > 0; --itcnt ) {
00510       int progress = 1;
00511       // select new position
00512       if ( mid.noGradient ) {
00513        if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) {  // subdivide left side
00514         x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
00515         MOVETO(x)
00516         if ( last.u <= mid.u ) {
00517           hi = mid; mid = last;
00518         } else {
00519           lo = last;
00520         }
00521        } else {  // subdivide right side
00522         x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
00523         MOVETO(x)
00524         if ( last.u <= mid.u ) {
00525           lo = mid; mid = last;
00526         } else {
00527           hi = last;
00528         }
00529        }
00530       } else {
00531        if ( mid.dudx > 0. ) {  // subdivide left side
00532         BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
00533         BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
00534         x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
00535         BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
00536         if ( xdiv ) x /= xdiv; else x = last.x;
00537         if ( x > altxhi ) x = altxhi;
00538         if ( x < altxlo ) x = altxlo;
00539         if ( x-last.x == 0 ) break;
00540         MOVETO(x)
00541         if ( last.u <= mid.u ) {
00542           hi = mid; mid = last;
00543         } else {
00544           if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
00545           lo = last;
00546         }
00547        } else {  // subdivide right side
00548         BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
00549         BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
00550         x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
00551         BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
00552         if ( xdiv ) x /= xdiv; else x = last.x;
00553         if ( x < altxlo ) x = altxlo;
00554         if ( x > altxhi ) x = altxhi;
00555         if ( x-last.x == 0 ) break;
00556         MOVETO(x)
00557         if ( last.u <= mid.u ) {
00558           lo = mid; mid = last;
00559         } else {
00560           if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
00561           hi = last;
00562         }
00563        }
00564       }
00565       PRINT_BRACKET
00566       if ( ! progress ) {
00567         MOVETO(mid.x);
00568         break;
00569       }
00570       if ( fabs(last.dudx) < tol ) break;
00571       if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
00572       if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
00573     }
00574     // new direction
00575     broadcast->minimizeCoefficient.publish(minSeq++,0.);
00576     BigReal c = min_f_dot_f / old_f_dot_f;
00577     c = ( c > 1.5 ? 1.5 : c );
00578     if ( atStart ) { c = 0; --atStart; }
00579     if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
00580       c = 0;
00581       if ( errorFactor < 100 ) errorFactor += 10;
00582     }
00583     if ( c == 0 ) {
00584       iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00585     }
00586     broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
00587     newDir = 1;
00588     old_f_dot_f = min_f_dot_f;
00589     min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
00590     min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
00591   }
00592 
00593 }
00594 
00595 #undef MOVETO
00596 #undef CALCULATE
00597 
00598 void Controller::berendsenPressure(int step)
00599 {
00600   if ( simParams->berendsenPressureOn ) {
00601    berendsenPressure_count += 1;
00602    berendsenPressure_avg += controlPressure;
00603    const int freq = simParams->berendsenPressureFreq;
00604    if ( ! (berendsenPressure_count % freq) ) {
00605     Tensor factor = berendsenPressure_avg / berendsenPressure_count;
00606     berendsenPressure_avg = 0;
00607     berendsenPressure_count = 0;
00608     // We only use on-diagonal terms (for now)
00609     factor = Tensor::diagonal(diagonal(factor));
00610     factor -= Tensor::identity(simParams->berendsenPressureTarget);
00611     factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
00612        simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
00613     factor += Tensor::identity(1.0);
00614 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
00615          if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
00616          if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
00617     int limited = 0;
00618     LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
00619     LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
00620     LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
00621 #undef LIMIT_SCALING
00622     if ( limited ) {
00623       iout << iERROR << "Step " << step <<
00624         " cell rescaling factor limited.\n" << endi;
00625     }
00626     broadcast->positionRescaleFactor.publish(step,factor);
00627     state->lattice.rescale(factor);
00628    }
00629   } else {
00630     berendsenPressure_avg = 0;
00631     berendsenPressure_count = 0;
00632   }
00633 }
00634 
00635 void Controller::langevinPiston1(int step)
00636 {
00637   if ( simParams->langevinPistonOn )
00638   {
00639     Tensor &strainRate = langevinPiston_strainRate;
00640     int cellDims = simParams->useFlexibleCell ? 1 : 3;
00641     BigReal dt = simParams->dt;
00642     BigReal dt_long = slowFreq * dt;
00643     BigReal kT = BOLTZMAN * simParams->langevinPistonTemp;
00644     BigReal tau = simParams->langevinPistonPeriod;
00645     BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00646 
00647 #ifdef DEBUG_PRESSURE
00648     iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
00649 #endif
00650 
00651     if ( ! ( (step-1) % slowFreq ) )
00652     {
00653       BigReal gamma = 1 / simParams->langevinPistonDecay;
00654       BigReal f1 = exp( -0.5 * dt_long * gamma );
00655       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00656       strainRate *= f1;
00657       if ( simParams->useFlexibleCell ) {
00658         // We only use on-diagonal terms (for now)
00659         if ( simParams->useConstantRatio ) {
00660           BigReal r = f2 * random->gaussian();
00661           strainRate.xx += r;
00662           strainRate.yy += r;
00663           strainRate.zz += f2 * random->gaussian();
00664         } else {
00665           strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00666         }
00667       } else {
00668         strainRate += f2 * Tensor::identity(random->gaussian());
00669       }
00670 
00671 #ifdef DEBUG_PRESSURE
00672       iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00673 #endif
00674     }
00675 
00676     // Apply surface tension.  If surfaceTensionTarget is zero, we get
00677     // the default (isotropic pressure) case.
00678     
00679     Tensor ptarget;
00680     ptarget.zz = simParams->langevinPistonTarget;
00681     ptarget.xx = ptarget.yy = simParams->langevinPistonTarget - 
00682         simParams->surfaceTensionTarget / state->lattice.c().z;
00683 
00684     strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00685       ( controlPressure - ptarget );
00686 
00687     if (simParams->fixCellDims) {
00688       if (simParams->fixCellDimX) strainRate.xx = 0;
00689       if (simParams->fixCellDimY) strainRate.yy = 0;
00690       if (simParams->fixCellDimZ) strainRate.zz = 0;
00691     }
00692 
00693 
00694 #ifdef DEBUG_PRESSURE
00695     iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00696 #endif
00697 
00698     if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
00699     {
00700       // We only use on-diagonal terms (for now)
00701       Tensor factor;
00702       if ( !simParams->useConstantArea ) {
00703         factor.xx = exp( dt_long * strainRate.xx );
00704         factor.yy = exp( dt_long * strainRate.yy );
00705       } else {
00706         factor.xx = factor.yy = 1;
00707       }
00708       factor.zz = exp( dt_long * strainRate.zz );
00709       broadcast->positionRescaleFactor.publish(step,factor);
00710       state->lattice.rescale(factor);
00711 #ifdef DEBUG_PRESSURE
00712       iout << iINFO << "rescaling by: " << factor << "\n";
00713 #endif
00714     }
00715 
00716     // corrections to integrator
00717     if ( ! ( step % nbondFreq ) )
00718     {
00719 #ifdef DEBUG_PRESSURE
00720       iout << iINFO << "correcting strain rate for nbond, ";
00721 #endif
00722       strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00723                 ( (nbondFreq - 1) * controlPressure_nbond );
00724 #ifdef DEBUG_PRESSURE
00725       iout << "strain rate: " << strainRate << "\n";
00726 #endif
00727     }
00728     if ( ! ( step % slowFreq ) )
00729     {
00730 #ifdef DEBUG_PRESSURE
00731       iout << iINFO << "correcting strain rate for slow, ";
00732 #endif
00733       strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00734                 ( (slowFreq - 1) * controlPressure_slow );
00735 #ifdef DEBUG_PRESSURE
00736       iout << "strain rate: " << strainRate << "\n";
00737 #endif
00738     }
00739     if (simParams->fixCellDims) {
00740       if (simParams->fixCellDimX) strainRate.xx = 0;
00741       if (simParams->fixCellDimY) strainRate.yy = 0;
00742       if (simParams->fixCellDimZ) strainRate.zz = 0;
00743     }
00744 
00745   }
00746 
00747 }
00748 
00749 void Controller::langevinPiston2(int step)
00750 {
00751   if ( simParams->langevinPistonOn )
00752   {
00753     Tensor &strainRate = langevinPiston_strainRate;
00754     int cellDims = simParams->useFlexibleCell ? 1 : 3;
00755     BigReal dt = simParams->dt;
00756     BigReal dt_long = slowFreq * dt;
00757     BigReal kT = BOLTZMAN * simParams->langevinPistonTemp;
00758     BigReal tau = simParams->langevinPistonPeriod;
00759     BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00760 
00761     // corrections to integrator
00762     if ( ! ( step % nbondFreq ) )
00763     {
00764 #ifdef DEBUG_PRESSURE
00765       iout << iINFO << "correcting strain rate for nbond, ";
00766 #endif
00767       strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00768                 ( (nbondFreq - 1) * controlPressure_nbond );
00769 #ifdef DEBUG_PRESSURE
00770       iout << "strain rate: " << strainRate << "\n";
00771 #endif
00772     }
00773     if ( ! ( step % slowFreq ) )
00774     {
00775 #ifdef DEBUG_PRESSURE
00776       iout << iINFO << "correcting strain rate for slow, ";
00777 #endif
00778       strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00779                 ( (slowFreq - 1) * controlPressure_slow );
00780 #ifdef DEBUG_PRESSURE
00781       iout << "strain rate: " << strainRate << "\n";
00782 #endif
00783     }
00784 
00785     // Apply surface tension.  If surfaceTensionTarget is zero, we get
00786     // the default (isotropic pressure) case.
00787    
00788     Tensor ptarget;
00789     ptarget.zz = simParams->langevinPistonTarget;
00790     ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
00791         simParams->surfaceTensionTarget / state->lattice.c().z;
00792 
00793     strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00794       ( controlPressure - ptarget );
00795 
00796  
00797 #ifdef DEBUG_PRESSURE
00798     iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00799 #endif
00800 
00801     if ( ! ( step % slowFreq ) )
00802     {
00803       BigReal gamma = 1 / simParams->langevinPistonDecay;
00804       BigReal f1 = exp( -0.5 * dt_long * gamma );
00805       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00806       strainRate *= f1;
00807       if ( simParams->useFlexibleCell ) {
00808         // We only use on-diagonal terms (for now)
00809         if ( simParams->useConstantRatio ) {
00810           BigReal r = f2 * random->gaussian();
00811           strainRate.xx += r;
00812           strainRate.yy += r;
00813           strainRate.zz += f2 * random->gaussian();
00814         } else {
00815           strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00816         }
00817       } else {
00818         strainRate += f2 * Tensor::identity(random->gaussian());
00819       }
00820 #ifdef DEBUG_PRESSURE
00821       iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00822 #endif
00823     }
00824 
00825 #ifdef DEBUG_PRESSURE
00826     iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
00827 #endif
00828     if (simParams->fixCellDims) {
00829       if (simParams->fixCellDimX) strainRate.xx = 0;
00830       if (simParams->fixCellDimY) strainRate.yy = 0;
00831       if (simParams->fixCellDimZ) strainRate.zz = 0;
00832     }
00833   }
00834 
00835 
00836 }
00837 
00838 void Controller::rescaleVelocities(int step)
00839 {
00840   const int rescaleFreq = simParams->rescaleFreq;
00841   if ( rescaleFreq > 0 ) {
00842     rescaleVelocities_sumTemps += temperature;  ++rescaleVelocities_numTemps;
00843     if ( rescaleVelocities_numTemps == rescaleFreq ) {
00844       BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
00845       const BigReal rescaleTemp = simParams->rescaleTemp;
00846       BigReal factor = sqrt(rescaleTemp/avgTemp);
00847       broadcast->velocityRescaleFactor.publish(step,factor);
00848       iout << "RESCALING VELOCITIES AT STEP " << step
00849            << " FROM AVERAGE TEMPERATURE OF " << avgTemp
00850            << " TO " << rescaleTemp << " KELVIN.\n" << endi;
00851       rescaleVelocities_sumTemps = 0;  rescaleVelocities_numTemps = 0;
00852     }
00853   }
00854 }
00855 
00856 void Controller::correctMomentum(int step) {
00857 
00858     Vector momentum;
00859     momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
00860     momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
00861     momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
00862     const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
00863 
00864     if ( momentum.length2() == 0. )
00865       iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
00866     if ( mass == 0. )
00867       NAMD_die("Total mass is zero in Controller::correctMomentum");
00868 
00869     momentum *= (-1./mass);
00870 
00871     broadcast->momentumCorrection.publish(step+slowFreq,momentum);
00872 }
00873 
00874 //Modifications for alchemical fep
00875 void Controller::printFepMessage(int step)
00876 {
00877   if (simParams->alchFepOn) {
00878     const BigReal alchLambda = simParams->alchLambda;
00879     const BigReal alchLambda2 = simParams->alchLambda2;
00880     const BigReal alchTemp = simParams->alchTemp;
00881     const int alchEquilSteps = simParams->alchEquilSteps;
00882     iout << "FEP: RESETTING FOR NEW FEP WINDOW "
00883          << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2
00884          << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
00885          << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
00886          << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp 
00887          << " K FOR FEP CALCULATION\n" << endi;
00888   }
00889 } 
00890 void Controller::printTiMessage(int step)
00891 {
00892   if (simParams->alchThermIntOn) {
00893     const BigReal alchLambda = simParams->alchLambda;
00894     const BigReal alchTemp = simParams->alchTemp;
00895     const int alchEquilSteps = simParams->alchEquilSteps;
00896     iout << "TI: RESETTING FOR NEW WINDOW "
00897          << "LAMBDA SET TO " << alchLambda 
00898          << "\nTI: WINDOW TO HAVE " << alchEquilSteps 
00899          << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
00900   }
00901 } 
00902 //fepe
00903 
00904 void Controller::reassignVelocities(int step)
00905 {
00906   const int reassignFreq = simParams->reassignFreq;
00907   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
00908     BigReal newTemp = simParams->reassignTemp;
00909     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
00910     if ( simParams->reassignIncr > 0.0 ) {
00911       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
00912         newTemp = simParams->reassignHold;
00913     } else {
00914       if ( newTemp < simParams->reassignHold )
00915         newTemp = simParams->reassignHold;
00916     }
00917     iout << "REASSIGNING VELOCITIES AT STEP " << step
00918          << " TO " << newTemp << " KELVIN.\n" << endi;
00919   }
00920 }
00921 
00922 void Controller::tcoupleVelocities(int step)
00923 {
00924   if ( simParams->tCoupleOn )
00925   {
00926     const BigReal tCoupleTemp = simParams->tCoupleTemp;
00927     BigReal coefficient = 1.;
00928     if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
00929     broadcast->tcoupleCoefficient.publish(step,coefficient);
00930   }
00931 }
00932 
00933 static char *FORMAT(BigReal X)
00934 {
00935   static char tmp_string[25];
00936   const double maxnum = 9999999999.9999;
00937   if ( X > maxnum ) X = maxnum;
00938   if ( X < -maxnum ) X = -maxnum;
00939   sprintf(tmp_string," %14.4f",X); 
00940   return tmp_string;
00941 }
00942 
00943 static char *FORMAT(const char *X)
00944 {
00945   static char tmp_string[25];
00946   sprintf(tmp_string," %14s",X); 
00947   return tmp_string;
00948 }
00949 
00950 static char *ETITLE(int X)
00951 {
00952   static char tmp_string[21];
00953   sprintf(tmp_string,"ENERGY: %7d",X); 
00954   return  tmp_string;
00955 }
00956 
00957 void Controller::receivePressure(int step, int minimize)
00958 {
00959     Node *node = Node::Object();
00960     Molecule *molecule = node->molecule;
00961     SimParameters *simParameters = node->simParameters;
00962     Lattice &lattice = state->lattice;
00963 
00964     reduction->require();
00965 
00966     Tensor virial;
00967     Tensor virial_normal;
00968     Tensor virial_nbond;
00969     Tensor virial_slow;
00970 #ifdef ALTVIRIAL
00971     Tensor altVirial_normal;
00972     Tensor altVirial_nbond;
00973     Tensor altVirial_slow;
00974 #endif
00975     Tensor intVirial;
00976     Tensor intVirial_normal;
00977     Tensor intVirial_nbond;
00978     Tensor intVirial_slow;
00979     Vector extForce_normal;
00980     Vector extForce_nbond;
00981     Vector extForce_slow;
00982     BigReal volume;
00983 
00984 #if 1
00985     numDegFreedom = molecule->num_deg_freedom();
00986     int numGroupDegFreedom = molecule->num_group_deg_freedom();
00987     int numFixedGroups = molecule->num_fixed_groups();
00988     int numFixedAtoms = molecule->num_fixed_atoms();
00989 #endif
00990 #if 0
00991     int numAtoms = molecule->numAtoms;
00992     numDegFreedom = 3 * numAtoms;
00993     int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
00994     int numFixedAtoms =
00995         ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
00996     int numLonepairs = molecule->numLonepairs;
00997     int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
00998     if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
00999     if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01000     if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01001     if ( ! ( numFixedAtoms || molecule->numConstraints
01002         || simParameters->comMove || simParameters->langevinOn ) ) {
01003       numDegFreedom -= 3;
01004       numGroupDegFreedom -= 3;
01005     }
01006     if (simParameters->pairInteractionOn) {
01007       // this doesn't attempt to deal with fixed atoms or constraints
01008       numDegFreedom = 3 * molecule->numFepInitial;
01009     }
01010     int numRigidBonds = molecule->numRigidBonds;
01011     int numFixedRigidBonds =
01012         ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01013 
01014     // numLonepairs is subtracted here because all lonepairs have a rigid bond
01015     // to oxygen, but all of the LP degrees of freedom are dealt with above
01016     numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01017 #endif
01018 
01019     kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01020     kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01021 
01022     BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01023         reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01024     BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01025         reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01026 
01027     BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01028                                         / ( numDegFreedom * BOLTZMAN );
01029     BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01030                                         / ( numDegFreedom * BOLTZMAN );
01031     BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01032                                         / ( numGroupDegFreedom * BOLTZMAN );
01033     BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01034                                         / ( numGroupDegFreedom * BOLTZMAN );
01035 
01036     /*  test code for comparing different temperatures  
01037     iout << "TEMPTEST: " << step << " " << 
01038         atomTempHalfstep << " " <<
01039         atomTempCentered << " " <<
01040         groupTempHalfstep << " " <<
01041         groupTempCentered << "\n" << endi;
01042   iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
01043     numGroupDegFreedom << "\n" << endi;
01044      */
01045 
01046     GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01047     GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01048     GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01049 
01050 #ifdef ALTVIRIAL
01051     GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01052     GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01053     GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01054 #endif
01055 
01056     GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01057     GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01058     GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01059 
01060     GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01061     GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01062     GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01063     Vector extPosition = lattice.origin();
01064     virial_normal -= outer(extForce_normal,extPosition);
01065     virial_nbond -= outer(extForce_nbond,extPosition);
01066     virial_slow -= outer(extForce_slow,extPosition);
01067 
01068 
01069     kineticEnergy = kineticEnergyCentered;
01070     temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMAN );
01071 
01072     if (simParameters->drudeOn) {
01073       BigReal drudeComKE
01074         = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01075       BigReal drudeBondKE
01076         = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01077       int g_bond = 3 * molecule->numDrudeAtoms;
01078       int g_com = numDegFreedom - g_bond;
01079       drudeComTemp = 2.0 * drudeComKE / (g_com * BOLTZMAN);
01080       drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMAN)) : 0.);
01081     }
01082 
01083     if ( (volume=lattice.volume()) != 0. )
01084     {
01085 
01086       if (simParameters->LJcorrection && volume) {
01087         // Apply tail correction to pressure
01088         //printf("Volume is %f\n", volume);
01089         //printf("Applying tail correction of %f to virial\n", molecule->tail_corr_virial / volume);
01090         virial_normal += Tensor::identity(molecule->tail_corr_virial / volume);
01091       }
01092 
01093       // kinetic energy component included in virials
01094       pressure_normal = virial_normal / volume;
01095       groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01096 
01097       if ( minimize || ! ( step % nbondFreq ) )
01098       {
01099         pressure_nbond = virial_nbond / volume;
01100         groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01101       }
01102 
01103       if ( minimize || ! ( step % slowFreq ) )
01104       {
01105         pressure_slow = virial_slow / volume;
01106         groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01107       }
01108 
01109 /*
01110       iout << "VIRIALS: " << virial_normal << " " << virial_nbond << " " <<
01111         virial_slow << " " << ( virial_normal - intVirial_normal ) << " " <<
01112         ( virial_nbond - intVirial_nbond ) << " " <<
01113         ( virial_slow - intVirial_slow ) << "\n";
01114 */
01115 
01116       pressure = pressure_normal + pressure_nbond + pressure_slow; 
01117       groupPressure = groupPressure_normal + groupPressure_nbond +
01118                                                 groupPressure_slow;
01119     }
01120     else
01121     {
01122       pressure = Tensor();
01123       groupPressure = Tensor();
01124     }
01125 
01126     if ( simParameters->useGroupPressure )
01127     {
01128       controlPressure_normal = groupPressure_normal;
01129       controlPressure_nbond = groupPressure_nbond;
01130       controlPressure_slow = groupPressure_slow;
01131       controlPressure = groupPressure;
01132       controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01133       if ( ! ( numFixedAtoms || molecule->numConstraints
01134         || simParameters->comMove || simParameters->langevinOn ) ) {
01135         controlNumDegFreedom -= 1;
01136       }
01137     }
01138     else
01139     {
01140       controlPressure_normal = pressure_normal;
01141       controlPressure_nbond = pressure_nbond;
01142       controlPressure_slow = pressure_slow;
01143       controlPressure = pressure;
01144       controlNumDegFreedom = numDegFreedom / 3;
01145     }
01146 
01147     if (simParameters->fixCellDims) {
01148       if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01149       if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01150       if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01151     }
01152 
01153     if ( simParameters->useFlexibleCell ) {
01154       // use symmetric pressure to control rotation
01155       // controlPressure_normal = symmetric(controlPressure_normal);
01156       // controlPressure_nbond = symmetric(controlPressure_nbond);
01157       // controlPressure_slow = symmetric(controlPressure_slow);
01158       // controlPressure = symmetric(controlPressure);
01159       // only use on-diagonal components for now
01160       controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01161       controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01162       controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01163       controlPressure = Tensor::diagonal(diagonal(controlPressure));
01164       if ( simParameters->useConstantRatio ) {
01165 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01166                  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01167         AVGXY(controlPressure_normal);
01168         AVGXY(controlPressure_nbond);
01169         AVGXY(controlPressure_slow);
01170         AVGXY(controlPressure);
01171 #undef AVGXY
01172       }
01173     } else {
01174       controlPressure_normal =
01175                 Tensor::identity(trace(controlPressure_normal)/3.);
01176       controlPressure_nbond =
01177                 Tensor::identity(trace(controlPressure_nbond)/3.);
01178       controlPressure_slow =
01179                 Tensor::identity(trace(controlPressure_slow)/3.);
01180       controlPressure =
01181                 Tensor::identity(trace(controlPressure)/3.);
01182     }
01183 
01184 #ifdef DEBUG_PRESSURE
01185     iout << iINFO << "Control pressure = " << controlPressure <<
01186       " = " << controlPressure_normal << " + " <<
01187       controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01188 #endif
01189 
01190 }
01191 
01192 void Controller::compareChecksums(int step, int forgiving) {
01193     Node *node = Node::Object();
01194     Molecule *molecule = node->molecule;
01195 
01196     // Some basic correctness checking
01197     BigReal checksum, checksum_b;
01198 
01199     checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
01200     if ( ((int)checksum) != molecule->numAtoms )
01201       NAMD_bug("Bad global atom count!\n");
01202 
01203     checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
01204     if ( ((int)checksum) != computeChecksum ) {
01205       if ( computeChecksum )
01206         NAMD_bug("Bad global compute count!\n");
01207       else
01208         computeChecksum = ((int)checksum);
01209     }
01210 
01211     checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
01212     if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
01213       if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
01214         iout << iWARN << "Bad global bond count!\n" << endi;
01215       else NAMD_bug("Bad global bond count!\n");
01216     }
01217 
01218     checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
01219     if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
01220       if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
01221         iout << iWARN << "Bad global angle count!\n" << endi;
01222       else NAMD_bug("Bad global angle count!\n");
01223     }
01224 
01225     checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
01226     if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
01227       if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
01228         iout << iWARN << "Bad global dihedral count!\n" << endi;
01229       else NAMD_bug("Bad global dihedral count!\n");
01230     }
01231 
01232     checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
01233     if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
01234       if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
01235         iout << iWARN << "Bad global improper count!\n" << endi;
01236       else NAMD_bug("Bad global improper count!\n");
01237     }
01238 
01239     checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
01240     if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
01241       if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
01242         iout << iWARN << "Bad global crossterm count!\n" << endi;
01243       else NAMD_bug("Bad global crossterm count!\n");
01244     }
01245 
01246     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01247     if ( ((int)checksum) > molecule->numCalcExclusions &&
01248          ( ! simParams->mollyOn || step % slowFreq ) )
01249       iout << iWARN << "Not all atoms have unique coordinates.\n" << endi;
01250 
01251     checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
01252     if ( ((int)checksum) && ! marginViolations ) {
01253       iout << iERROR << "Margin is too small for " << ((int)checksum) <<
01254         " atoms during timestep " << step << ".\n" << iERROR <<
01255       "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
01256     }
01257     marginViolations += (int)checksum;
01258 
01259     checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
01260     if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
01261       iout << iINFO <<
01262         "Pairlistdist is too small for " << ((int)checksum) <<
01263         " computes during timestep " << step << ".\n" << endi;
01264     }
01265     if ( simParams->outputPairlists )  pairlistWarnings += (int)checksum;
01266 
01267     checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
01268     if ( checksum ) {
01269       if ( forgiving )
01270         iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
01271       else NAMD_bug("Stray PME grid charges detected!\n");
01272     }
01273 }
01274 
01275 void Controller::printTiming(int step) {
01276 
01277     if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
01278     {
01279       const double endWTime = CmiWallTimer() - firstWTime;
01280       const double endCTime = CmiTimer() - firstCTime;
01281 
01282       const double elapsedW = 
01283         (endWTime - startWTime) / simParams->outputTiming;
01284       const double elapsedC = 
01285         (endCTime - startCTime) / simParams->outputTiming;
01286 
01287       const double remainingW = elapsedW * (simParams->N - step);
01288       const double remainingW_hours = remainingW / 3600;
01289 
01290       startWTime = endWTime;
01291       startCTime = endCTime;
01292 
01293 #ifdef NAMD_CUDA
01294       if ( simParams->outputEnergies < 60 &&
01295            step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
01296         iout << iWARN << "Energy evaluation is done on CPU, increase outputEnergies to improve performance.\n" << endi;
01297       }
01298 #endif
01299       if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
01300         CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
01301                   ", %g hours remaining, %f MB of memory in use.\n",
01302                   step, endCTime, elapsedC, endWTime, elapsedW,
01303                   remainingW_hours, memusage_MB());
01304       }
01305     }
01306 }
01307 
01308 void Controller::printMinimizeEnergies(int step) {
01309 
01310     receivePressure(step,1);
01311 
01312     Node *node = Node::Object();
01313     Molecule *molecule = node->molecule;
01314     BigReal checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01315     if ( ((int)checksum) && ((int)checksum) < molecule->numCalcExclusions )
01316       iout << iWARN << "Bad global exclusion count, possible error!\n" << iWARN
01317         << "Increasing cutoff during minimization may avoid this.\n" << endi;
01318     compareChecksums(step,1);
01319 
01320     printEnergies(step,1);
01321 
01322     min_energy = totalEnergy;
01323     min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
01324     min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
01325     min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
01326     min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
01327 
01328 }
01329 
01330 void Controller::printDynamicsEnergies(int step) {
01331 
01332     Node *node = Node::Object();
01333     Molecule *molecule = node->molecule;
01334     SimParameters *simParameters = node->simParameters;
01335     Lattice &lattice = state->lattice;
01336 
01337     BigReal checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01338     if ( ((int)checksum) && ((int)checksum) < molecule->numCalcExclusions )
01339       NAMD_bug("Bad global exclusion count!\n");
01340     compareChecksums(step);
01341 
01342     printEnergies(step,0);
01343 }
01344 
01345 void Controller::printEnergies(int step, int minimize)
01346 {
01347     Node *node = Node::Object();
01348     Molecule *molecule = node->molecule;
01349     SimParameters *simParameters = node->simParameters;
01350     Lattice &lattice = state->lattice;
01351 
01352     BigReal bondEnergy;
01353     BigReal angleEnergy;
01354     BigReal dihedralEnergy;
01355     BigReal improperEnergy;
01356     BigReal crosstermEnergy;
01357     BigReal boundaryEnergy;
01358     BigReal miscEnergy;
01359     BigReal potentialEnergy;
01360     BigReal flatEnergy;
01361     BigReal smoothEnergy;
01362     Vector momentum;
01363     Vector angularMomentum;
01364     BigReal volume = lattice.volume();
01365 
01366     bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
01367     angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
01368     dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
01369     improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
01370     crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
01371     boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
01372     miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
01373 
01374     if ( minimize || ! ( step % nbondFreq ) )
01375     {
01376       electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
01377       ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
01378 //fepb
01379       electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
01380       ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
01381 
01382       electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
01383       ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
01384       electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
01385       ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
01386 //fepe
01387     }
01388 
01389     if ( minimize || ! ( step % slowFreq ) )
01390     {
01391       electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
01392 //fepb
01393       electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
01394 
01395       electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
01396       electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
01397       electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
01398       electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
01399 //fepe
01400     }
01401 
01402     if (simParameters->LJcorrection && volume) {
01403       // Apply tail correction to energy
01404       //printf("Volume is %f\n", volume);
01405       //printf("Applying tail correction of %f to energy\n", molecule->tail_corr_ener / volume);
01406       ljEnergy += molecule->tail_corr_ener / volume;
01407     }
01408 
01409 
01410     momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
01411     momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
01412     momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
01413     angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
01414     angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
01415     angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
01416 
01417     potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
01418         improperEnergy + electEnergy + electEnergySlow + ljEnergy +
01419         crosstermEnergy + boundaryEnergy + miscEnergy;
01420     totalEnergy = potentialEnergy + kineticEnergy;
01421     flatEnergy = totalEnergy +
01422         (1.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
01423     if ( !(step%slowFreq) ) {
01424       // only adjust based on most accurate energies
01425       BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
01426       if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
01427       smooth2_avg *= 0.9375;
01428       smooth2_avg += 0.0625 * s;
01429     }
01430     smoothEnergy = flatEnergy + smooth2_avg -
01431         (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
01432 
01433     if ( simParameters->outputMomenta && ! minimize &&
01434          ! ( step % simParameters->outputMomenta ) )
01435     {
01436       iout << "MOMENTUM: " << step 
01437            << " P: " << momentum
01438            << " L: " << angularMomentum
01439            << "\n" << endi;
01440     }
01441 
01442     if ( simParameters->outputPressure ) {
01443       pressure_tavg += pressure;
01444       groupPressure_tavg += groupPressure;
01445       tavg_count += 1;
01446       if ( minimize || ! ( step % simParameters->outputPressure ) ) {
01447 #ifdef NAMD_CUDA
01448         if ( 
01449           step < (simParams->firstTimestep + 10*simParams->outputPressure) ) {
01450           iout << iWARN << "Pressure tensor is incorrect.  CUDA currently evaluates scalar pressure only.\n" << endi;
01451         }
01452 #endif
01453         iout << "PRESSURE: " << step << " "
01454            << PRESSUREFACTOR * pressure << "\n"
01455            << "GPRESSURE: " << step << " "
01456            << PRESSUREFACTOR * groupPressure << "\n";
01457         if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
01458            << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
01459            << "GPRESSAVG: " << step << " "
01460            << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n" << endi;
01461         pressure_tavg = 0;
01462         groupPressure_tavg = 0;
01463         tavg_count = 0;
01464       }
01465     }
01466 
01467     // pressure profile reductions
01468     if (pressureProfileSlabs) {
01469       const int freq = simParams->pressureProfileFreq;
01470       const int arraysize = 3*pressureProfileSlabs;
01471       
01472       BigReal *total = new BigReal[arraysize];
01473       memset(total, 0, arraysize*sizeof(BigReal));
01474       const int first = simParams->firstTimestep;
01475 
01476       if (ppbonded)    ppbonded->getData(first, step, lattice, total);
01477       if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
01478       if (ppint)       ppint->getData(first, step, lattice, total);
01479       for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
01480       pressureProfileCount++;
01481 
01482       if (!(step % freq)) {
01483         // convert NAMD internal virial to pressure in units of bar 
01484         BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
01485    
01486         iout << "PRESSUREPROFILE: " << step << " ";
01487         if (step == first) {
01488           // output pressure profile for this step
01489           for (int i=0; i<arraysize; i++) {
01490             iout << total[i] * scalefac << " ";
01491           }
01492         } else {
01493           // output pressure profile averaged over the last count steps.
01494           scalefac /= pressureProfileCount;
01495           for (int i=0; i<arraysize; i++) 
01496             iout << pressureProfileAverage[i]*scalefac << " ";
01497         }
01498         iout << "\n" << endi; 
01499        
01500         // Clear the average for the next block
01501         memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
01502         pressureProfileCount = 0;
01503       }
01504       delete [] total;
01505     }
01506   
01507     int stepInRun = step - simParams->firstTimestep;
01508     if ( stepInRun % simParams->firstLdbStep == 0 ) {
01509      int benchPhase = stepInRun / simParams->firstLdbStep;
01510      if ( benchPhase > 0 && benchPhase < 7 ) {
01511 #ifdef NAMD_CUDA
01512       if ( simParams->outputEnergies < 60 ) {
01513         iout << iWARN << "Energy evaluation is done on CPU, increase outputEnergies to improve performance.\n";
01514       }
01515 #endif
01516       iout << iINFO;
01517       if ( benchPhase < 4 ) iout << "Initial time: ";
01518       else iout << "Benchmark time: ";
01519       iout << CkNumPes() << " CPUs ";
01520       {
01521         BigReal wallTime = CmiWallTimer() - startBenchTime;
01522         BigReal wallPerStep =
01523                 (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
01524         BigReal ns = simParams->dt / 1000000.0;
01525         BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
01526         BigReal daysPerNano = wallPerStep * days / ns;
01527         iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
01528         iout << memusage_MB() << " MB memory\n" << endi;
01529       }
01530      }
01531      startBenchTime = CmiWallTimer();
01532     }
01533 
01534     printTiming(step);
01535 
01536     Vector pairVDWForce, pairElectForce;
01537     if ( simParameters->pairInteractionOn ) {
01538       GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
01539       GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
01540     }
01541 
01542     // callback to Tcl with whatever we can
01543 #ifdef NAMD_TCL
01544 #define CALLBACKDATA(LABEL,VALUE) \
01545                 labels << (LABEL) << " "; values << (VALUE) << " ";
01546 #define CALLBACKLIST(LABEL,VALUE) \
01547                 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
01548     if (node->getScript() && node->getScript()->doCallback()) {
01549       std::ostringstream labels, values;
01550 #if CMK_BLUEGENEL
01551       // the normal version below gives a compiler error
01552       values.precision(16);
01553 #else
01554       values << std::setprecision(16);
01555 #endif
01556       CALLBACKDATA("TS",step);
01557       CALLBACKDATA("BOND",bondEnergy);
01558       CALLBACKDATA("ANGLE",angleEnergy);
01559       CALLBACKDATA("DIHED",dihedralEnergy);
01560       CALLBACKDATA("CROSS",crosstermEnergy);
01561       CALLBACKDATA("IMPRP",improperEnergy);
01562       CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
01563       CALLBACKDATA("VDW",ljEnergy);
01564       CALLBACKDATA("BOUNDARY",boundaryEnergy);
01565       CALLBACKDATA("MISC",miscEnergy);
01566       CALLBACKDATA("POTENTIAL",potentialEnergy);
01567       CALLBACKDATA("KINETIC",kineticEnergy);
01568       CALLBACKDATA("TOTAL",totalEnergy);
01569       CALLBACKDATA("TEMP",temperature);
01570       CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
01571       CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
01572       CALLBACKDATA("VOLUME",lattice.volume());
01573       CALLBACKLIST("CELL_A",lattice.a());
01574       CALLBACKLIST("CELL_B",lattice.b());
01575       CALLBACKLIST("CELL_C",lattice.c());
01576       CALLBACKLIST("CELL_O",lattice.origin());
01577       labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
01578                 << lattice.b_p() << " " << lattice.c_p() << "} ";
01579       if ( simParameters->pairInteractionOn ) {
01580         CALLBACKLIST("VDW_FORCE",pairVDWForce);
01581         CALLBACKLIST("ELECT_FORCE",pairElectForce);
01582       }
01583 
01584       labels << '\0';  values << '\0';  // insane but makes Linux work
01585       std::string labelstring = labels.str();
01586       std::string valuestring = values.str();
01587       node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
01588     }
01589 #undef CALLBACKDATA
01590 #endif
01591 
01592     drudeComTempAvg += drudeComTemp;
01593     drudeBondTempAvg += drudeBondTemp;
01594 
01595     temp_avg += temperature;
01596     pressure_avg += trace(pressure)/3.;
01597     groupPressure_avg += trace(groupPressure)/3.;
01598     avg_count += 1;
01599 
01600     if ( simParams->outputPairlists && pairlistWarnings &&
01601                                 ! (step % simParams->outputPairlists) ) {
01602       iout << iINFO << pairlistWarnings <<
01603         " pairlist warnings in past " << simParams->outputPairlists <<
01604         " steps.\n" << endi;
01605       pairlistWarnings = 0;
01606     }
01607     
01608     // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
01609     if ( ! minimize &&  step % simParameters->outputEnergies ) return;
01610     // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
01611 
01612     if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
01613       IMDEnergies energies;
01614       energies.tstep = step;
01615       energies.T = temp_avg/avg_count;
01616       energies.Etot = totalEnergy;
01617       energies.Epot = potentialEnergy;
01618       energies.Evdw = ljEnergy;
01619       energies.Eelec = electEnergy + electEnergySlow;
01620       energies.Ebond = bondEnergy;
01621       energies.Eangle = angleEnergy;
01622       energies.Edihe = dihedralEnergy + crosstermEnergy;
01623       energies.Eimpr = improperEnergy;
01624       Node::Object()->imd->gather_energies(&energies);
01625     }
01626 
01627     if ( marginViolations ) {
01628       iout << iERROR << marginViolations <<
01629         " margin violations detected since previous energy output.\n" << endi;
01630     }
01631     marginViolations = 0;
01632 
01633     if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
01634     {
01635         iout << "ETITLE:      TS";
01636         iout << FORMAT("BOND");
01637         iout << FORMAT("ANGLE");
01638         iout << FORMAT("DIHED");
01639         if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
01640         iout << FORMAT("IMPRP");
01641         iout << "     ";
01642         iout << FORMAT("ELECT");
01643         iout << FORMAT("VDW");
01644         iout << FORMAT("BOUNDARY");
01645         iout << FORMAT("MISC");
01646         iout << FORMAT("KINETIC");
01647         iout << "     ";
01648         iout << FORMAT("TOTAL");
01649         iout << FORMAT("TEMP");
01650         iout << FORMAT("POTENTIAL");
01651         // iout << FORMAT("TOTAL2");
01652         iout << FORMAT("TOTAL3");
01653         iout << FORMAT("TEMPAVG");
01654         if ( volume != 0. ) {
01655           iout << "     ";
01656           iout << FORMAT("PRESSURE");
01657           iout << FORMAT("GPRESSURE");
01658           iout << FORMAT("VOLUME");
01659           iout << FORMAT("PRESSAVG");
01660           iout << FORMAT("GPRESSAVG");
01661         }
01662         if (simParameters->drudeOn) {
01663           iout << "     ";
01664           iout << FORMAT("DRUDECOM");
01665           iout << FORMAT("DRUDEBOND");
01666           iout << FORMAT("DRCOMAVG");
01667           iout << FORMAT("DRBONDAVG");
01668         }
01669         iout << "\n\n" << endi;
01670     }
01671 
01672     // N.B.  HP's aCC compiler merges FORMAT calls in the same expression.
01673     //       Need separate statements because data returned in static array.
01674 
01675     iout << ETITLE(step);
01676     iout << FORMAT(bondEnergy);
01677     iout << FORMAT(angleEnergy);
01678     if ( simParameters->mergeCrossterms ) {
01679       iout << FORMAT(dihedralEnergy+crosstermEnergy);
01680     } else {
01681       iout << FORMAT(dihedralEnergy);
01682       iout << FORMAT(crosstermEnergy);
01683     }
01684     iout << FORMAT(improperEnergy);
01685     iout << "     ";
01686     iout << FORMAT(electEnergy+electEnergySlow);
01687     iout << FORMAT(ljEnergy);
01688     iout << FORMAT(boundaryEnergy);
01689     iout << FORMAT(miscEnergy);
01690     iout << FORMAT(kineticEnergy);
01691     iout << "     ";
01692     iout << FORMAT(totalEnergy);
01693     iout << FORMAT(temperature);
01694     iout << FORMAT(potentialEnergy);
01695     // iout << FORMAT(flatEnergy);
01696     iout << FORMAT(smoothEnergy);
01697     iout << FORMAT(temp_avg/avg_count);
01698     if ( volume != 0. )
01699     {
01700         iout << "     ";
01701         iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
01702         iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
01703         iout << FORMAT(volume);
01704         iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
01705         iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
01706     }
01707     if (simParameters->drudeOn) {
01708         iout << "     ";
01709         iout << FORMAT(drudeComTemp);
01710         iout << FORMAT(drudeBondTemp);
01711         iout << FORMAT(drudeComTempAvg/avg_count);
01712         iout << FORMAT(drudeBondTempAvg/avg_count);
01713     }
01714     iout << "\n\n" << endi;
01715 
01716 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
01717      char webout[80];
01718      sprintf(webout,"%d %d %d %d",(int)totalEnergy,
01719              (int)(potentialEnergy),
01720              (int)kineticEnergy,(int)temperature);
01721      CApplicationDepositNode0Data(webout);
01722 #endif
01723 
01724     if (simParameters->pairInteractionOn) {
01725       iout << "PAIR INTERACTION:";
01726       iout << " STEP: " << step;
01727       iout << " VDW_FORCE: ";
01728       iout << FORMAT(pairVDWForce.x);
01729       iout << FORMAT(pairVDWForce.y);
01730       iout << FORMAT(pairVDWForce.z);
01731       iout << " ELECT_FORCE: ";
01732       iout << FORMAT(pairElectForce.x);
01733       iout << FORMAT(pairElectForce.y);
01734       iout << FORMAT(pairElectForce.z);
01735       iout << "\n" << endi;
01736     }
01737     drudeComTempAvg = 0;
01738     drudeBondTempAvg = 0;
01739     temp_avg = 0;
01740     pressure_avg = 0;
01741     groupPressure_avg = 0;
01742     avg_count = 0;
01743 
01744     if ( fflush_count ) {
01745       --fflush_count;
01746       fflush(stdout);
01747     }
01748 }
01749 
01750 void Controller::writeExtendedSystemLabels(std::ofstream &file) {
01751   Lattice &lattice = state->lattice;
01752   file << "#$LABELS step";
01753   if ( lattice.a_p() ) file << " a_x a_y a_z";
01754   if ( lattice.b_p() ) file << " b_x b_y b_z";
01755   if ( lattice.c_p() ) file << " c_x c_y c_z";
01756   file << " o_x o_y o_z";
01757   if ( simParams->langevinPistonOn ) {
01758     file << " s_x s_y s_z s_u s_v s_w";
01759   }
01760   file << std::endl;
01761 }
01762 
01763 void Controller::writeExtendedSystemData(int step, std::ofstream &file) {
01764   Lattice &lattice = state->lattice;
01765   file << step;
01766     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
01767     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
01768     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
01769     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
01770   if ( simParams->langevinPistonOn ) {
01771     Vector strainRate = diagonal(langevinPiston_strainRate);
01772     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
01773     file << " " << strainRate.x;
01774     file << " " << strainRate.y;
01775     file << " " << strainRate.z;
01776     file << " " << strainRate2.x;
01777     file << " " << strainRate2.y;
01778     file << " " << strainRate2.z;
01779   }
01780   file << std::endl;
01781 }
01782 
01783 void Controller::enqueueCollections(int timestep)
01784 {
01785   if ( Output::coordinateNeeded(timestep) ){
01786     #ifdef MEM_OPT_VERSION
01787     EnqueueDataMsg *msg = new EnqueueDataMsg;
01788     msg->timestep = timestep;
01789     msg->l = state->lattice;
01790     collection->enqueuePositions(msg);
01791     #else
01792     collection->enqueuePositions(timestep,state->lattice);
01793     #endif
01794   }
01795   if ( Output::velocityNeeded(timestep) )
01796     collection->enqueueVelocities(timestep);
01797 }
01798 
01799 //Modifications for alchemical fep
01800 static char *FEPTITLE(int X)
01801 { 
01802   static char tmp_string[21];
01803   sprintf(tmp_string, "FepEnergy: %6d ",X);
01804   return tmp_string;
01805 }
01806 
01807 static char *TITITLE(int X)
01808 { 
01809   static char tmp_string[21];
01810   sprintf(tmp_string, "TI   %6d ",X);
01811   return tmp_string;
01812 }
01813 
01814 void Controller::outputFepEnergy(int step) {
01815  if (simParams->alchFepOn) {
01816   const int stepInRun = step - simParams->firstTimestep;
01817   const int alchEquilSteps = simParams->alchEquilSteps;
01818   const BigReal alchLambda = simParams->alchLambda;
01819   const BigReal alchLambda2 = simParams->alchLambda2;
01820   if (stepInRun == 0 || stepInRun == alchEquilSteps) {
01821     FepNo = 0;
01822     exp_dE_ByRT = 0.0;
01823     net_dE = 0.0;
01824   }
01825   BigReal dE = electEnergy_f + electEnergySlow_f + ljEnergy_f
01826                 - (electEnergy + electEnergySlow + ljEnergy);
01827   BigReal RT = BOLTZMAN * simParams->alchTemp;
01828   FepNo++;
01829   exp_dE_ByRT += exp(-dE/RT);
01830   net_dE += dE;
01831  
01832   if (stepInRun == 0) {
01833     if (!fepFile.rdbuf()->is_open()) {
01834       fepSum = 0.0;
01835       NAMD_backup_file(simParams->alchOutFile);
01836       fepFile.open(simParams->alchOutFile);
01837       iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
01838       fepFile << "#            STEP                 Elec                            "
01839               << "vdW                    dE           dE_avg         Temp             dG\n"
01840               << "#                           l             l+dl      "
01841               << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
01842     }
01843     fepFile << "#NEW FEP WINDOW: "
01844             << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
01845             << alchLambda2 << std::endl;
01846   }
01847   if (stepInRun == alchEquilSteps) {
01848     fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
01849             << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
01850             << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
01851   }
01852   if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) {
01853     writeFepEnergyData(step, fepFile);
01854     fepFile.flush();
01855   }
01856   if (step == simParams->N) {
01857     fepSum = fepSum + dG;
01858     fepFile << "#Free energy change for lambda window [ " << alchLambda
01859             << " " << alchLambda2 << " ] is " << dG << " ; net change until now is " << fepSum << std::endl;
01860   }
01861  }
01862 }
01863 
01864 void Controller::outputTiEnergy(int step) {
01865  if (simParams->alchThermIntOn) {
01866   const int stepInRun = step - simParams->firstTimestep;
01867   const int alchEquilSteps = simParams->alchEquilSteps;
01868   const BigReal alchLambda = simParams->alchLambda;
01869   
01870   if (stepInRun == 0 || stepInRun == alchEquilSteps) {
01871     TiNo = 0;
01872     net_dEdl_elec_1 = 0;
01873     net_dEdl_elec_2 = 0;
01874     net_dEdl_lj_1 = 0;
01875     net_dEdl_lj_2 = 0;
01876   }
01877   if (stepInRun == 0 || (! ((step - 1) % simParams->alchOutFreq))) {
01878     // output of instantaneous dU/dl now replaced with running average
01879     // over last alchOutFreq steps (except for step 0)
01880     recent_TiNo = 0;
01881     recent_dEdl_elec_1 = 0;
01882     recent_dEdl_elec_2 = 0;
01883     recent_dEdl_lj_1 = 0;
01884     recent_dEdl_lj_2 = 0;
01885   }
01886   TiNo++;
01887   recent_TiNo++;
01888   // FB - PME is no longer scaled by global lambda, but by the respective
01889   // lambda as dictated by elecLambdaStart. All electrostatics now go together.
01890   net_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1;
01891   net_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2;
01892   net_dEdl_lj_1 += ljEnergy_ti_1;
01893   net_dEdl_lj_2 += ljEnergy_ti_2;
01894   recent_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1; 
01895   recent_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2; 
01896   recent_dEdl_lj_1 += ljEnergy_ti_1;
01897   recent_dEdl_lj_2 += ljEnergy_ti_2;
01898 
01899   if (stepInRun == 0) {
01900     BigReal alchElecLambdaStart = simParams->alchElecLambdaStart;
01901     BigReal alchVdwLambdaEnd = simParams->alchVdwLambdaEnd;
01902     BigReal elec_lambda_1 = (alchLambda <= alchElecLambdaStart)? 0. : \
01903             (alchLambda - alchElecLambdaStart) / (1. - alchElecLambdaStart);
01904     BigReal elec_lambda_2 = ((1-alchLambda) <= alchElecLambdaStart)? 0. : \
01905             ((1-alchLambda) - alchElecLambdaStart) / (1. - alchElecLambdaStart);
01906     BigReal vdw_lambda_1 =  (alchLambda >= alchVdwLambdaEnd)? 1. : \
01907             alchLambda / alchVdwLambdaEnd; 
01908     BigReal vdw_lambda_2 =  ((1-alchLambda) >= alchVdwLambdaEnd)? 1. : \
01909             (1-alchLambda) / alchVdwLambdaEnd; 
01910     if (!tiFile.rdbuf()->is_open()) {
01911       //tiSum = 0.0;
01912       NAMD_backup_file(simParams->alchOutFile);
01913       tiFile.open(simParams->alchOutFile);
01914       iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
01915       tiFile << "#       STEP      Elec_dU/dl      Elec_avg        vdW_dU/dl      vdw_avg       Elec_dU/dl      Elec_avg      vdW_dU/dl       vdw_avg       PME_dU/dl      PME_avg\n"
01916               << "#               <---------------------PARTITION 1------------------------>    <---------------------PARTITION 2--------------------->" 
01917               << std::endl;
01918     }
01919     tiFile << "#NEW TI WINDOW: "
01920             << "LAMBDA " << alchLambda 
01921             << "\n#PARTITION 1 VDW LAMBDA " << vdw_lambda_1 
01922             << "\n#PARTITION 1 ELEC LAMBDA " << elec_lambda_1 
01923             << "\n#PARTITION 2 VDW LAMBDA " << vdw_lambda_2 
01924             << "\n#PARTITION 2 ELEC LAMBDA " << elec_lambda_2 
01925             << "\n" << std::endl;
01926   }
01927   if (stepInRun == alchEquilSteps) {
01928     tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
01929             << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
01930             << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
01931   }
01932   if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) {
01933     writeTiEnergyData(step, tiFile);
01934     tiFile.flush();
01935   }
01936  }
01937 }
01938 
01939 void Controller::writeFepEnergyData(int step, std::ofstream &file) {
01940   BigReal eeng = electEnergy+electEnergySlow;
01941   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
01942   BigReal dE = eeng_f + ljEnergy_f - eeng - ljEnergy;
01943   BigReal RT = BOLTZMAN * simParams->alchTemp;
01944   dG = -(RT * log(exp_dE_ByRT/FepNo));
01945   BigReal dE_avg = net_dE/FepNo;
01946   fepFile << FEPTITLE(step);
01947   fepFile << FORMAT(eeng);
01948   fepFile << FORMAT(eeng_f);
01949   fepFile << FORMAT(ljEnergy);
01950   fepFile << FORMAT(ljEnergy_f);
01951   fepFile << FORMAT(dE);
01952   fepFile << FORMAT(dE_avg);
01953   fepFile << FORMAT(temperature);
01954   fepFile << FORMAT(dG);
01955   fepFile << std::endl;
01956 }
01957 //fepe
01958 void Controller::writeTiEnergyData(int step, std::ofstream &file) {
01959   tiFile << TITITLE(step);
01960   tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
01961   tiFile << FORMAT(net_dEdl_elec_1/TiNo);
01962   tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
01963   tiFile << FORMAT(net_dEdl_lj_1/TiNo);
01964   tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
01965   tiFile << FORMAT(net_dEdl_elec_2/TiNo);
01966   tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
01967   tiFile << FORMAT(net_dEdl_lj_2/TiNo);
01968   tiFile << std::endl;
01969 }
01970 
01971 void Controller::outputExtendedSystem(int step)
01972 {
01973 
01974   if ( step >= 0 ) {
01975 
01976     // Write out eXtended System Trajectory (XST) file
01977     if ( simParams->xstFrequency &&
01978          ((step % simParams->xstFrequency) == 0) )
01979     {
01980       if ( ! xstFile.rdbuf()->is_open() )
01981       {
01982         NAMD_backup_file(simParams->xstFilename);
01983         xstFile.open(simParams->xstFilename);
01984         iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
01985         xstFile << "# NAMD extended system trajectory file" << std::endl;
01986         writeExtendedSystemLabels(xstFile);
01987       }
01988       writeExtendedSystemData(step,xstFile);
01989       xstFile.flush();
01990     }
01991 
01992     // Write out eXtended System Configuration (XSC) files
01993     //  Output a restart file
01994     if ( simParams->restartFrequency &&
01995          ((step % simParams->restartFrequency) == 0) &&
01996          (step != simParams->firstTimestep) )
01997     {
01998       iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
01999                 << step << "\n" << endi;
02000       char fname[140];
02001       strcpy(fname, simParams->restartFilename);
02002       if ( simParams->restartSave ) {
02003         char timestepstr[20];
02004         sprintf(timestepstr,".%d",step);
02005         strcat(fname, timestepstr);
02006       }
02007       strcat(fname, ".xsc");
02008       NAMD_backup_file(fname,".old");
02009       std::ofstream xscFile(fname);
02010       if (!xscFile) {
02011         char err_msg[257];
02012         sprintf(err_msg, "Error opening XSC restart file %s",fname);
02013         NAMD_err(err_msg);
02014       } 
02015       xscFile << "# NAMD extended system configuration restart file" << std::endl;
02016       writeExtendedSystemLabels(xscFile);
02017       writeExtendedSystemData(step,xscFile);
02018       if (!xscFile) {
02019         char err_msg[257];
02020         sprintf(err_msg, "Error writing XSC restart file %s",fname);
02021         NAMD_err(err_msg);
02022       } 
02023     }
02024 
02025   }
02026 
02027   //  Output final coordinates
02028   if (step == FILE_OUTPUT || step == END_OF_RUN)
02029   {
02030     iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
02031                 << simParams->N << "\n" << endi;
02032     static char fname[140];
02033     strcpy(fname, simParams->outputFilename);
02034     strcat(fname, ".xsc");
02035     NAMD_backup_file(fname);
02036     std::ofstream xscFile(fname);
02037     if (!xscFile) {
02038       char err_msg[257];
02039       sprintf(err_msg, "Error opening XSC output file %s",fname);
02040       NAMD_err(err_msg);
02041     } 
02042     xscFile << "# NAMD extended system configuration output file" << std::endl;
02043     writeExtendedSystemLabels(xscFile);
02044     writeExtendedSystemData(simParams->N,xscFile);
02045     if (!xscFile) {
02046       char err_msg[257];
02047       sprintf(err_msg, "Error writing XSC output file %s",fname);
02048       NAMD_err(err_msg);
02049     } 
02050   }
02051 
02052   //  Close trajectory file
02053   if (step == END_OF_RUN) {
02054     if ( xstFile.rdbuf()->is_open() ) {
02055       xstFile.close();
02056       iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
02057     }
02058   }
02059 
02060 }
02061 
02062 void Controller::rebalanceLoad(int)
02063 {
02064   if ( ! ldbSteps ) { 
02065     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
02066   }
02067   if ( ! --ldbSteps ) {
02068     startBenchTime -= CmiWallTimer();
02069     LdbCoordinator::Object()->rebalance(this);
02070     startBenchTime += CmiWallTimer();
02071     fflush_count = 3;
02072   }
02073 }
02074 
02075 void Controller::cycleBarrier(int doBarrier, int step) {
02076 #if USE_BARRIER
02077         if (doBarrier) {
02078           broadcast->cycleBarrier.publish(step,1);
02079           CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
02080                   CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02081         }
02082 #endif
02083 }
02084 
02085 void Controller::terminate(void) {
02086   BackEnd::awaken();
02087   CthFree(thread);
02088   CthSuspend();
02089 }
02090