Controller Class Reference

#include <Controller.h>

Inheritance diagram for Controller:
ControllerState

List of all members.

Classes

struct  checkpoint

Public Member Functions

 Controller (NamdState *s)
virtual ~Controller (void)
void run (void)
void awaken (void)
void resumeAfterTraceBarrier (int)

Protected Member Functions

virtual void algorithm (void)
void integrate (int)
void minimize ()
void receivePressure (int step, int minimize=0)
void calcPressure (int step, int minimize, const Tensor &virial_normal_in, const Tensor &virial_nbond_in, const Tensor &virial_slow_in, const Tensor &intVirial_normal, const Tensor &intVirial_nbond, const Tensor &intVirial_slow, const Vector &extForce_normal, const Vector &extForce_nbond, const Vector &extForce_slow)
void compareChecksums (int, int=0)
void printTiming (int)
void printMinimizeEnergies (int)
void printDynamicsEnergies (int)
void printEnergies (int step, int minimize)
void printFepMessage (int)
void printTiMessage (int)
void enqueueCollections (int)
void correctMomentum (int step)
void rescaleVelocities (int)
void reassignVelocities (int)
void tcoupleVelocities (int)
void stochRescaleVelocities (int)
void berendsenPressure (int)
void langevinPiston1 (int)
void langevinPiston2 (int)
void multigratorPressure (int step, int callNumber)
void multigratorTemperature (int step, int callNumber)
BigReal multigatorCalcEnthalpy (BigReal potentialEnergy, int step, int minimize)
void rebalanceLoad (int)
void cycleBarrier (int, int)
void traceBarrier (int, int)
void terminate (void)
void outputExtendedSystem (int step)
void writeExtendedSystemLabels (ofstream_namd &file)
void writeExtendedSystemData (int step, ofstream_namd &file)
void outputFepEnergy (int step)
void writeFepEnergyData (int step, ofstream_namd &file)
void outputTiEnergy (int step)
BigReal computeAlchWork (const int step)
void writeTiEnergyData (int step, ofstream_namd &file)
void recvCheckpointReq (const char *key, int task, checkpoint &cp)
void recvCheckpointAck (checkpoint &cp)
void calc_accelMDG_mean_std (BigReal testV, int step_n, BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV)
void calc_accelMDG_E_k (int iE, int V_n, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal *k0, BigReal *k, BigReal *E, int *iEused, char *warn)
void calc_accelMDG_force_factor (BigReal k, BigReal E, BigReal testV, Tensor vir_orig, BigReal *dV, BigReal *factor, Tensor *vir)
void write_accelMDG_rest_file (int step_n, char type, int V_n, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal M2, BigReal E, BigReal k, bool write_topic, bool lasttime)
void rescaleaccelMD (int step, int minimize=0)
void adaptTempInit (int step)
void adaptTempUpdate (int step, int minimize=0)
void adaptTempWriteRestart (int step)

Protected Attributes

RequireReductionmin_reduction
Tensor pressure_normal
Tensor pressure_nbond
Tensor pressure_slow
Tensor pressure_amd
Tensor virial_amd
Tensor groupPressure_normal
Tensor groupPressure_nbond
Tensor groupPressure_slow
Tensor controlPressure_normal
Tensor controlPressure_nbond
Tensor controlPressure_slow
int nbondFreq
int slowFreq
BigReal temp_avg
BigReal pressure_avg
BigReal groupPressure_avg
int avg_count
Tensor pressure_tavg
Tensor groupPressure_tavg
int tavg_count
int computeChecksum
int marginViolations
int pairlistWarnings
BigReal min_energy
BigReal min_f_dot_f
BigReal min_f_dot_v
BigReal min_v_dot_v
int min_huge_count
int64_t numDegFreedom
int stepInFullRun
BigReal totalEnergy
BigReal electEnergy
BigReal electEnergySlow
BigReal ljEnergy
BigReal groLJEnergy
BigReal groGaussEnergy
BigReal goNativeEnergy
BigReal goNonnativeEnergy
BigReal goTotalEnergy
BigReal bondedEnergyDiff_f
BigReal electEnergy_f
BigReal electEnergySlow_f
BigReal ljEnergy_f
BigReal ljEnergy_f_left
BigReal exp_dE_ByRT
BigReal dE
BigReal net_dE
BigReal dG
int FepNo
BigReal fepSum
BigReal bondedEnergy_ti_1
BigReal bondedEnergy_ti_2
BigReal electEnergy_ti_1
BigReal electEnergySlow_ti_1
BigReal ljEnergy_ti_1
BigReal electEnergy_ti_2
BigReal electEnergySlow_ti_2
BigReal ljEnergy_ti_2
BigReal net_dEdl_bond_1
BigReal net_dEdl_bond_2
BigReal net_dEdl_elec_1
BigReal net_dEdl_elec_2
BigReal net_dEdl_lj_1
BigReal net_dEdl_lj_2
BigReal cumAlchWork
BigReal electEnergyPME_ti_1
BigReal electEnergyPME_ti_2
int TiNo
BigReal recent_dEdl_bond_1
BigReal recent_dEdl_bond_2
BigReal recent_dEdl_elec_1
BigReal recent_dEdl_elec_2
BigReal recent_dEdl_lj_1
BigReal recent_dEdl_lj_2
BigReal recent_alchWork
BigReal alchWork
int recent_TiNo
BigReal drudeBondTemp
BigReal drudeBondTempAvg
BigReal kineticEnergy
BigReal kineticEnergyHalfstep
BigReal kineticEnergyCentered
BigReal temperature
BigReal heat
BigReal totalEnergy0
BigReal smooth2_avg2
Tensor pressure
Tensor groupPressure
int controlNumDegFreedom
Tensor controlPressure
BigReal rescaleVelocities_sumTemps
int rescaleVelocities_numTemps
int stochRescale_count
BigReal stochRescaleTimefactor
Tensor langevinPiston_origStrainRate
Tensor strainRate_old
Tensor positionRescaleFactor
BigReal multigratorXi
BigReal multigratorXiT
Tensor momentumSqrSum
std::vector< BigRealmultigratorNu
std::vector< BigRealmultigratorNuT
std::vector< BigRealmultigratorOmega
std::vector< BigRealmultigratorZeta
RequireReductionmultigratorReduction
int ldbSteps
int fflush_count
Randomrandom
SimParameters *const simParams
NamdState *const state
RequireReductionreduction
RequireReductionamd_reduction
SubmitReductionsubmit_reduction
PressureProfileReductionppbonded
PressureProfileReductionppnonbonded
PressureProfileReductionppint
int pressureProfileSlabs
int pressureProfileCount
BigRealpressureProfileAverage
CollectionMaster *const collection
ControllerBroadcastsbroadcast
ofstream_namd xstFile
ofstream_namd fepFile
ofstream_namd tiFile
int checkpoint_stored
Lattice checkpoint_lattice
ControllerState checkpoint_state
std::map< std::string,
checkpoint * > 
checkpoints
int checkpoint_task
Lattice origLattice
BigReal accelMDdVAverage
BigRealadaptTempPotEnergyAveNum
BigRealadaptTempPotEnergyAveDen
BigRealadaptTempPotEnergyVarNum
BigRealadaptTempPotEnergyAve
BigRealadaptTempPotEnergyVar
int * adaptTempPotEnergySamples
BigRealadaptTempBetaN
BigReal adaptTempT
BigReal adaptTempDTave
BigReal adaptTempDTavenum
BigReal adaptTempBetaMin
BigReal adaptTempBetaMax
int adaptTempBin
int adaptTempBins
BigReal adaptTempDBeta
BigReal adaptTempCg
BigReal adaptTempDt
Bool adaptTempAutoDt
BigReal adaptTempDtMin
BigReal adaptTempDtMax
ofstream_namd adaptTempRestartFile

Friends

class ScriptTcl
class Node
class CheckpointMsg

Detailed Description

Definition at line 39 of file Controller.h.


Constructor & Destructor Documentation

Controller::Controller ( NamdState s  ) 

Definition at line 128 of file Controller.C.

References SimParameters::accelMDOn, amd_reduction, avg_count, AVGXY, ControllerState::berendsenPressure_avg, ControllerState::berendsenPressure_count, BOLTZMANN, broadcast, checkpoint_stored, cumAlchWork, drudeBondTemp, drudeBondTempAvg, SimParameters::dt, groupPressure_avg, groupPressure_tavg, heat, Tensor::identity(), langevinPiston_origStrainRate, ControllerState::langevinPiston_strainRate, min_reduction, Node::molecule, MULTIGRATOR_REDUCTION_MAX_RESERVED, SimParameters::multigratorNoseHooverChainLength, multigratorNu, multigratorNuT, multigratorOmega, SimParameters::multigratorOn, multigratorReduction, SimParameters::multigratorTemperatureRelaxationTime, SimParameters::multigratorTemperatureTarget, multigratorXi, multigratorZeta, Molecule::num_deg_freedom(), numPatches, Node::Object(), PatchMap::Object(), ReductionMgr::Object(), origLattice, ppbonded, ppint, ppnonbonded, pressure_avg, pressure_tavg, SimParameters::pressureProfileAtomTypes, pressureProfileAverage, pressureProfileCount, SimParameters::pressureProfileEwaldOn, SimParameters::pressureProfileFreq, SimParameters::pressureProfileOn, SimParameters::pressureProfileSlabs, pressureProfileSlabs, random, SimParameters::randomSeed, reduction, REDUCTIONS_AMD, REDUCTIONS_BASIC, REDUCTIONS_MINIMIZER, REDUCTIONS_MULTIGRATOR, REDUCTIONS_PPROF_BONDED, REDUCTIONS_PPROF_INTERNAL, REDUCTIONS_PPROF_NONBONDED, rescaleVelocities_numTemps, rescaleVelocities_sumTemps, simParams, ControllerState::smooth2_avg, Random::split(), state, stochRescale_count, SimParameters::stochRescaleFreq, SimParameters::stochRescaleOn, SimParameters::stochRescalePeriod, stochRescaleTimefactor, SimParameters::strainRate, SimParameters::strainRate2, submit_reduction, Tensor::symmetric(), tavg_count, temp_avg, totalEnergy0, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, ReductionMgr::willRequire(), ReductionMgr::willSubmit(), and XXXBIGREAL.

00128                                    :
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     heat = totalEnergy0 = 0.0;
00187     rescaleVelocities_sumTemps = 0;  
00188     rescaleVelocities_numTemps = 0;
00189     stochRescale_count = 0;
00190     if (simParams->stochRescaleOn) {
00191       stochRescaleTimefactor = exp(-simParams->stochRescaleFreq *
00192           simParams->dt * 0.001 / simParams->stochRescalePeriod);
00193     }
00194     berendsenPressure_avg = 0; berendsenPressure_count = 0;
00195     // strainRate tensor is symmetric to avoid rotation
00196     langevinPiston_strainRate =
00197         Tensor::symmetric(simParams->strainRate,simParams->strainRate2);
00198     if ( ! simParams->useFlexibleCell ) {
00199       BigReal avg = trace(langevinPiston_strainRate) / 3.;
00200       langevinPiston_strainRate = Tensor::identity(avg);
00201     } else if ( simParams->useConstantRatio ) {
00202 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
00203                  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
00204       AVGXY(langevinPiston_strainRate);
00205 #undef AVGXY
00206     }
00207     langevinPiston_origStrainRate = langevinPiston_strainRate;
00208     if (simParams->multigratorOn) {
00209       multigratorXi = 0.0;
00210       int n = simParams->multigratorNoseHooverChainLength;
00211       BigReal tau = simParams->multigratorTemperatureRelaxationTime;
00212       Node *node = Node::Object();
00213       Molecule *molecule = node->molecule;
00214       BigReal Nf = molecule->num_deg_freedom();
00215       BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00216       multigratorNu.resize(n);
00217       multigratorNuT.resize(n);
00218       multigratorZeta.resize(n);
00219       multigratorOmega.resize(n);
00220       for (int i=0;i < n;i++) {
00221         multigratorNu[i] = 0.0;
00222         multigratorZeta[i] = 0.0;
00223         if (i == 0) {
00224           multigratorOmega[i] = Nf*kT0*tau*tau;
00225         } else {
00226           multigratorOmega[i] = kT0*tau*tau;
00227         }
00228       }
00229       multigratorReduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
00230     } else {
00231       multigratorReduction = NULL;
00232     }
00233     origLattice = state->lattice;
00234     smooth2_avg = XXXBIGREAL;
00235     temp_avg = 0;
00236     pressure_avg = 0;
00237     groupPressure_avg = 0;
00238     avg_count = 0;
00239     pressure_tavg = 0;
00240     groupPressure_tavg = 0;
00241     tavg_count = 0;
00242     checkpoint_stored = 0;
00243     drudeBondTemp = 0;
00244     drudeBondTempAvg = 0;
00245     cumAlchWork = 0;
00246 }

Controller::~Controller ( void   )  [virtual]

Definition at line 248 of file Controller.C.

References amd_reduction, broadcast, min_reduction, multigratorReduction, ppbonded, ppint, ppnonbonded, pressureProfileAverage, random, reduction, and submit_reduction.

00249 {
00250     delete broadcast;
00251     delete reduction;
00252     delete min_reduction;
00253     delete amd_reduction;
00254     delete submit_reduction;
00255     delete ppbonded;
00256     delete ppnonbonded;
00257     delete ppint;
00258     delete [] pressureProfileAverage;
00259     delete random;
00260     if (multigratorReduction) delete multigratorReduction;
00261 }


Member Function Documentation

void Controller::adaptTempInit ( int  step  )  [protected]

Definition at line 2331 of file Controller.C.

References SimParameters::adaptTempAutoDt, adaptTempAutoDt, adaptTempBetaMax, adaptTempBetaMin, adaptTempBetaN, adaptTempBins, SimParameters::adaptTempBins, adaptTempCg, SimParameters::adaptTempCgamma, adaptTempDBeta, SimParameters::adaptTempDt, adaptTempDt, adaptTempDTave, adaptTempDTavenum, adaptTempDtMax, adaptTempDtMin, SimParameters::adaptTempInFile, SimParameters::adaptTempOn, adaptTempPotEnergyAve, adaptTempPotEnergyAveDen, adaptTempPotEnergyAveNum, adaptTempPotEnergySamples, adaptTempPotEnergyVar, adaptTempPotEnergyVarNum, adaptTempRestartFile, SimParameters::adaptTempRestartFile, SimParameters::adaptTempRestartFreq, adaptTempT, SimParameters::adaptTempTmax, SimParameters::adaptTempTmin, endi(), iINFO(), SimParameters::initialTemp, iout, j, SimParameters::langevinOn, SimParameters::langevinTemp, NAMD_backup_file(), NAMD_die(), ofstream_namd::open(), SimParameters::rescaleFreq, SimParameters::rescaleTemp, and simParams.

Referenced by adaptTempUpdate().

02331                                        {
02332     if (!simParams->adaptTempOn) return;
02333     iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
02334     adaptTempDtMin = 0;
02335     adaptTempDtMax = 0;
02336     adaptTempAutoDt = false;
02337     if (simParams->adaptTempBins == 0) {
02338       iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
02339       std::ifstream adaptTempRead(simParams->adaptTempInFile);
02340       if (adaptTempRead) {
02341         int readInt;
02342         BigReal readReal;
02343         bool readBool;
02344         // step
02345         adaptTempRead >> readInt;
02346         // Start with min and max temperatures
02347         adaptTempRead >> adaptTempT;     // KELVIN
02348         adaptTempRead >> adaptTempBetaMin;  // KELVIN
02349         adaptTempRead >> adaptTempBetaMax;  // KELVIN
02350         adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
02351         adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
02352         // In case file is manually edited
02353         if (adaptTempBetaMin > adaptTempBetaMax){
02354             readReal = adaptTempBetaMax;
02355             adaptTempBetaMax = adaptTempBetaMin;
02356             adaptTempBetaMin = adaptTempBetaMax;
02357         }
02358         adaptTempRead >> adaptTempBins;     
02359         adaptTempRead >> adaptTempCg;
02360         adaptTempRead >> adaptTempDt;
02361         adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02362         adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02363         adaptTempPotEnergySamples = new int[adaptTempBins];
02364         adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02365         adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02366         adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02367         adaptTempBetaN           = new BigReal[adaptTempBins];
02368         adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02369         for(int j = 0; j < adaptTempBins; ++j) {
02370           adaptTempRead >> adaptTempPotEnergyAve[j];
02371           adaptTempRead >> adaptTempPotEnergyVar[j];
02372           adaptTempRead >> adaptTempPotEnergySamples[j];
02373           adaptTempRead >> adaptTempPotEnergyAveNum[j];
02374           adaptTempRead >> adaptTempPotEnergyVarNum[j];
02375           adaptTempRead >> adaptTempPotEnergyAveDen[j];
02376           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02377         } 
02378         adaptTempRead.close();
02379       }
02380       else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
02381     } 
02382     else {
02383       adaptTempBins = simParams->adaptTempBins;
02384       adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02385       adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02386       adaptTempPotEnergySamples = new int[adaptTempBins];
02387       adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02388       adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02389       adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02390       adaptTempBetaN           = new BigReal[adaptTempBins];
02391       adaptTempBetaMax = 1./simParams->adaptTempTmin;
02392       adaptTempBetaMin = 1./simParams->adaptTempTmax;
02393       adaptTempCg = simParams->adaptTempCgamma;   
02394       adaptTempDt  = simParams->adaptTempDt;
02395       adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02396       adaptTempT = simParams->initialTemp; 
02397       if (simParams->langevinOn)
02398         adaptTempT = simParams->langevinTemp;
02399       else if (simParams->rescaleFreq > 0)
02400         adaptTempT = simParams->rescaleTemp;
02401       for(int j = 0; j < adaptTempBins; ++j){
02402           adaptTempPotEnergyAveNum[j] = 0.;
02403           adaptTempPotEnergyAveDen[j] = 0.;
02404           adaptTempPotEnergySamples[j] = 0;
02405           adaptTempPotEnergyVarNum[j] = 0.;
02406           adaptTempPotEnergyVar[j] = 0.;
02407           adaptTempPotEnergyAve[j] = 0.;
02408           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02409       }
02410     }
02411     if (simParams->adaptTempAutoDt > 0.0) {
02412        adaptTempAutoDt = true;
02413        adaptTempDtMin =  simParams->adaptTempAutoDt - 0.01;
02414        adaptTempDtMax =  simParams->adaptTempAutoDt + 0.01;
02415     }
02416     adaptTempDTave = 0;
02417     adaptTempDTavenum = 0;
02418     iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
02419     iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
02420     iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
02421     iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
02422     iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
02423     if (simParams->adaptTempRestartFreq > 0) {
02424       NAMD_backup_file(simParams->adaptTempRestartFile);
02425       adaptTempRestartFile.open(simParams->adaptTempRestartFile);
02426        if(!adaptTempRestartFile)
02427         NAMD_die("Error opening restart file");
02428     }
02429 }

void Controller::adaptTempUpdate ( int  step,
int  minimize = 0 
) [protected]

Definition at line 2457 of file Controller.C.

References adaptTempAutoDt, adaptTempBetaMax, adaptTempBetaMin, adaptTempBetaN, adaptTempBin, adaptTempBins, adaptTempCg, adaptTempDBeta, SimParameters::adaptTempDebug, adaptTempDt, adaptTempDTave, adaptTempDTavenum, adaptTempDtMax, adaptTempDtMin, ControllerBroadcasts::adaptTemperature, SimParameters::adaptTempFreq, adaptTempInit(), SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, SimParameters::adaptTempOutFreq, adaptTempPotEnergyAve, adaptTempPotEnergyAveDen, adaptTempPotEnergyAveNum, adaptTempPotEnergySamples, adaptTempPotEnergyVar, adaptTempPotEnergyVarNum, SimParameters::adaptTempRandom, adaptTempT, adaptTempWriteRestart(), BOLTZMANN, broadcast, endi(), SimParameters::firstTimestep, Random::gaussian(), iout, iWARN(), j, kineticEnergy, SimpleBroadcastObject< T >::publish(), random, simParams, totalEnergy, and Random::uniform().

Referenced by integrate().

02458 {
02459     //Beta = 1./T
02460     if ( !simParams->adaptTempOn ) return;
02461     int j = 0;
02462     if (step == simParams->firstTimestep) {
02463         adaptTempInit(step);
02464         return;
02465     }
02466     if ( minimize || (step < simParams->adaptTempFirstStep ) || 
02467         ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
02468     const int adaptTempOutFreq  = simParams->adaptTempOutFreq;
02469     const bool adaptTempDebug  = simParams->adaptTempDebug;
02470     //Calculate Current inverse temperature and bin 
02471     BigReal adaptTempBeta = 1./adaptTempT;
02472     adaptTempBin   = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
02473 
02474     if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
02475         iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin  
02476                                << " adaptTempBeta: " << adaptTempBeta 
02477                               << " adaptTempDBeta: " << adaptTempDBeta 
02478                                << " betaMin:" << adaptTempBetaMin 
02479                                << " betaMax: " << adaptTempBetaMax << "\n";
02480     adaptTempPotEnergySamples[adaptTempBin] += 1;
02481     BigReal gammaAve = 1.-adaptTempCg/adaptTempPotEnergySamples[adaptTempBin];
02482 
02483     BigReal potentialEnergy;
02484     BigReal potEnergyAverage;
02485     BigReal potEnergyVariance;
02486     potentialEnergy = totalEnergy-kineticEnergy;
02487 
02488     //calculate new bin average and variance using adaptive averaging
02489     adaptTempPotEnergyAveNum[adaptTempBin] = adaptTempPotEnergyAveNum[adaptTempBin]*gammaAve + potentialEnergy;
02490     adaptTempPotEnergyAveDen[adaptTempBin] = adaptTempPotEnergyAveDen[adaptTempBin]*gammaAve + 1;
02491     adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
02492     
02493     potEnergyAverage = adaptTempPotEnergyAveNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
02494     potEnergyVariance = adaptTempPotEnergyVarNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
02495     potEnergyVariance -= potEnergyAverage*potEnergyAverage;
02496 
02497     adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
02498     adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
02499     
02500     // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
02501     // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
02502     // is constant for each bin. This is to estimate <E(beta)> where beta \in
02503     // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
02504     if ( ! ( step % simParams->adaptTempFreq ) ) {
02505       // If adaptTempBin not at the edge, calculate improved average:
02506       if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
02507           // Get Averaging Limits:
02508           BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
02509           BigReal betaPlus;
02510           BigReal betaMinus;
02511           int     nMinus =0;
02512           int     nPlus = 0;
02513           if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
02514           if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
02515           betaMinus = adaptTempBeta - deltaBeta;
02516           betaPlus =  adaptTempBeta + deltaBeta;
02517           BigReal betaMinus2 = betaMinus*betaMinus;
02518           BigReal betaPlus2  = betaPlus*betaPlus;
02519           nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
02520           nPlus  = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
02521           // Variables for <E(beta)> estimate:
02522           BigReal potEnergyAve0 = 0.0;
02523           BigReal potEnergyAve1 = 0.0;
02524           // Integral terms
02525           BigReal A0 = 0;
02526           BigReal A1 = 0;
02527           BigReal A2 = 0;
02528           //A0 phi_s integral for beta_minus < beta < beta_{i+1}
02529           BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1]; 
02530           BigReal DeltaE2Ave;
02531           BigReal DeltaE2AveSum = 0;
02532           BigReal betaj;
02533           for (j = nMinus+1; j <= adaptTempBin; ++j) {
02534               potEnergyAve0 += adaptTempPotEnergyAve[j]; 
02535               DeltaE2Ave = adaptTempPotEnergyVar[j];
02536               DeltaE2AveSum += DeltaE2Ave;
02537               A0 += j*DeltaE2Ave;
02538           }
02539           A0 *= adaptTempDBeta;
02540           A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
02541           A0 *= adaptTempDBeta;
02542           betaj = adaptTempBetaN[nMinus+1]-betaMinus; 
02543           betaj *= betaj;
02544           A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
02545           A0 /= (betaNp1-betaMinus);
02546 
02547           //A1 phi_s integral for beta_{i+1} < beta < beta_plus
02548           DeltaE2AveSum = 0;
02549           for (j = adaptTempBin+1; j < nPlus; j++) {
02550               potEnergyAve1 += adaptTempPotEnergyAve[j];
02551               DeltaE2Ave = adaptTempPotEnergyVar[j];
02552               DeltaE2AveSum += DeltaE2Ave;
02553               A1 += j*DeltaE2Ave;
02554           }
02555           A1 *= adaptTempDBeta;   
02556           A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
02557           A1 *= adaptTempDBeta;
02558           if ( nPlus < adaptTempBins ) {
02559             betaj = betaPlus-adaptTempBetaN[nPlus];
02560             betaj *= betaj;
02561             A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
02562           }
02563           A1 /= (betaPlus-betaNp1);
02564 
02565           //A2 phi_t integral for beta_i
02566           A2 = 0.5*adaptTempDBeta*potEnergyVariance;
02567 
02568           // Now calculate a+ and a-
02569           DeltaE2Ave = A0-A1;
02570           BigReal aplus = 0;
02571           BigReal aminus = 0;
02572           if (DeltaE2Ave != 0) {
02573             aplus = (A0-A2)/(A0-A1);
02574             if (aplus < 0) {
02575                     aplus = 0;
02576             }
02577             if (aplus > 1)  {
02578                     aplus = 1;
02579             }
02580             aminus = 1-aplus;
02581             potEnergyAve0 *= adaptTempDBeta;
02582             potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
02583             potEnergyAve0 /= (betaNp1-betaMinus);
02584             potEnergyAve1 *= adaptTempDBeta;
02585             if ( nPlus < adaptTempBins ) {
02586                 potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
02587             }
02588             potEnergyAve1 /= (betaPlus-betaNp1);
02589             potEnergyAverage = aminus*potEnergyAve0;
02590             potEnergyAverage += aplus*potEnergyAve1;
02591           }
02592           if (simParams->adaptTempDebug) {
02593        iout << "ADAPTEMP DEBUG:"  << "\n"
02594             << "     adaptTempBin:    " << adaptTempBin << "\n"
02595             << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
02596             << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
02597             << "     adaptTemp:   " << adaptTempT<< "\n"
02598             << "     betaMin:   " << adaptTempBetaMin << "\n"
02599             << "     betaMax:   " << adaptTempBetaMax << "\n"
02600             << "     gammaAve:  " << gammaAve << "\n"
02601             << "     deltaBeta: " << deltaBeta << "\n"
02602             << "     betaMinus: " << betaMinus << "\n"
02603             << "     betaPlus:  " << betaPlus << "\n"
02604             << "     nMinus:    " << nMinus << "\n"
02605             << "     nPlus:     " << nPlus << "\n"
02606             << "     A0:        " << A0 << "\n"
02607             << "     A1:        " << A1 << "\n"
02608             << "     A2:        " << A2 << "\n"
02609             << "     a+:        " << aplus << "\n"
02610             << "     a-:        " << aminus << "\n"
02611             << endi;
02612           }
02613       }
02614       else {
02615           if (simParams->adaptTempDebug) {
02616        iout << "ADAPTEMP DEBUG:"  << "\n"
02617             << "     adaptTempBin:    " << adaptTempBin << "\n"
02618             << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
02619             << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
02620             << "     adaptTemp:   " << adaptTempT<< "\n"
02621             << "     betaMin:   " << adaptTempBetaMin << "\n"
02622             << "     betaMax:   " << adaptTempBetaMax << "\n"
02623             << "     gammaAve:  " << gammaAve << "\n"
02624             << "     deltaBeta:  N/A\n"
02625             << "     betaMinus:  N/A\n"
02626             << "     betaPlus:   N/A\n"
02627             << "     nMinus:     N/A\n"
02628             << "     nPlus:      N/A\n"
02629             << "     A0:         N/A\n"
02630             << "     A1:         N/A\n"
02631             << "     A2:         N/A\n"
02632             << "     a+:         N/A\n"
02633             << "     a-:         N/A\n"
02634             << endi;
02635           }
02636       }
02637       
02638       //dT is new temperature
02639       BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
02640       dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
02641       dT += adaptTempT;
02642       
02643      // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
02644      // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
02645       if ( dT > 1./adaptTempBetaMin || dT  < 1./adaptTempBetaMax ) {
02646         dT = ((potentialEnergy-adaptTempPotEnergyAve[adaptTempBin])/BOLTZMANN+adaptTempT)*adaptTempDt;
02647         dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
02648         dT += adaptTempT;
02649         // Check again, if not then keep original adaptTempTor assign random.
02650         if ( dT > 1./adaptTempBetaMin ) {
02651           if (!simParams->adaptTempRandom) {             
02652              //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
02653              //     << " K higher than adaptTempTmax."
02654              //     << " Keeping temperature at " 
02655              //     << adaptTempT<< "\n"<< endi;             
02656              dT = adaptTempT;
02657           }
02658           else {             
02659              //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
02660              //     << " K higher than adaptTempTmax."
02661              //     << " Assigning random temperature in range\n" << endi;
02662              dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);             
02663              dT = 1./dT;
02664           }
02665         } 
02666         else if ( dT  < 1./adaptTempBetaMax ) {
02667           if (!simParams->adaptTempRandom) {            
02668             //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
02669             //     << " K lower than adaptTempTmin."
02670             //     << " Keeping temperature at " << adaptTempT<< "\n" << endi; 
02671             dT = adaptTempT;
02672           }
02673           else {
02674             //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
02675             //     << " K lower than adaptTempTmin."
02676             //     << " Assigning random temperature in range\n" << endi;
02677             dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
02678             dT = 1./dT;
02679           }
02680         }
02681         else if (adaptTempAutoDt) {
02682           //update temperature step size counter
02683           //FOR "TRUE" ADAPTIVE TEMPERING 
02684           BigReal adaptTempTdiff = fabs(dT-adaptTempT);
02685           if (adaptTempTdiff > 0) {
02686             adaptTempDTave += adaptTempTdiff; 
02687             adaptTempDTavenum++;
02688 //            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
02689           }
02690           if(adaptTempDTavenum == 100){
02691                 BigReal Frac;
02692                 adaptTempDTave /= adaptTempDTavenum;
02693                 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
02694                 Frac = adaptTempDTave/Frac;
02695                 //if average temperature jump is > 3% of temperature range,
02696                 //modify jump size to match 3%
02697                 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
02698                 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
02699                     Frac = adaptTempDtMax/Frac;
02700                     iout << "ADAPTEMP: Updating adaptTempDt to ";
02701                     adaptTempDt *= Frac;
02702                     iout << adaptTempDt << "\n" << endi;
02703                 }
02704                 adaptTempDTave = 0;
02705                 adaptTempDTavenum = 0;
02706           }
02707         }
02708       }
02709       else if (adaptTempAutoDt) {
02710           //update temperature step size counter
02711           // FOR "TRUE" ADAPTIVE TEMPERING
02712           BigReal adaptTempTdiff = fabs(dT-adaptTempT);
02713           if (adaptTempTdiff > 0) {
02714             adaptTempDTave += adaptTempTdiff; 
02715             adaptTempDTavenum++;
02716 //            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
02717           }
02718           if(adaptTempDTavenum == 100){
02719                 BigReal Frac;
02720                 adaptTempDTave /= adaptTempDTavenum;
02721                 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
02722                 Frac = adaptTempDTave/Frac;
02723                 //if average temperature jump is > 3% of temperature range,
02724                 //modify jump size to match 3%
02725                 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
02726                 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
02727                     Frac = adaptTempDtMax/Frac;
02728                     iout << "ADAPTEMP: Updating adaptTempDt to ";
02729                     adaptTempDt *= Frac;
02730                     iout << adaptTempDt << "\n" << endi;
02731                 }
02732                 adaptTempDTave = 0;
02733                 adaptTempDTavenum = 0;
02734 
02735           }
02736           
02737       }
02738       
02739       adaptTempT = dT; 
02740       broadcast->adaptTemperature.publish(step,adaptTempT);
02741     }
02742     adaptTempWriteRestart(step);
02743     if ( ! (step % adaptTempOutFreq) ) {
02744         iout << "ADAPTEMP: STEP " << step
02745              << " TEMP "   << adaptTempT
02746              << " ENERGY " << std::setprecision(10) << potentialEnergy   
02747              << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
02748              << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
02749         iout << "\n" << endi;
02750    }
02751    
02752 }

void Controller::adaptTempWriteRestart ( int  step  )  [protected]

Definition at line 2431 of file Controller.C.

References adaptTempBetaMax, adaptTempBetaMin, adaptTempBins, adaptTempCg, adaptTempDt, SimParameters::adaptTempOn, adaptTempPotEnergyAve, adaptTempPotEnergyAveDen, adaptTempPotEnergyAveNum, adaptTempPotEnergySamples, adaptTempPotEnergyVar, adaptTempPotEnergyVarNum, adaptTempRestartFile, SimParameters::adaptTempRestartFreq, adaptTempT, endi(), ofstream_namd::flush(), iout, j, and simParams.

Referenced by adaptTempUpdate().

02431                                                {
02432     if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
02433         adaptTempRestartFile.seekp(std::ios::beg);        
02434         iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
02435         adaptTempRestartFile << step << " ";
02436         // Start with min and max temperatures
02437         adaptTempRestartFile << adaptTempT<< " ";     // KELVIN
02438         adaptTempRestartFile << 1./adaptTempBetaMin << " ";  // KELVIN
02439         adaptTempRestartFile << 1./adaptTempBetaMax << " ";  // KELVIN
02440         adaptTempRestartFile << adaptTempBins << " ";     
02441         adaptTempRestartFile << adaptTempCg << " ";
02442         adaptTempRestartFile << adaptTempDt ;
02443         adaptTempRestartFile << "\n" ;
02444         for(int j = 0; j < adaptTempBins; ++j) {
02445           adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
02446           adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
02447           adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
02448           adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
02449           adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
02450           adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
02451           adaptTempRestartFile << "\n";          
02452         }
02453         adaptTempRestartFile.flush(); 
02454     }
02455 }    

void Controller::algorithm ( void   )  [protected, virtual]

Definition at line 281 of file Controller.C.

References awaken(), broadcast, checkpoint_lattice, checkpoint_state, checkpoint_stored, checkpoints, END_OF_RUN, endi(), enqueueCollections(), EVAL_MEASURE, FILE_OUTPUT, SimParameters::firstTimestep, FORCE_OUTPUT, SimpleBroadcastObject< T >::get(), if(), SimParameters::initialTemp, integrate(), iout, Controller::checkpoint::lattice, minimize(), NAMD_bug(), NAMD_die(), Node::Object(), outputExtendedSystem(), SCRIPT_ATOMRECV, SCRIPT_ATOMSEND, SCRIPT_ATOMSENDRECV, SCRIPT_CHECKPOINT, SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_STORE, SCRIPT_CHECKPOINT_SWAP, SCRIPT_CONTINUE, SCRIPT_END, SCRIPT_FORCEOUTPUT, SCRIPT_MEASURE, SCRIPT_MINIMIZE, SCRIPT_OUTPUT, SCRIPT_REINITVELS, SCRIPT_RESCALESOLUTECHARGES, SCRIPT_RESCALEVELS, SCRIPT_REVERT, SCRIPT_RUN, SimParameters::scriptArg1, ControllerBroadcasts::scriptBarrier, SimParameters::scriptIntArg1, SimParameters::scriptStringArg1, Node::sendCheckpointReq(), simParams, Controller::checkpoint::state, state, msm::swap(), and terminate().

00282 {
00283   int scriptTask;
00284   int scriptSeq = 0;
00285   BackEnd::awaken();
00286   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
00287     switch ( scriptTask ) {
00288       case SCRIPT_OUTPUT:
00289         enqueueCollections(FILE_OUTPUT);
00290         outputExtendedSystem(FILE_OUTPUT);
00291         break;
00292       case SCRIPT_FORCEOUTPUT:
00293         enqueueCollections(FORCE_OUTPUT);
00294         break;
00295       case SCRIPT_MEASURE:
00296         enqueueCollections(EVAL_MEASURE);
00297         break;
00298       case SCRIPT_REINITVELS:
00299         iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
00300           << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
00301         break;
00302       case SCRIPT_RESCALEVELS:
00303         iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
00304           << " BY " << simParams->scriptArg1 << "\n" << endi;
00305         break;
00306       case SCRIPT_RESCALESOLUTECHARGES:
00307         // Parameter setting already reported in ScriptTcl
00308         // Nothing to do!
00309         break;
00310       case SCRIPT_CHECKPOINT:
00311         iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
00312           << "\n" << endi;
00313         checkpoint_stored = 1;
00314         checkpoint_lattice = state->lattice;
00315         checkpoint_state = *(ControllerState*)this;
00316         break;
00317       case SCRIPT_REVERT:
00318         iout << "REVERTING AT STEP " << simParams->firstTimestep
00319           << "\n" << endi;
00320         if ( ! checkpoint_stored )
00321           NAMD_die("Unable to revert, checkpoint was never called!");
00322         state->lattice = checkpoint_lattice;
00323         *(ControllerState*)this = checkpoint_state;
00324         break;
00325       case SCRIPT_CHECKPOINT_STORE:
00326         iout << "STORING CHECKPOINT AT STEP " << simParams->firstTimestep
00327           << " TO KEY " << simParams->scriptStringArg1;
00328         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00329         iout << "\n" << endi;
00330         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00331           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00332             checkpoints[simParams->scriptStringArg1] = new checkpoint;
00333           }
00334           checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
00335           cp.lattice = state->lattice;
00336           cp.state = *(ControllerState*)this;
00337         } else {
00338           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00339                                             scriptTask, state->lattice, *(ControllerState*)this);
00340         }
00341         break;
00342       case SCRIPT_CHECKPOINT_LOAD:
00343         iout << "LOADING CHECKPOINT AT STEP " << simParams->firstTimestep
00344           << " FROM KEY " << simParams->scriptStringArg1;
00345         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00346         iout << "\n" << endi;
00347         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00348           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00349             NAMD_die("Unable to load checkpoint, requested key was never stored.");
00350           }
00351           checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
00352           state->lattice = cp.lattice;
00353           *(ControllerState*)this = cp.state;
00354         } else {
00355           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00356                                             scriptTask, state->lattice, *(ControllerState*)this);
00357         }
00358         break;
00359       case SCRIPT_CHECKPOINT_SWAP:
00360         iout << "SWAPPING CHECKPOINT AT STEP " << simParams->firstTimestep
00361           << " FROM KEY " << simParams->scriptStringArg1;
00362         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00363         iout << "\n" << endi;
00364         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00365           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00366             NAMD_die("Unable to swap checkpoint, requested key was never stored.");
00367           }
00368           checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
00369           std::swap(state->lattice,cp.lattice);
00370           std::swap(*(ControllerState*)this,cp.state);
00371         } else {
00372           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00373                                             scriptTask, state->lattice, *(ControllerState*)this);
00374         }
00375         break;
00376       case SCRIPT_CHECKPOINT_FREE:
00377         iout << "FREEING CHECKPOINT AT STEP " << simParams->firstTimestep
00378           << " FROM KEY " << simParams->scriptStringArg1;
00379         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
00380         iout << "\n" << endi;
00381         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
00382           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
00383             NAMD_die("Unable to free checkpoint, requested key was never stored.");
00384           }
00385           delete checkpoints[simParams->scriptStringArg1];
00386           checkpoints.erase(simParams->scriptStringArg1);
00387         } else {
00388           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
00389                                             scriptTask, state->lattice, *(ControllerState*)this);
00390         }
00391         break;
00392       case SCRIPT_ATOMSENDRECV:
00393       case SCRIPT_ATOMSEND:
00394       case SCRIPT_ATOMRECV:
00395         break;
00396       case SCRIPT_MINIMIZE:
00397         minimize();
00398         break;
00399       case SCRIPT_RUN:
00400       case SCRIPT_CONTINUE:
00401         integrate(scriptTask);
00402         break;
00403       default:
00404         NAMD_bug("Unknown task in Controller::algorithm");
00405     }
00406     BackEnd::awaken();
00407   }
00408   enqueueCollections(END_OF_RUN);
00409   outputExtendedSystem(END_OF_RUN);
00410   terminate();
00411 }

void Controller::awaken ( void   )  [inline]

Definition at line 45 of file Controller.h.

Referenced by algorithm(), LdbCoordinator::awakenSequencers(), resumeAfterTraceBarrier(), run(), and terminate().

00045 { CthAwaken(thread); };

void Controller::berendsenPressure ( int  step  )  [protected]

Definition at line 950 of file Controller.C.

References ControllerState::berendsenPressure_avg, ControllerState::berendsenPressure_count, SimParameters::berendsenPressureCompressibility, SimParameters::berendsenPressureFreq, SimParameters::berendsenPressureOn, SimParameters::berendsenPressureRelaxationTime, SimParameters::berendsenPressureTarget, broadcast, controlPressure, Tensor::diagonal(), SimParameters::dt, endi(), Tensor::identity(), iERROR(), iout, LIMIT_SCALING, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), Lattice::rescale(), simParams, state, Tensor::xx, Tensor::yy, and Tensor::zz.

Referenced by integrate().

00951 {
00952   if ( simParams->berendsenPressureOn ) {
00953    berendsenPressure_count += 1;
00954    berendsenPressure_avg += controlPressure;
00955    const int freq = simParams->berendsenPressureFreq;
00956    if ( ! (berendsenPressure_count % freq) ) {
00957     Tensor factor = berendsenPressure_avg / berendsenPressure_count;
00958     berendsenPressure_avg = 0;
00959     berendsenPressure_count = 0;
00960     // We only use on-diagonal terms (for now)
00961     factor = Tensor::diagonal(diagonal(factor));
00962     factor -= Tensor::identity(simParams->berendsenPressureTarget);
00963     factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
00964        simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
00965     factor += Tensor::identity(1.0);
00966 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
00967          if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
00968          if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
00969     int limited = 0;
00970     LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
00971     LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
00972     LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
00973 #undef LIMIT_SCALING
00974     if ( limited ) {
00975       iout << iERROR << "Step " << step <<
00976         " cell rescaling factor limited.\n" << endi;
00977     }
00978     broadcast->positionRescaleFactor.publish(step,factor);
00979     state->lattice.rescale(factor);
00980    }
00981   } else {
00982     berendsenPressure_avg = 0;
00983     berendsenPressure_count = 0;
00984   }
00985 }

void Controller::calc_accelMDG_E_k ( int  iE,
int  V_n,
BigReal  sigma0,
BigReal  Vmax,
BigReal  Vmin,
BigReal  Vavg,
BigReal  sigmaV,
BigReal k0,
BigReal k,
BigReal E,
int *  iEused,
char *  warn 
) [inline, protected]

Definition at line 1720 of file Controller.C.

Referenced by rescaleaccelMD().

01721                                                               {
01722    switch(iE){
01723       case 2:
01724          *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
01725          // if k0 <=0 OR k0 > 1, switch to iE=1 mode
01726          if(*k0 > 1.)
01727             *warn |= 1;
01728          else if(*k0 <= 0.)
01729             *warn |= 2;
01730          //else stay in iE=2 mode
01731          else{
01732             *E = Vmin + (Vmax-Vmin)/(*k0);
01733             *iEused = 2;
01734             break;
01735          }
01736       case 1:
01737          *E = Vmax;
01738          *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
01739          if(*k0 > 1.)
01740             *k0 = 1.;
01741          *iEused = 1;
01742          break;
01743    }
01744 
01745    *k = *k0 / (Vmax - Vmin);
01746 }

void Controller::calc_accelMDG_force_factor ( BigReal  k,
BigReal  E,
BigReal  testV,
Tensor  vir_orig,
BigReal dV,
BigReal factor,
Tensor vir 
) [inline, protected]

Definition at line 1749 of file Controller.C.

Referenced by rescaleaccelMD().

01750                                            {
01751    BigReal VE = testV - E;
01752    //if V < E, apply boost
01753    if(VE < 0){
01754       *factor = k * VE;
01755       *vir += vir_orig * (*factor);
01756       *dV += (*factor) * VE/2.;
01757       *factor += 1.;
01758    }
01759    else{
01760       *factor = 1.;
01761    }
01762 }

void Controller::calc_accelMDG_mean_std ( BigReal  testV,
int  step_n,
BigReal Vmax,
BigReal Vmin,
BigReal Vavg,
BigReal M2,
BigReal sigmaV 
) [inline, protected]

Definition at line 1700 of file Controller.C.

Referenced by rescaleaccelMD().

01701                                                                            {
01702    BigReal delta; 
01703 
01704    if(testV > *Vmax){
01705       *Vmax = testV;
01706    }
01707    else if(testV < *Vmin){
01708       *Vmin = testV;
01709    }
01710 
01711    //mean and std calculated by Online Algorithm
01712    delta = testV - *Vavg;
01713    *Vavg += delta / (BigReal)n;
01714    *M2 += delta * (testV - (*Vavg));
01715 
01716    *sigmaV = sqrt(*M2/n);
01717 }

void Controller::calcPressure ( int  step,
int  minimize,
const Tensor virial_normal_in,
const Tensor virial_nbond_in,
const Tensor virial_slow_in,
const Tensor intVirial_normal,
const Tensor intVirial_nbond,
const Tensor intVirial_slow,
const Vector extForce_normal,
const Vector extForce_nbond,
const Vector extForce_slow 
) [protected]

Definition at line 1582 of file Controller.C.

References SimParameters::accelMDOn, AVGXY, controlPressure, controlPressure_nbond, controlPressure_normal, controlPressure_slow, Tensor::diagonal(), SimParameters::getCurrentLambda(), Molecule::getVirialTailCorr(), groupPressure, groupPressure_nbond, groupPressure_normal, groupPressure_slow, Tensor::identity(), SimParameters::LJcorrection, Node::molecule, NAMD_bug(), nbondFreq, Node::Object(), Lattice::origin(), outer(), pressure, pressure_amd, pressure_nbond, pressure_normal, pressure_slow, Node::simParameters, slowFreq, state, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, SimParameters::useGroupPressure, virial_amd, and Lattice::volume().

Referenced by multigratorPressure(), and receivePressure().

01585                                                                                             {
01586 
01587   Tensor virial_normal = virial_normal_in;
01588   Tensor virial_nbond = virial_nbond_in;
01589   Tensor virial_slow = virial_slow_in;
01590 
01591   // Tensor tmp = virial_normal;
01592   // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
01593   //   tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
01594 
01595   Node *node = Node::Object();
01596   Molecule *molecule = node->molecule;
01597   SimParameters *simParameters = node->simParameters;
01598   Lattice &lattice = state->lattice;
01599 
01600   BigReal volume;
01601 
01602   Vector extPosition = lattice.origin();
01603   virial_normal -= outer(extForce_normal,extPosition);
01604   virial_nbond -= outer(extForce_nbond,extPosition);
01605   virial_slow -= outer(extForce_slow,extPosition);
01606 
01607   if ( (volume=lattice.volume()) != 0. )
01608   {
01609 
01610     if (simParameters->LJcorrection && volume) {
01611 #ifdef MEM_OPT_VERSION
01612       NAMD_bug("LJcorrection not supported in memory optimized build.");
01613 #else
01614       // Apply tail correction to pressure
01615       BigReal alchLambda = simParameters->getCurrentLambda(step);
01616       virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
01617 #endif
01618     }
01619 
01620     // kinetic energy component included in virials
01621     pressure_normal = virial_normal / volume;
01622     groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01623 
01624     if (simParameters->accelMDOn) {
01625       pressure_amd = virial_amd / volume;
01626       pressure_normal += pressure_amd;
01627       groupPressure_normal +=  pressure_amd;
01628     }
01629 
01630     if ( minimize || ! ( step % nbondFreq ) )
01631     {
01632       pressure_nbond = virial_nbond / volume;
01633       groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01634     }
01635 
01636     if ( minimize || ! ( step % slowFreq ) )
01637     {
01638       pressure_slow = virial_slow / volume;
01639       groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01640     }
01641 
01642     pressure = pressure_normal + pressure_nbond + pressure_slow; 
01643     groupPressure = groupPressure_normal + groupPressure_nbond +
01644           groupPressure_slow;
01645   }
01646   else
01647   {
01648     pressure = Tensor();
01649     groupPressure = Tensor();
01650   }
01651 
01652   if ( simParameters->useGroupPressure )
01653   {
01654     controlPressure_normal = groupPressure_normal;
01655     controlPressure_nbond = groupPressure_nbond;
01656     controlPressure_slow = groupPressure_slow;
01657     controlPressure = groupPressure;
01658   }
01659   else
01660   {
01661     controlPressure_normal = pressure_normal;
01662     controlPressure_nbond = pressure_nbond;
01663     controlPressure_slow = pressure_slow;
01664     controlPressure = pressure;
01665   }
01666 
01667   if ( simParameters->useFlexibleCell ) {
01668     // use symmetric pressure to control rotation
01669     // controlPressure_normal = symmetric(controlPressure_normal);
01670     // controlPressure_nbond = symmetric(controlPressure_nbond);
01671     // controlPressure_slow = symmetric(controlPressure_slow);
01672     // controlPressure = symmetric(controlPressure);
01673     // only use on-diagonal components for now
01674     controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01675     controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01676     controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01677     controlPressure = Tensor::diagonal(diagonal(controlPressure));
01678     if ( simParameters->useConstantRatio ) {
01679 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01680    T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01681       AVGXY(controlPressure_normal);
01682       AVGXY(controlPressure_nbond);
01683       AVGXY(controlPressure_slow);
01684       AVGXY(controlPressure);
01685 #undef AVGXY
01686     }
01687   } else {
01688     controlPressure_normal =
01689   Tensor::identity(trace(controlPressure_normal)/3.);
01690     controlPressure_nbond =
01691   Tensor::identity(trace(controlPressure_nbond)/3.);
01692     controlPressure_slow =
01693   Tensor::identity(trace(controlPressure_slow)/3.);
01694     controlPressure =
01695   Tensor::identity(trace(controlPressure)/3.);
01696   }
01697 }

void Controller::compareChecksums ( int  step,
int  forgiving = 0 
) [protected]

Definition at line 2755 of file Controller.C.

References computeChecksum, endi(), SimParameters::fullElectFrequency, SimParameters::goGroPair, iERROR(), iINFO(), iout, RequireReduction::item(), iWARN(), marginViolations, Node::molecule, SimParameters::mollyOn, NAMD_bug(), Molecule::numAtoms, Molecule::numCalcAngles, Molecule::numCalcAnisos, Molecule::numCalcBonds, Molecule::numCalcCrossterms, Molecule::numCalcDihedrals, Molecule::numCalcExclusions, Molecule::numCalcFullExclusions, Molecule::numCalcImpropers, Molecule::numCalcTholes, Node::Object(), SimParameters::outputPairlists, pairlistWarnings, SimParameters::qmForcesOn, reduction, REDUCTION_ANGLE_CHECKSUM, REDUCTION_ANISO_CHECKSUM, REDUCTION_ATOM_CHECKSUM, REDUCTION_BOND_CHECKSUM, REDUCTION_COMPUTE_CHECKSUM, REDUCTION_CROSSTERM_CHECKSUM, REDUCTION_DIHEDRAL_CHECKSUM, REDUCTION_EXCLUSION_CHECKSUM, REDUCTION_EXCLUSION_CHECKSUM_CUDA, REDUCTION_IMPROPER_CHECKSUM, REDUCTION_MARGIN_VIOLATIONS, REDUCTION_PAIRLIST_WARNINGS, REDUCTION_STRAY_CHARGE_ERRORS, REDUCTION_THOLE_CHECKSUM, simParams, and slowFreq.

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

02755                                                          {
02756     Node *node = Node::Object();
02757     Molecule *molecule = node->molecule;
02758 
02759     // Some basic correctness checking
02760     BigReal checksum, checksum_b;
02761     char errmsg[256];
02762 
02763     checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
02764     if ( ((int)checksum) != molecule->numAtoms ) {
02765       sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
02766               (int)checksum, molecule->numAtoms);
02767       NAMD_bug(errmsg);
02768     }
02769 
02770     checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
02771     if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
02772       if ( ! computeChecksum ) {
02773         computesPartitioned = 0;
02774       } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
02775         computesPartitioned = 1;
02776       } else {
02777         NAMD_bug("Bad global compute count!\n");
02778       }
02779       computeChecksum = ((int)checksum);
02780     }
02781 
02782     checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
02783     if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
02784       sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
02785               (int)checksum, molecule->numCalcBonds);
02786       if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
02787         iout << iWARN << errmsg << endi;
02788       else NAMD_bug(errmsg);
02789     }
02790 
02791     checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
02792     if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
02793       sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
02794               (int)checksum, molecule->numCalcAngles);
02795       if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
02796         iout << iWARN << errmsg << endi;
02797       else NAMD_bug(errmsg);
02798     }
02799 
02800     checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
02801     if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
02802       sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
02803               (int)checksum, molecule->numCalcDihedrals);
02804       if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
02805         iout << iWARN << errmsg << endi;
02806       else NAMD_bug(errmsg);
02807     }
02808 
02809     checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
02810     if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
02811       sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
02812               (int)checksum, molecule->numCalcImpropers);
02813       if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
02814         iout << iWARN << errmsg << endi;
02815       else NAMD_bug(errmsg);
02816     }
02817 
02818     checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
02819     if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
02820       sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
02821               (int)checksum, molecule->numCalcTholes);
02822       if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
02823         iout << iWARN << errmsg << endi;
02824       else NAMD_bug(errmsg);
02825     }
02826 
02827     checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
02828     if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
02829       sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
02830               (int)checksum, molecule->numCalcAnisos);
02831       if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
02832         iout << iWARN << errmsg << endi;
02833       else NAMD_bug(errmsg);
02834     }
02835 
02836     checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
02837     if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
02838       sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
02839               (int)checksum, molecule->numCalcCrossterms);
02840       if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
02841         iout << iWARN << errmsg << endi;
02842       else NAMD_bug(errmsg);
02843     }
02844 
02845     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
02846     if ( ((int64)checksum) > molecule->numCalcExclusions &&
02847          ( ! simParams->mollyOn || step % slowFreq ) ) {
02848       if ( forgiving )
02849         iout << iWARN << "High global exclusion count ("
02850                       << ((int64)checksum) << " vs "
02851                       << molecule->numCalcExclusions << "), possible error!\n"
02852           << iWARN << "This warning is not unusual during minimization.\n"
02853           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02854       else {
02855         char errmsg[256];
02856         sprintf(errmsg, "High global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02857                 (long long)checksum, (long long)molecule->numCalcExclusions);
02858         NAMD_bug(errmsg);
02859       }
02860     }
02861 
02862     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
02863         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02864       if ( forgiving || ! simParams->fullElectFrequency ) {
02865         iout << iWARN << "Low global exclusion count!  ("
02866           << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
02867         if ( forgiving ) {
02868           iout << "\n"
02869             << iWARN << "This warning is not unusual during minimization.\n"
02870             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02871         } else {
02872           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02873         }
02874         iout << endi;
02875       } else {
02876         char errmsg[256];
02877         sprintf(errmsg, "Low global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
02878                 (long long)checksum, (long long)molecule->numCalcExclusions);
02879         NAMD_bug(errmsg);
02880       }
02881     }
02882 
02883 #ifdef NAMD_CUDA
02884     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
02885     if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
02886          ( ! simParams->mollyOn || step % slowFreq ) ) {
02887       if ( forgiving )
02888         iout << iWARN << "High global CUDA exclusion count ("
02889                       << ((int64)checksum) << " vs "
02890                       << molecule->numCalcFullExclusions << "), possible error!\n"
02891           << iWARN << "This warning is not unusual during minimization.\n"
02892           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02893       else {
02894         char errmsg[256];
02895         sprintf(errmsg, "High global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02896                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02897         NAMD_bug(errmsg);
02898       }
02899     }
02900 
02901     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
02902         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02903       if ( forgiving || ! simParams->fullElectFrequency ) {
02904         iout << iWARN << "Low global CUDA exclusion count!  ("
02905           << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
02906         if ( forgiving ) {
02907           iout << "\n"
02908             << iWARN << "This warning is not unusual during minimization.\n"
02909             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02910         } else {
02911           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02912         }
02913         iout << endi;
02914       } else {
02915         char errmsg[256];
02916         sprintf(errmsg, "Low global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
02917                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02918         NAMD_bug(errmsg);
02919       }
02920     }
02921 #endif
02922 
02923     checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
02924     if ( ((int)checksum) && ! marginViolations ) {
02925       iout << iERROR << "Margin is too small for " << ((int)checksum) <<
02926         " atoms during timestep " << step << ".\n" << iERROR <<
02927       "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
02928     }
02929     marginViolations += (int)checksum;
02930 
02931     checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
02932     if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
02933       iout << iINFO <<
02934         "Pairlistdist is too small for " << ((int)checksum) <<
02935         " computes during timestep " << step << ".\n" << endi;
02936     }
02937     if ( simParams->outputPairlists )  pairlistWarnings += (int)checksum;
02938 
02939     checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
02940     if ( checksum ) {
02941       if ( forgiving )
02942         iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
02943       else NAMD_bug("Stray PME grid charges detected!\n");
02944     }
02945 }

BigReal Controller::computeAlchWork ( const int  step  )  [protected]

Definition at line 3858 of file Controller.C.

References bondedEnergy_ti_1, bondedEnergy_ti_2, electEnergy_ti_1, electEnergy_ti_2, electEnergyPME_ti_1, electEnergyPME_ti_2, electEnergySlow_ti_1, electEnergySlow_ti_2, SimParameters::getBondLambda(), SimParameters::getCurrentLambda(), SimParameters::getElecLambda(), SimParameters::getVdwLambda(), ljEnergy_ti_1, ljEnergy_ti_2, and simParams.

Referenced by printEnergies().

03858                                                   {
03859   // alchemical scaling factors for groups 1/2 at the previous lambda
03860   const BigReal oldLambda = simParams->getCurrentLambda(step-1);
03861   const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
03862   const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
03863   const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
03864   const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
03865   const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
03866   const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
03867   // alchemical scaling factors for groups 1/2 at the new lambda
03868   const BigReal alchLambda = simParams->getCurrentLambda(step);
03869   const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
03870   const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
03871   const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
03872   const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
03873   const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
03874   const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda); 
03875 
03876   return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
03877           (elec_lambda_12 - elec_lambda_1)*
03878           (electEnergy_ti_1 + electEnergySlow_ti_1 +  electEnergyPME_ti_1) +
03879           (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
03880           (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
03881           (elec_lambda_22 - elec_lambda_2)*
03882           (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) + 
03883           (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
03884   );
03885 }

void Controller::correctMomentum ( int  step  )  [protected]

Definition at line 1255 of file Controller.C.

References broadcast, endi(), iERROR(), iout, RequireReduction::item(), Vector::length2(), ControllerBroadcasts::momentumCorrection, NAMD_die(), SimpleBroadcastObject< T >::publish(), reduction, REDUCTION_MOMENTUM_MASS, slowFreq, Vector::x, Vector::y, and Vector::z.

Referenced by integrate().

01255                                          {
01256 
01257     Vector momentum;
01258     momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
01259     momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
01260     momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
01261     const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
01262 
01263     if ( momentum.length2() == 0. )
01264       iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
01265     if ( mass == 0. )
01266       NAMD_die("Total mass is zero in Controller::correctMomentum");
01267 
01268     momentum *= (-1./mass);
01269 
01270     broadcast->momentumCorrection.publish(step+slowFreq,momentum);
01271 }

void Controller::cycleBarrier ( int  doBarrier,
int  step 
) [protected]

Definition at line 4119 of file Controller.C.

References broadcast.

Referenced by integrate().

04119                                                      {
04120 #if USE_BARRIER
04121         if (doBarrier) {
04122           broadcast->cycleBarrier.publish(step,1);
04123           CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
04124                   CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04125         }
04126 #endif
04127 }

void Controller::enqueueCollections ( int  timestep  )  [protected]
void Controller::integrate ( int  scriptTask  )  [protected]

Definition at line 429 of file Controller.C.

References adaptTempUpdate(), berendsenPressure(), correctMomentum(), cycleBarrier(), enqueueCollections(), SimParameters::firstTimestep, SimParameters::fullElectFrequency, langevinPiston1(), langevinPiston2(), multigratorPressure(), multigratorTemperature(), SimParameters::N, nbondFreq, SimParameters::nonbondedFrequency, SimParameters::numTraceSteps, Node::Object(), outputExtendedSystem(), outputFepEnergy(), outputTiEnergy(), printDynamicsEnergies(), printFepMessage(), printTiMessage(), reassignVelocities(), rebalanceLoad(), receivePressure(), rescaleaccelMD(), rescaleVelocities(), SCRIPT_RUN, simParams, slowFreq, SimParameters::statsOn, SimParameters::stepsPerCycle, stochRescaleVelocities(), tcoupleVelocities(), traceBarrier(), SimParameters::traceStartStep, and SimParameters::zeroMomentum.

Referenced by algorithm().

00429                                          {
00430     char traceNote[24];
00431   
00432     int step = simParams->firstTimestep;
00433 
00434     const int numberOfSteps = simParams->N;
00435     const int stepsPerCycle = simParams->stepsPerCycle;
00436 
00437     const int zeroMomentum = simParams->zeroMomentum;
00438 
00439     nbondFreq = simParams->nonbondedFrequency;
00440     const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00441     if (dofull)
00442       slowFreq = simParams->fullElectFrequency;
00443     else
00444       slowFreq = simParams->nonbondedFrequency;
00445     if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
00446 
00447   if ( scriptTask == SCRIPT_RUN ) {
00448 
00449     reassignVelocities(step);  // only for full-step velecities
00450     rescaleaccelMD(step);
00451     adaptTempUpdate(step); // Init adaptive tempering;
00452 
00453     receivePressure(step);
00454     if ( zeroMomentum && dofull && ! (step % slowFreq) )
00455                                                 correctMomentum(step);
00456     printFepMessage(step);
00457     printTiMessage(step);
00458     printDynamicsEnergies(step);
00459     outputFepEnergy(step);
00460     outputTiEnergy(step);
00461     if(traceIsOn()){
00462         traceUserEvent(eventEndOfTimeStep);
00463         sprintf(traceNote, "s:%d", step);
00464         traceUserSuppliedNote(traceNote);
00465     }
00466     outputExtendedSystem(step);
00467     rebalanceLoad(step);
00468 
00469   }
00470 
00471     // Handling SIGINT doesn't seem to be working on Lemieux, and it
00472     // sometimes causes the net-xxx versions of NAMD to segfault on exit, 
00473     // so disable it for now.
00474     // namd_sighandler_t oldhandler = signal(SIGINT, 
00475     //  (namd_sighandler_t)my_sigint_handler);
00476     for ( ++step ; step <= numberOfSteps; ++step )
00477     {
00478 
00479         adaptTempUpdate(step);
00480         rescaleVelocities(step);
00481         tcoupleVelocities(step);
00482         stochRescaleVelocities(step);
00483         berendsenPressure(step);
00484         langevinPiston1(step);
00485         rescaleaccelMD(step);
00486         enqueueCollections(step);  // after lattice scaling!
00487         receivePressure(step);
00488         if ( zeroMomentum && dofull && ! (step % slowFreq) )
00489                                                 correctMomentum(step);
00490         langevinPiston2(step);
00491         reassignVelocities(step);
00492 
00493       multigratorTemperature(step, 1);
00494       multigratorPressure(step, 1);
00495       multigratorPressure(step, 2);
00496       multigratorTemperature(step, 2);
00497 
00498         printDynamicsEnergies(step);
00499         outputFepEnergy(step);
00500         outputTiEnergy(step);
00501         if(traceIsOn()){
00502             traceUserEvent(eventEndOfTimeStep);
00503             sprintf(traceNote, "s:%d", step);
00504             traceUserSuppliedNote(traceNote);
00505         }
00506   // if (gotsigint) {
00507   //   iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
00508   //   NAMD_quit();
00509   // }
00510         outputExtendedSystem(step);
00511 #if CYCLE_BARRIER
00512         cycleBarrier(!((step+1) % stepsPerCycle),step);
00513 #elif  PME_BARRIER
00514         cycleBarrier(dofull && !(step%slowFreq),step);
00515 #elif  STEP_BARRIER
00516         cycleBarrier(1, step);
00517 #endif
00518                 
00519         if(Node::Object()->specialTracing || simParams->statsOn){               
00520                  int bstep = simParams->traceStartStep;
00521                  int estep = bstep + simParams->numTraceSteps;          
00522                  if(step == bstep){
00523                          traceBarrier(1, step);
00524                  }
00525                  if(step == estep){                     
00526                          traceBarrier(0, step);
00527                  }
00528          }
00529 
00530 #ifdef MEASURE_NAMD_WITH_PAPI
00531         if(simParams->papiMeasure) {
00532                 int bstep = simParams->papiMeasureStartStep;
00533                 int estep = bstep + simParams->numPapiMeasureSteps;
00534                 if(step == bstep) {
00535                         papiMeasureBarrier(1, step);
00536                 }
00537                 if(step == estep) {
00538                         papiMeasureBarrier(0, step);
00539                 }
00540         }
00541 #endif
00542          
00543         rebalanceLoad(step);
00544 
00545 #if  PME_BARRIER
00546         cycleBarrier(dofull && !((step+1)%slowFreq),step);   // step before PME
00547 #endif
00548     }
00549     // signal(SIGINT, oldhandler);
00550 }

void Controller::langevinPiston1 ( int  step  )  [protected]

Definition at line 987 of file Controller.C.

References BOLTZMANN, broadcast, Lattice::c(), controlNumDegFreedom, controlPressure, controlPressure_nbond, controlPressure_slow, Tensor::diagonal(), SimParameters::dt, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, Random::gaussian(), Random::gaussian_vector(), Tensor::identity(), iINFO(), iout, ControllerState::langevinPiston_strainRate, SimParameters::langevinPistonBarrier, SimParameters::langevinPistonDecay, SimParameters::langevinPistonOn, SimParameters::langevinPistonPeriod, SimParameters::langevinPistonTarget, SimParameters::langevinPistonTemp, nbondFreq, positionRescaleFactor, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), random, Lattice::rescale(), simParams, slowFreq, state, strainRate_old, SimParameters::surfaceTensionTarget, SimParameters::useConstantArea, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Tensor::xx, Tensor::yy, Vector::z, and Tensor::zz.

Referenced by integrate().

00988 {
00989   if ( simParams->langevinPistonOn )
00990   {
00991     Tensor &strainRate = langevinPiston_strainRate;
00992     int cellDims = simParams->useFlexibleCell ? 1 : 3;
00993     BigReal dt = simParams->dt;
00994     BigReal dt_long = slowFreq * dt;
00995     BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
00996     BigReal tau = simParams->langevinPistonPeriod;
00997     BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00998 
00999 #ifdef DEBUG_PRESSURE
01000     iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
01001 #endif
01002 
01003     if ( ! ( (step-1) % slowFreq ) )
01004     {
01005       BigReal gamma = 1 / simParams->langevinPistonDecay;
01006       BigReal f1 = exp( -0.5 * dt_long * gamma );
01007       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
01008       strainRate *= f1;
01009       if ( simParams->useFlexibleCell ) {
01010         // We only use on-diagonal terms (for now)
01011         if ( simParams->useConstantRatio ) {
01012           BigReal r = f2 * random->gaussian();
01013           strainRate.xx += r;
01014           strainRate.yy += r;
01015           strainRate.zz += f2 * random->gaussian();
01016         } else {
01017           strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
01018         }
01019       } else {
01020         strainRate += f2 * Tensor::identity(random->gaussian());
01021       }
01022 
01023 #ifdef DEBUG_PRESSURE
01024       iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
01025 #endif
01026     }
01027 
01028     // Apply surface tension.  If surfaceTensionTarget is zero, we get
01029     // the default (isotropic pressure) case.
01030     
01031     Tensor ptarget;
01032     ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
01033     if ( simParams->surfaceTensionTarget != 0. ) {
01034       ptarget.xx = ptarget.yy = simParams->langevinPistonTarget - 
01035         simParams->surfaceTensionTarget / state->lattice.c().z;
01036     }
01037 
01038     strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01039       ( controlPressure - ptarget );
01040 
01041     if (simParams->fixCellDims) {
01042       if (simParams->fixCellDimX) strainRate.xx = 0;
01043       if (simParams->fixCellDimY) strainRate.yy = 0;
01044       if (simParams->fixCellDimZ) strainRate.zz = 0;
01045     }
01046 
01047 
01048 #ifdef DEBUG_PRESSURE
01049     iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
01050 #endif
01051 
01052     if ( simParams->langevinPistonBarrier ) {
01053     if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
01054     {
01055       // We only use on-diagonal terms (for now)
01056       Tensor factor;
01057       if ( !simParams->useConstantArea ) {
01058         factor.xx = exp( dt_long * strainRate.xx );
01059         factor.yy = exp( dt_long * strainRate.yy );
01060       } else {
01061         factor.xx = factor.yy = 1;
01062       }
01063       factor.zz = exp( dt_long * strainRate.zz );
01064       broadcast->positionRescaleFactor.publish(step,factor);
01065       state->lattice.rescale(factor);
01066 #ifdef DEBUG_PRESSURE
01067       iout << iINFO << "rescaling by: " << factor << "\n";
01068 #endif
01069     }
01070     } else { // langevinPistonBarrier
01071     if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
01072     {
01073       if ( ! positionRescaleFactor.xx ) {  // first pass without MTS
01074       // We only use on-diagonal terms (for now)
01075       Tensor factor;
01076       if ( !simParams->useConstantArea ) {
01077         factor.xx = exp( dt_long * strainRate.xx );
01078         factor.yy = exp( dt_long * strainRate.yy );
01079       } else {
01080         factor.xx = factor.yy = 1;
01081       }
01082       factor.zz = exp( dt_long * strainRate.zz );
01083       broadcast->positionRescaleFactor.publish(step,factor);
01084       positionRescaleFactor = factor;
01085       strainRate_old = strainRate;
01086       }
01087       state->lattice.rescale(positionRescaleFactor);
01088 #ifdef DEBUG_PRESSURE
01089       iout << iINFO << "rescaling by: " << positionRescaleFactor << "\n";
01090 #endif
01091     }
01092     if ( ! ( (step-slowFreq/2) % slowFreq ) )
01093     {
01094       // We only use on-diagonal terms (for now)
01095       Tensor factor;
01096       if ( !simParams->useConstantArea ) {
01097         factor.xx = exp( (dt+dt_long) * strainRate.xx - dt * strainRate_old.xx );
01098         factor.yy = exp( (dt+dt_long) * strainRate.yy - dt * strainRate_old.yy );
01099       } else {
01100         factor.xx = factor.yy = 1;
01101       }
01102       factor.zz = exp( (dt+dt_long) * strainRate.zz - dt * strainRate_old.zz );
01103       broadcast->positionRescaleFactor.publish(step+1,factor);
01104       positionRescaleFactor = factor;
01105       strainRate_old = strainRate;
01106     }
01107     }
01108 
01109     // corrections to integrator
01110     if ( ! ( step % nbondFreq ) )
01111     {
01112 #ifdef DEBUG_PRESSURE
01113       iout << iINFO << "correcting strain rate for nbond, ";
01114 #endif
01115       strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01116                 ( (nbondFreq - 1) * controlPressure_nbond );
01117 #ifdef DEBUG_PRESSURE
01118       iout << "strain rate: " << strainRate << "\n";
01119 #endif
01120     }
01121     if ( ! ( step % slowFreq ) )
01122     {
01123 #ifdef DEBUG_PRESSURE
01124       iout << iINFO << "correcting strain rate for slow, ";
01125 #endif
01126       strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01127                 ( (slowFreq - 1) * controlPressure_slow );
01128 #ifdef DEBUG_PRESSURE
01129       iout << "strain rate: " << strainRate << "\n";
01130 #endif
01131     }
01132     if (simParams->fixCellDims) {
01133       if (simParams->fixCellDimX) strainRate.xx = 0;
01134       if (simParams->fixCellDimY) strainRate.yy = 0;
01135       if (simParams->fixCellDimZ) strainRate.zz = 0;
01136     }
01137 
01138   }
01139 
01140 }

void Controller::langevinPiston2 ( int  step  )  [protected]

Definition at line 1142 of file Controller.C.

References BOLTZMANN, Lattice::c(), controlNumDegFreedom, controlPressure, controlPressure_nbond, controlPressure_slow, Tensor::diagonal(), SimParameters::dt, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, Random::gaussian(), Random::gaussian_vector(), Tensor::identity(), iINFO(), iout, ControllerState::langevinPiston_strainRate, SimParameters::langevinPistonDecay, SimParameters::langevinPistonOn, SimParameters::langevinPistonPeriod, SimParameters::langevinPistonTarget, SimParameters::langevinPistonTemp, nbondFreq, random, simParams, slowFreq, state, SimParameters::surfaceTensionTarget, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Tensor::xx, Tensor::yy, Vector::z, and Tensor::zz.

Referenced by integrate().

01143 {
01144   if ( simParams->langevinPistonOn )
01145   {
01146     Tensor &strainRate = langevinPiston_strainRate;
01147     int cellDims = simParams->useFlexibleCell ? 1 : 3;
01148     BigReal dt = simParams->dt;
01149     BigReal dt_long = slowFreq * dt;
01150     BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
01151     BigReal tau = simParams->langevinPistonPeriod;
01152     BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
01153 
01154     // corrections to integrator
01155     if ( ! ( step % nbondFreq ) )
01156     {
01157 #ifdef DEBUG_PRESSURE
01158       iout << iINFO << "correcting strain rate for nbond, ";
01159 #endif
01160       strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01161                 ( (nbondFreq - 1) * controlPressure_nbond );
01162 #ifdef DEBUG_PRESSURE
01163       iout << "strain rate: " << strainRate << "\n";
01164 #endif
01165     }
01166     if ( ! ( step % slowFreq ) )
01167     {
01168 #ifdef DEBUG_PRESSURE
01169       iout << iINFO << "correcting strain rate for slow, ";
01170 #endif
01171       strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01172                 ( (slowFreq - 1) * controlPressure_slow );
01173 #ifdef DEBUG_PRESSURE
01174       iout << "strain rate: " << strainRate << "\n";
01175 #endif
01176     }
01177 
01178     // Apply surface tension.  If surfaceTensionTarget is zero, we get
01179     // the default (isotropic pressure) case.
01180    
01181     Tensor ptarget;
01182     ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
01183     if ( simParams->surfaceTensionTarget != 0. ) {
01184       ptarget.xx = ptarget.yy = simParams->langevinPistonTarget - 
01185         simParams->surfaceTensionTarget / state->lattice.c().z;
01186     }
01187 
01188     strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
01189       ( controlPressure - ptarget );
01190 
01191  
01192 #ifdef DEBUG_PRESSURE
01193     iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
01194 #endif
01195 
01196     if ( ! ( step % slowFreq ) )
01197     {
01198       BigReal gamma = 1 / simParams->langevinPistonDecay;
01199       BigReal f1 = exp( -0.5 * dt_long * gamma );
01200       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
01201       strainRate *= f1;
01202       if ( simParams->useFlexibleCell ) {
01203         // We only use on-diagonal terms (for now)
01204         if ( simParams->useConstantRatio ) {
01205           BigReal r = f2 * random->gaussian();
01206           strainRate.xx += r;
01207           strainRate.yy += r;
01208           strainRate.zz += f2 * random->gaussian();
01209         } else {
01210           strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
01211         }
01212       } else {
01213         strainRate += f2 * Tensor::identity(random->gaussian());
01214       }
01215 #ifdef DEBUG_PRESSURE
01216       iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
01217 #endif
01218     }
01219 
01220 #ifdef DEBUG_PRESSURE
01221     iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
01222 #endif
01223     if (simParams->fixCellDims) {
01224       if (simParams->fixCellDimX) strainRate.xx = 0;
01225       if (simParams->fixCellDimY) strainRate.yy = 0;
01226       if (simParams->fixCellDimZ) strainRate.zz = 0;
01227     }
01228   }
01229 
01230 
01231 }

void Controller::minimize (  )  [protected]

Definition at line 594 of file Controller.C.

References broadcast, CALCULATE, minpoint::dudx, endi(), enqueueCollections(), SimParameters::firstTimestep, iout, RequireReduction::item(), min_energy, min_f_dot_f, min_f_dot_v, min_huge_count, min_reduction, min_v_dot_v, ControllerBroadcasts::minimizeCoefficient, SimParameters::minLineGoal, SimParameters::minVerbose, MOVETO, SimParameters::N, nbondFreq, minpoint::noGradient, PRINT_BRACKET, SimpleBroadcastObject< T >::publish(), RequireReduction::require(), simParams, slowFreq, minpoint::u, and minpoint::x.

Referenced by algorithm().

00594                           {
00595   // iout << "Controller::minimize() called.\n" << endi;
00596 
00597   const int minVerbose = simParams->minVerbose;
00598   const int numberOfSteps = simParams->N;
00599   int step = simParams->firstTimestep;
00600   slowFreq = nbondFreq = 1;
00601   BigReal linegoal = simParams->minLineGoal;  // 1.0e-3
00602   const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
00603 
00604   CALCULATE
00605 
00606   int minSeq = 0;
00607 
00608   // just move downhill until initial bad contacts go away or 100 steps
00609   int old_min_huge_count = 2000000000;
00610   int steps_at_same_min_huge_count = 0;
00611   for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
00612     if ( min_huge_count >= old_min_huge_count ) {
00613       if ( ++steps_at_same_min_huge_count > 10 ) break;
00614     } else {
00615       old_min_huge_count = min_huge_count;
00616       steps_at_same_min_huge_count = 0;
00617     }
00618     iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
00619     broadcast->minimizeCoefficient.publish(minSeq++,1.);
00620     enqueueCollections(step);
00621     CALCULATE
00622   }
00623   if ( min_huge_count ) {
00624     iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
00625   }
00626   iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00627 
00628   int atStart = 2;
00629   int errorFactor = 10;
00630   BigReal old_f_dot_f = min_f_dot_f;
00631   broadcast->minimizeCoefficient.publish(minSeq++,0.);
00632   broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
00633   int newDir = 1;
00634   min_f_dot_v = min_f_dot_f;
00635   min_v_dot_v = min_f_dot_f;
00636   while ( 1 ) {
00637     // line minimization
00638     // bracket minimum on line
00639     minpoint lo,hi,mid,last;
00640     BigReal x = 0;
00641     lo.x = x;
00642     lo.u = min_energy;
00643     lo.dudx = -1. * min_f_dot_v;
00644     lo.noGradient = min_huge_count;
00645     mid = lo;
00646     last = mid;
00647     if ( minVerbose ) {
00648       iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
00649       if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
00650       iout << "\n" << endi;
00651     }
00652     BigReal tol = fabs( linegoal * min_f_dot_v );
00653     iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
00654             fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
00655     int start_with_huge = last.noGradient;
00656     min_reduction->require();
00657     BigReal maxstep = 0.1 / sqrt(min_reduction->item(0));
00658     x = maxstep; MOVETO(x);
00659     // bracket minimum on line
00660     while ( last.u < mid.u ) {
00661       if ( last.dudx < mid.dudx * (5.5 - x/maxstep)/5. ) {
00662         x = 6 * maxstep; break;
00663       }
00664       lo = mid; mid = last;
00665       x += maxstep;
00666       if ( x > 5.5 * maxstep ) break;
00667       MOVETO(x)
00668     }
00669     if ( x > 5.5 * maxstep ) {
00670       iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" << endi;
00671       broadcast->minimizeCoefficient.publish(minSeq++,0.);
00672       broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
00673       newDir = 1;
00674       old_f_dot_f = min_f_dot_f;
00675       min_f_dot_v = min_f_dot_f;
00676       min_v_dot_v = min_f_dot_f;
00677       continue;
00678     }
00679     hi = last;
00680 #define PRINT_BRACKET \
00681     iout << "LINE MINIMIZER BRACKET: DX " \
00682          << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
00683         " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
00684         lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
00685     PRINT_BRACKET
00686     // converge on minimum on line
00687     int itcnt;
00688     for ( itcnt = 10; itcnt > 0; --itcnt ) {
00689       int progress = 1;
00690       // select new position
00691       if ( mid.noGradient ) {
00692        if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) {  // subdivide left side
00693         x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
00694         MOVETO(x)
00695         if ( last.u <= mid.u ) {
00696           hi = mid; mid = last;
00697         } else {
00698           lo = last;
00699         }
00700        } else {  // subdivide right side
00701         x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
00702         MOVETO(x)
00703         if ( last.u <= mid.u ) {
00704           lo = mid; mid = last;
00705         } else {
00706           hi = last;
00707         }
00708        }
00709       } else {
00710        if ( mid.dudx > 0. ) {  // subdivide left side
00711         BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
00712         BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
00713         x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
00714         BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
00715         if ( xdiv ) x /= xdiv; else x = last.x;
00716         if ( x > altxhi ) x = altxhi;
00717         if ( x < altxlo ) x = altxlo;
00718         if ( x-last.x == 0 ) break;
00719         MOVETO(x)
00720         if ( last.u <= mid.u ) {
00721           hi = mid; mid = last;
00722         } else {
00723           if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
00724           lo = last;
00725         }
00726        } else {  // subdivide right side
00727         BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
00728         BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
00729         x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
00730         BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
00731         if ( xdiv ) x /= xdiv; else x = last.x;
00732         if ( x < altxlo ) x = altxlo;
00733         if ( x > altxhi ) x = altxhi;
00734         if ( x-last.x == 0 ) break;
00735         MOVETO(x)
00736         if ( last.u <= mid.u ) {
00737           lo = mid; mid = last;
00738         } else {
00739           if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
00740           hi = last;
00741         }
00742        }
00743       }
00744       PRINT_BRACKET
00745       if ( ! progress ) {
00746         MOVETO(mid.x);
00747         break;
00748       }
00749       if ( fabs(last.dudx) < tol ) break;
00750       if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
00751       if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
00752     }
00753     // new direction
00754     broadcast->minimizeCoefficient.publish(minSeq++,0.);
00755     BigReal c = min_f_dot_f / old_f_dot_f;
00756     c = ( c > 1.5 ? 1.5 : c );
00757     if ( atStart ) { c = 0; --atStart; }
00758     if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
00759       c = 0;
00760       if ( errorFactor < 100 ) errorFactor += 10;
00761     }
00762     if ( c == 0 ) {
00763       iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00764     }
00765     broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
00766     newDir = 1;
00767     old_f_dot_f = min_f_dot_f;
00768     min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
00769     min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
00770   }
00771 
00772 }

BigReal Controller::multigatorCalcEnthalpy ( BigReal  potentialEnergy,
int  step,
int  minimize 
) [protected]

Definition at line 919 of file Controller.C.

References BOLTZMANN, controlNumDegFreedom, kineticEnergy, Node::molecule, SimParameters::multigratorNoseHooverChainLength, multigratorNu, multigratorOmega, SimParameters::multigratorPressureRelaxationTime, SimParameters::multigratorPressureTarget, SimParameters::multigratorTemperatureTarget, multigratorXi, multigratorZeta, numDegFreedom, Node::Object(), Node::simParameters, simParams, state, and Lattice::volume().

Referenced by printEnergies().

00919                                                                                           {
00920   Node *node = Node::Object();
00921   Molecule *molecule = node->molecule;
00922   SimParameters *simParameters = node->simParameters;
00923 
00924   BigReal V = state->lattice.volume();
00925   BigReal P0 = simParams->multigratorPressureTarget;
00926   BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00927   BigReal NG = controlNumDegFreedom;
00928   BigReal Nf = numDegFreedom;
00929   BigReal tauV = simParams->multigratorPressureRelaxationTime;
00930   BigReal sumZeta = 0.0;
00931   for (int i=1;i < simParams->multigratorNoseHooverChainLength;i++) {
00932     sumZeta += multigratorZeta[i];
00933   }
00934   BigReal nuOmegaSum = 0.0;
00935   for (int i=0;i < simParams->multigratorNoseHooverChainLength;i++) {
00936     nuOmegaSum += multigratorNu[i]*multigratorNu[i]/(2.0*multigratorOmega[i]);
00937   }
00938   BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
00939   BigReal eta = sqrt(kT0*W)*multigratorXi;
00940 
00941   BigReal enthalpy = kineticEnergy + potentialEnergy + eta*eta/(2.0*W) + P0*V + nuOmegaSum + kT0*(Nf*multigratorZeta[0] + sumZeta);
00942 
00943 //  if (!(step % 100))
00944     // fprintf(stderr, "enthalpy %lf %lf %lf %lf %lf %lf %lf\n", enthalpy,
00945     //   kineticEnergy, potentialEnergy, eta*eta/(2.0*W), P0*V, nuOmegaSum, kT0*(Nf*multigratorZeta[0] + sumZeta));
00946 
00947   return enthalpy;
00948 }

void Controller::multigratorPressure ( int  step,
int  callNumber 
) [protected]

Definition at line 778 of file Controller.C.

References BOLTZMANN, broadcast, calcPressure(), controlNumDegFreedom, controlPressure, SimParameters::dt, GET_TENSOR, GET_VECTOR, Tensor::identity(), RequireReduction::item(), kineticEnergy, momentumSqrSum, SimParameters::multigratorOn, SimParameters::multigratorPressureFreq, SimParameters::multigratorPressureRelaxationTime, SimParameters::multigratorPressureTarget, SimParameters::multigratorTemperatureTarget, multigratorXi, multigratorXiT, numDegFreedom, ControllerBroadcasts::positionRescaleFactor, ControllerBroadcasts::positionRescaleFactor2, SimpleBroadcastObject< T >::publish(), reduction, REDUCTION_CENTERED_KINETIC_ENERGY, RequireReduction::require(), Lattice::rescale(), simParams, state, temperature, ControllerBroadcasts::velocityRescaleTensor, ControllerBroadcasts::velocityRescaleTensor2, and Lattice::volume().

Referenced by integrate().

00778                                                              {
00779   if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
00780     BigReal P0 = simParams->multigratorPressureTarget;
00781     BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00782     BigReal NG = controlNumDegFreedom;
00783     BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
00784     BigReal tau = simParams->multigratorPressureRelaxationTime;
00785     BigReal s = 0.5*simParams->multigratorPressureFreq*simParams->dt/tau;
00786     {
00787       // Compute new scaling factors and send them to Sequencer
00788       BigReal V = state->lattice.volume();
00789       BigReal Pinst = trace(controlPressure)/3.0;
00790       BigReal PGsum = trace(momentumSqrSum);
00791       //
00792       multigratorXiT = multigratorXi + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
00793       BigReal scale = exp(s*NGfac*multigratorXiT);
00794       BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
00795       // fprintf(stderr, "%d | T %lf P %lf V %1.3lf\n", step, temperature, Pinst*PRESSUREFACTOR, V);
00796       Tensor scaleTensor = Tensor::identity(scale);
00797       Tensor volScaleTensor = Tensor::identity(scale);
00798       Tensor velScaleTensor = Tensor::identity(velScale);
00799       state->lattice.rescale(volScaleTensor);
00800       if (callNumber == 1) {
00801         broadcast->positionRescaleFactor.publish(step,scaleTensor);
00802         broadcast->velocityRescaleTensor.publish(step,velScaleTensor);
00803       } else {
00804         broadcast->positionRescaleFactor2.publish(step,scaleTensor);
00805         broadcast->velocityRescaleTensor2.publish(step,velScaleTensor);      
00806       }
00807     }
00808 
00809     {
00810       // Wait here for Sequencer to finish scaling and force computation
00811       reduction->require();
00812       Tensor virial_normal;
00813       Tensor virial_nbond;
00814       Tensor virial_slow;
00815       Tensor intVirial_normal;
00816       Tensor intVirial_nbond;
00817       Tensor intVirial_slow;
00818       Vector extForce_normal;
00819       Vector extForce_nbond;
00820       Vector extForce_slow;
00821       GET_TENSOR(momentumSqrSum, reduction, REDUCTION_MOMENTUM_SQUARED);
00822       GET_TENSOR(virial_normal, reduction, REDUCTION_VIRIAL_NORMAL);
00823       GET_TENSOR(virial_nbond, reduction, REDUCTION_VIRIAL_NBOND);
00824       GET_TENSOR(virial_slow, reduction, REDUCTION_VIRIAL_SLOW);
00825       GET_TENSOR(intVirial_normal, reduction, REDUCTION_INT_VIRIAL_NORMAL);
00826       GET_TENSOR(intVirial_nbond, reduction, REDUCTION_INT_VIRIAL_NBOND);
00827       GET_TENSOR(intVirial_slow, reduction, REDUCTION_INT_VIRIAL_SLOW);
00828       GET_VECTOR(extForce_normal, reduction, REDUCTION_EXT_FORCE_NORMAL);
00829       GET_VECTOR(extForce_nbond, reduction, REDUCTION_EXT_FORCE_NBOND);
00830       GET_VECTOR(extForce_slow, reduction, REDUCTION_EXT_FORCE_SLOW);
00831       calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
00832         intVirial_normal, intVirial_nbond, intVirial_slow,
00833         extForce_normal, extForce_nbond, extForce_slow);
00834       if (callNumber == 2) {
00835         // Update temperature for the Temperature Cycle that is coming next
00836         BigReal kineticEnergy = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
00837         temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
00838       }
00839     }
00840 
00841     {
00842       // Update pressure integrator
00843       BigReal V = state->lattice.volume();
00844       BigReal Pinst = trace(controlPressure)/3.0;
00845       BigReal PGsum = trace(momentumSqrSum);
00846       //
00847       multigratorXi = multigratorXiT + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
00848     }
00849 
00850   }
00851 }

void Controller::multigratorTemperature ( int  step,
int  callNumber 
) [protected]

Definition at line 853 of file Controller.C.

References BOLTZMANN, broadcast, SimParameters::dt, GET_TENSOR, RequireReduction::item(), kineticEnergy, momentumSqrSum, MULTIGRATOR_REDUCTION_KINETIC_ENERGY, SimParameters::multigratorNoseHooverChainLength, multigratorNu, multigratorNuT, multigratorOmega, SimParameters::multigratorOn, SimParameters::multigratorPressureFreq, multigratorReduction, SimParameters::multigratorTemperatureFreq, SimParameters::multigratorTemperatureRelaxationTime, SimParameters::multigratorTemperatureTarget, multigratorZeta, numDegFreedom, SimpleBroadcastObject< T >::publish(), RequireReduction::require(), simParams, temperature, ControllerBroadcasts::velocityRescaleFactor, and ControllerBroadcasts::velocityRescaleFactor2.

Referenced by integrate().

00853                                                                 {
00854   if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
00855     BigReal tau = simParams->multigratorTemperatureRelaxationTime;
00856     BigReal t = 0.5*simParams->multigratorTemperatureFreq*simParams->dt;
00857     BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00858     BigReal Nf = numDegFreedom;
00859     int n = simParams->multigratorNoseHooverChainLength;
00860     BigReal T1, T2, v;
00861     {
00862       T1 = temperature;
00863       BigReal kTinst = BOLTZMANN * temperature;
00864       for (int i=n-1;i >= 0;i--) {
00865         if (i == 0) {
00866           BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
00867           multigratorNuT[0] = exp(-0.5*t*NuOmega)*multigratorNu[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
00868         } else if (i == n-1) {
00869           multigratorNuT[n-1] = multigratorNu[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
00870         } else {
00871           BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
00872           multigratorNuT[i] = exp(-0.5*t*NuOmega)*multigratorNu[i] + 
00873           0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
00874         }
00875       }
00876       BigReal velScale = exp(-t*multigratorNuT[0]/multigratorOmega[0]);
00877       v = velScale;
00878       if (callNumber == 1)
00879         broadcast->velocityRescaleFactor.publish(step,velScale);
00880       else
00881         broadcast->velocityRescaleFactor2.publish(step,velScale);
00882     }
00883 
00884     {
00885       // Wait here for Sequencer to finish scaling and re-calculating kinetic energy
00886       multigratorReduction->require();
00887       BigReal kineticEnergy = multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY);
00888       temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
00889       T2 = temperature;
00890       if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
00891         // If this is pressure cycle, receive new momentum product
00892         GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
00893       }
00894     }
00895 
00896     // fprintf(stderr, "%d | T %lf scale %lf T' %lf\n", step, T1, v, T2);
00897 
00898     {
00899       BigReal kTinst = BOLTZMANN * temperature;
00900       for (int i=0;i < n;i++) {
00901         if (i == 0) {
00902           BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
00903           multigratorNu[0] = exp(-0.5*t*NuOmega)*multigratorNuT[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
00904         } else if (i == n-1) {
00905           multigratorNu[n-1] = multigratorNuT[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
00906         } else {
00907           BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
00908           multigratorNu[i] = exp(-0.5*t*NuOmega)*multigratorNuT[i] + 
00909           0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
00910         }
00911         multigratorZeta[i] += t*multigratorNuT[i]/multigratorOmega[i];
00912       }
00913     }
00914 
00915   }
00916 }

void Controller::outputExtendedSystem ( int  step  )  [protected]

Definition at line 3947 of file Controller.C.

References ofstream_namd::clear(), ofstream_namd::close(), END_OF_RUN, endi(), FILE_OUTPUT, SimParameters::firstTimestep, ofstream_namd::flush(), iout, ofstream_namd::is_open(), SimParameters::N, NAMD_backup_file(), NAMD_err(), ofstream_namd::open(), SimParameters::outputFilename, SimParameters::restartFilename, SimParameters::restartFrequency, SimParameters::restartSave, simParams, writeExtendedSystemData(), writeExtendedSystemLabels(), xstFile, SimParameters::xstFilename, and SimParameters::xstFrequency.

Referenced by algorithm(), and integrate().

03948 {
03949 
03950   if ( step >= 0 ) {
03951 
03952     // Write out eXtended System Trajectory (XST) file
03953     if ( simParams->xstFrequency &&
03954          ((step % simParams->xstFrequency) == 0) )
03955     {
03956       if ( ! xstFile.is_open() )
03957       {
03958         iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
03959         NAMD_backup_file(simParams->xstFilename);
03960         xstFile.open(simParams->xstFilename);
03961         while (!xstFile) {
03962           if ( errno == EINTR ) {
03963             CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
03964             xstFile.clear();
03965             xstFile.open(simParams->xstFilename);
03966             continue;
03967           }
03968           char err_msg[257];
03969           sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
03970           NAMD_err(err_msg);
03971         }
03972         xstFile << "# NAMD extended system trajectory file" << std::endl;
03973         writeExtendedSystemLabels(xstFile);
03974       }
03975       writeExtendedSystemData(step,xstFile);
03976       xstFile.flush();
03977     }
03978 
03979     // Write out eXtended System Configuration (XSC) files
03980     //  Output a restart file
03981     if ( simParams->restartFrequency &&
03982          ((step % simParams->restartFrequency) == 0) &&
03983          (step != simParams->firstTimestep) )
03984     {
03985       iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
03986                 << step << "\n" << endi;
03987       char fname[140];
03988       const char *bsuffix = ".old";
03989       strcpy(fname, simParams->restartFilename);
03990       if ( simParams->restartSave ) {
03991         char timestepstr[20];
03992         sprintf(timestepstr,".%d",step);
03993         strcat(fname, timestepstr);
03994         bsuffix = ".BAK";
03995       }
03996       strcat(fname, ".xsc");
03997       NAMD_backup_file(fname,bsuffix);
03998       ofstream_namd xscFile(fname);
03999       while (!xscFile) {
04000         if ( errno == EINTR ) {
04001           CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
04002           xscFile.clear();
04003           xscFile.open(fname);
04004           continue;
04005         }
04006         char err_msg[257];
04007         sprintf(err_msg, "Error opening XSC restart file %s",fname);
04008         NAMD_err(err_msg);
04009       } 
04010       xscFile << "# NAMD extended system configuration restart file" << std::endl;
04011       writeExtendedSystemLabels(xscFile);
04012       writeExtendedSystemData(step,xscFile);
04013       if (!xscFile) {
04014         char err_msg[257];
04015         sprintf(err_msg, "Error writing XSC restart file %s",fname);
04016         NAMD_err(err_msg);
04017       } 
04018     }
04019 
04020   }
04021 
04022   //  Output final coordinates
04023   if (step == FILE_OUTPUT || step == END_OF_RUN)
04024   {
04025     int realstep = ( step == FILE_OUTPUT ?
04026         simParams->firstTimestep : simParams->N );
04027     iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
04028                 << realstep << "\n" << endi;
04029     static char fname[140];
04030     strcpy(fname, simParams->outputFilename);
04031     strcat(fname, ".xsc");
04032     NAMD_backup_file(fname);
04033     ofstream_namd xscFile(fname);
04034     while (!xscFile) {
04035       if ( errno == EINTR ) {
04036         CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
04037         xscFile.clear();
04038         xscFile.open(fname);
04039         continue;
04040       }
04041       char err_msg[257];
04042       sprintf(err_msg, "Error opening XSC output file %s",fname);
04043       NAMD_err(err_msg);
04044     } 
04045     xscFile << "# NAMD extended system configuration output file" << std::endl;
04046     writeExtendedSystemLabels(xscFile);
04047     writeExtendedSystemData(realstep,xscFile);
04048     if (!xscFile) {
04049       char err_msg[257];
04050       sprintf(err_msg, "Error writing XSC output file %s",fname);
04051       NAMD_err(err_msg);
04052     } 
04053   }
04054 
04055   //  Close trajectory file
04056   if (step == END_OF_RUN) {
04057     if ( xstFile.is_open() ) {
04058       xstFile.close();
04059       iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
04060     }
04061   }
04062 
04063 }

void Controller::outputFepEnergy ( int  step  )  [protected]

Definition at line 3655 of file Controller.C.

References SimParameters::alchEnsembleAvg, SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchIDWSFreq, SimParameters::alchLambda, SimParameters::alchLambda2, SimParameters::alchLambdaIDWS, SimParameters::alchOn, SimParameters::alchOutFile, SimParameters::alchOutFreq, SimParameters::alchTemp, BOLTZMANN, bondedEnergyDiff_f, dE, dG, electEnergy, electEnergy_f, electEnergySlow, electEnergySlow_f, endi(), exp_dE_ByRT, fepFile, FepNo, fepSum, SimParameters::firstTimestep, ofstream_namd::flush(), iout, ofstream_namd::is_open(), ljEnergy, ljEnergy_f, SimParameters::N, NAMD_backup_file(), net_dE, ofstream_namd::open(), simParams, and writeFepEnergyData().

Referenced by integrate().

03655                                          {
03656   if (simParams->alchOn && simParams->alchFepOn) {
03657     const int stepInRun = step - simParams->firstTimestep;
03658     const int alchEquilSteps = simParams->alchEquilSteps;
03659     const BigReal alchLambda = simParams->alchLambda;
03660     const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03661 
03662     if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
03663       FepNo = 0;
03664       exp_dE_ByRT = 0.0;
03665       net_dE = 0.0;
03666     }
03667     dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
03668          electEnergy - electEnergySlow - ljEnergy;
03669     BigReal RT = BOLTZMANN * simParams->alchTemp;
03670 
03671     if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
03672         (step / simParams->alchIDWSFreq) % 2 == 1 )) {
03673       // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
03674       FepNo++;
03675       exp_dE_ByRT += exp(-dE/RT);
03676       net_dE += dE;
03677     }
03678 
03679     if (simParams->alchOutFreq) {
03680       if (stepInRun == 0) {
03681         if (!fepFile.is_open()) {
03682           NAMD_backup_file(simParams->alchOutFile);
03683           fepFile.open(simParams->alchOutFile);
03684           iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
03685           if(alchEnsembleAvg){
03686             fepSum = 0.0;
03687             fepFile << "#            STEP                 Elec                            "
03688                     << "vdW                    dE           dE_avg         Temp             dG\n"
03689                     << "#                           l             l+dl      "
03690                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03691           }
03692           else{
03693             fepFile << "#            STEP                 Elec                            "
03694                     << "vdW                    dE         Temp\n"
03695                     << "#                           l             l+dl      "
03696                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03697           }
03698         }
03699         if(!step){
03700           fepFile << "#NEW FEP WINDOW: "
03701                   << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
03702                   << simParams->alchLambda2;
03703           if ( simParams->alchLambdaIDWS >= 0. ) {
03704             fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
03705           }
03706           fepFile << std::endl;
03707         }
03708       }
03709       if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03710         fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03711                 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03712                 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03713       }
03714       if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
03715         writeFepEnergyData(step, fepFile);
03716         fepFile.flush();
03717       }
03718       if (alchEnsembleAvg && (step == simParams->N)) {
03719         fepSum = fepSum + dG;
03720         fepFile << "#Free energy change for lambda window [ " << alchLambda
03721                 << " " << simParams->alchLambda2 << " ] is " << dG
03722                 << " ; net change until now is " << fepSum << std::endl;
03723         fepFile.flush();
03724       }
03725     }
03726   }
03727 }

void Controller::outputTiEnergy ( int  step  )  [protected]

Definition at line 3729 of file Controller.C.

References SimParameters::alchEquilSteps, SimParameters::alchLambda, SimParameters::alchLambdaFreq, SimParameters::alchOn, SimParameters::alchOutFile, SimParameters::alchOutFreq, SimParameters::alchTemp, SimParameters::alchThermIntOn, alchWork, bondedEnergy_ti_1, bondedEnergy_ti_2, electEnergy_ti_1, electEnergy_ti_2, electEnergyPME_ti_1, electEnergyPME_ti_2, electEnergySlow_ti_1, electEnergySlow_ti_2, endi(), SimParameters::firstTimestep, ofstream_namd::flush(), FORMAT(), SimParameters::getBondLambda(), SimParameters::getCurrentLambda2(), SimParameters::getElecLambda(), SimParameters::getLambdaDelta(), SimParameters::getVdwLambda(), iout, ofstream_namd::is_open(), ljEnergy_ti_1, ljEnergy_ti_2, NAMD_backup_file(), net_dEdl_bond_1, net_dEdl_bond_2, net_dEdl_elec_1, net_dEdl_elec_2, net_dEdl_lj_1, net_dEdl_lj_2, ofstream_namd::open(), recent_alchWork, recent_dEdl_bond_1, recent_dEdl_bond_2, recent_dEdl_elec_1, recent_dEdl_elec_2, recent_dEdl_lj_1, recent_dEdl_lj_2, recent_TiNo, simParams, tiFile, TiNo, and writeTiEnergyData().

Referenced by integrate().

03729                                         {
03730  if (simParams->alchOn && simParams->alchThermIntOn) {
03731   const int stepInRun = step - simParams->firstTimestep;
03732   const int alchEquilSteps = simParams->alchEquilSteps;
03733   const int alchLambdaFreq = simParams->alchLambdaFreq;
03734 
03735   if (stepInRun == 0 || stepInRun == alchEquilSteps) {
03736     TiNo = 0;
03737     net_dEdl_bond_1 = 0;
03738     net_dEdl_bond_2 = 0;
03739     net_dEdl_elec_1 = 0;
03740     net_dEdl_elec_2 = 0;
03741     net_dEdl_lj_1 = 0;
03742     net_dEdl_lj_2 = 0;
03743   }
03744   TiNo++;
03745   net_dEdl_bond_1 += bondedEnergy_ti_1;
03746   net_dEdl_bond_2 += bondedEnergy_ti_2;
03747   net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
03748                       + electEnergyPME_ti_1);
03749   net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
03750                       + electEnergyPME_ti_2);
03751   net_dEdl_lj_1 += ljEnergy_ti_1;
03752   net_dEdl_lj_2 += ljEnergy_ti_2;
03753 
03754   // Don't accumulate block averages (you'll get a divide by zero!) or even 
03755   // open the TI output if alchOutFreq is zero.
03756   if (simParams->alchOutFreq) {
03757     if (stepInRun == 0 || stepInRun == alchEquilSteps 
03758         || (! ((step - 1) % simParams->alchOutFreq))) {
03759       // output of instantaneous dU/dl now replaced with running average
03760       // over last alchOutFreq steps (except for step 0)
03761       recent_TiNo = 0;
03762       recent_dEdl_bond_1 = 0;
03763       recent_dEdl_bond_2 = 0;
03764       recent_dEdl_elec_1 = 0;
03765       recent_dEdl_elec_2 = 0;
03766       recent_dEdl_lj_1 = 0;
03767       recent_dEdl_lj_2 = 0;
03768       recent_alchWork = 0;
03769     }
03770     recent_TiNo++;
03771     recent_dEdl_bond_1 += bondedEnergy_ti_1;
03772     recent_dEdl_bond_2 += bondedEnergy_ti_2;
03773     recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 
03774                            + electEnergyPME_ti_1); 
03775     recent_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 
03776                            + electEnergyPME_ti_2); 
03777     recent_dEdl_lj_1 += ljEnergy_ti_1;
03778     recent_dEdl_lj_2 += ljEnergy_ti_2;
03779     recent_alchWork += alchWork;
03780 
03781     if (stepInRun == 0) {
03782       if (!tiFile.is_open()) {
03783         NAMD_backup_file(simParams->alchOutFile);
03784         tiFile.open(simParams->alchOutFile);
03785         /* BKR - This has been rather drastically updated to better match 
03786            stdout. This was necessary for several reasons:
03787            1) PME global scaling is obsolete (now removed)
03788            2) scaling of bonded terms was added
03789            3) alchemical work is now accumulated when switching is active
03790          */
03791         iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
03792         tiFile << "#TITITLE:    TS";
03793         tiFile << FORMAT("BOND1");
03794         tiFile << FORMAT("AVGBOND1");
03795         tiFile << FORMAT("ELECT1");
03796         tiFile << FORMAT("AVGELECT1");
03797         tiFile << "     ";
03798         tiFile << FORMAT("VDW1");
03799         tiFile << FORMAT("AVGVDW1");
03800         tiFile << FORMAT("BOND2");
03801         tiFile << FORMAT("AVGBOND2");
03802         tiFile << FORMAT("ELECT2");
03803         tiFile << "     ";
03804         tiFile << FORMAT("AVGELECT2");
03805         tiFile << FORMAT("VDW2");
03806         tiFile << FORMAT("AVGVDW2");
03807         if (alchLambdaFreq > 0) {
03808           tiFile << FORMAT("ALCHWORK");
03809           tiFile << FORMAT("CUMALCHWORK");
03810         }
03811         tiFile << std::endl;
03812       }
03813 
03814       if (alchLambdaFreq > 0) {
03815         tiFile << "#ALCHEMICAL SWITCHING ACTIVE " 
03816                << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
03817                << "\n#LAMBDA SCHEDULE: " 
03818                << "dL: " << simParams->getLambdaDelta() 
03819                << " Freq: " << alchLambdaFreq;
03820       }
03821       else {
03822         const BigReal alchLambda = simParams->alchLambda;    
03823         const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
03824         const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
03825         const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
03826         const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
03827         const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
03828         const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
03829         tiFile << "#NEW TI WINDOW: "
03830                << "LAMBDA " << alchLambda 
03831                << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
03832                << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
03833                << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
03834                << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
03835       }
03836       tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
03837              << std::endl;
03838     }
03839 
03840     if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03841       tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03842              << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03843              << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03844     }
03845     if ((step%simParams->alchOutFreq) == 0) {
03846       writeTiEnergyData(step, tiFile);
03847       tiFile.flush();
03848     }
03849   }
03850  }
03851 }

void Controller::printDynamicsEnergies ( int  step  )  [protected]

Definition at line 3003 of file Controller.C.

References compareChecksums(), Node::molecule, Node::Object(), printEnergies(), Node::simParameters, and state.

Referenced by integrate().

03003                                                {
03004 
03005     Node *node = Node::Object();
03006     Molecule *molecule = node->molecule;
03007     SimParameters *simParameters = node->simParameters;
03008     Lattice &lattice = state->lattice;
03009 
03010     compareChecksums(step);
03011 
03012     printEnergies(step,0);
03013 }

void Controller::printEnergies ( int  step,
int  minimize 
) [protected]

Definition at line 3015 of file Controller.C.

References Lattice::a(), Lattice::a_p(), SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchLambdaFreq, SimParameters::alchOn, SimParameters::alchThermIntOn, alchWork, avg_count, Lattice::b(), Lattice::b_p(), bondedEnergy_ti_1, bondedEnergy_ti_2, bondedEnergyDiff_f, Lattice::c(), Lattice::c_p(), CALLBACKDATA, CALLBACKLIST, computeAlchWork(), cumAlchWork, drudeBondTemp, drudeBondTempAvg, SimParameters::drudeOn, SimParameters::dt, IMDEnergies::Eangle, IMDEnergies::Ebond, IMDEnergies::Edihe, IMDEnergies::Eelec, IMDEnergies::Eimpr, electEnergy, electEnergy_f, electEnergy_ti_1, electEnergy_ti_2, electEnergyPME_ti_1, electEnergyPME_ti_2, electEnergySlow, electEnergySlow_f, electEnergySlow_ti_1, electEnergySlow_ti_2, endi(), IMDEnergies::Epot, ETITLE(), IMDEnergies::Etot, IMDEnergies::Evdw, FEPTITLE2(), fflush_count, SimParameters::firstLdbStep, SimParameters::firstTimestep, FORMAT(), IMDOutput::gather_energies(), GET_VECTOR, SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), PressureProfileReduction::getData(), Molecule::getEnergyTailCorr(), Node::getScript(), SimParameters::goForcesOn, SimParameters::goGroPair, goNativeEnergy, goNonnativeEnergy, goTotalEnergy, groGaussEnergy, groLJEnergy, groupPressure, groupPressure_avg, groupPressure_tavg, heat, iERROR(), iINFO(), Node::imd, SimParameters::IMDfreq, SimParameters::IMDon, iout, RequireReduction::item(), iWARN(), kineticEnergy, kineticEnergyCentered, kineticEnergyHalfstep, SimParameters::LJcorrection, ljEnergy, ljEnergy_f, ljEnergy_ti_1, ljEnergy_ti_2, marginViolations, memusage_MB(), SimParameters::mergeCrossterms, Node::molecule, multigatorCalcEnthalpy(), SimParameters::multigratorOn, SimParameters::N, NAMD_bug(), nbondFreq, Node::Object(), Lattice::origin(), SimParameters::outputEnergies, SimParameters::outputMomenta, SimParameters::outputPairlists, SimParameters::outputPressure, SimParameters::pairInteractionOn, pairlistWarnings, SimParameters::PMEOn, ppbonded, ppint, ppnonbonded, pressure, pressure_avg, pressure_tavg, PRESSUREFACTOR, pressureProfileAverage, pressureProfileCount, SimParameters::pressureProfileFreq, pressureProfileSlabs, printTiming(), SimParameters::qmForcesOn, reduction, REDUCTION_ANGLE_ENERGY, REDUCTION_BC_ENERGY, REDUCTION_BOND_ENERGY, REDUCTION_BONDED_ENERGY_F, REDUCTION_BONDED_ENERGY_TI_1, REDUCTION_BONDED_ENERGY_TI_2, REDUCTION_CROSSTERM_ENERGY, REDUCTION_DIHEDRAL_ENERGY, REDUCTION_ELECT_ENERGY, REDUCTION_ELECT_ENERGY_F, REDUCTION_ELECT_ENERGY_PME_TI_1, REDUCTION_ELECT_ENERGY_PME_TI_2, REDUCTION_ELECT_ENERGY_SLOW, REDUCTION_ELECT_ENERGY_SLOW_F, REDUCTION_ELECT_ENERGY_SLOW_TI_1, REDUCTION_ELECT_ENERGY_SLOW_TI_2, REDUCTION_ELECT_ENERGY_TI_1, REDUCTION_ELECT_ENERGY_TI_2, REDUCTION_GO_NATIVE_ENERGY, REDUCTION_GO_NONNATIVE_ENERGY, REDUCTION_GRO_GAUSS_ENERGY, REDUCTION_GRO_LJ_ENERGY, REDUCTION_IMPROPER_ENERGY, REDUCTION_LJ_ENERGY, REDUCTION_LJ_ENERGY_F, REDUCTION_LJ_ENERGY_TI_1, REDUCTION_LJ_ENERGY_TI_2, REDUCTION_MISC_ENERGY, Node::simParameters, simParams, slowFreq, ControllerState::smooth2_avg, state, stepInFullRun, SimParameters::stochRescaleHeat, SimParameters::stochRescaleOn, IMDEnergies::T, tavg_count, temp_avg, temperature, TITITLE(), totalEnergy, totalEnergy0, IMDEnergies::tstep, values, Lattice::volume(), Vector::x, XXXBIGREAL, Vector::y, and Vector::z.

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

03016 {
03017     Node *node = Node::Object();
03018     Molecule *molecule = node->molecule;
03019     SimParameters *simParameters = node->simParameters;
03020     Lattice &lattice = state->lattice;
03021 
03022     // Drude model ANISO energy is added into BOND energy
03023     // and THOLE energy is added into ELECT energy
03024 
03025     BigReal bondEnergy;
03026     BigReal angleEnergy;
03027     BigReal dihedralEnergy;
03028     BigReal improperEnergy;
03029     BigReal crosstermEnergy;
03030     BigReal boundaryEnergy;
03031     BigReal miscEnergy;
03032     BigReal potentialEnergy;
03033     BigReal flatEnergy;
03034     BigReal smoothEnergy;
03035     BigReal work;
03036 
03037     Vector momentum;
03038     Vector angularMomentum;
03039     BigReal volume = lattice.volume();
03040 
03041     bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
03042     angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
03043     dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
03044     improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
03045     crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
03046     boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
03047     miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
03048 
03049     if ( minimize || ! ( step % nbondFreq ) )
03050     {
03051       electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
03052       ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
03053 
03054       // JLai
03055       groLJEnergy = reduction->item(REDUCTION_GRO_LJ_ENERGY);
03056       groGaussEnergy = reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
03057 
03058       goNativeEnergy = reduction->item(REDUCTION_GO_NATIVE_ENERGY);
03059       goNonnativeEnergy = reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
03060       goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
03061 
03062 //fepb
03063       bondedEnergyDiff_f = reduction->item(REDUCTION_BONDED_ENERGY_F);
03064       electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
03065       ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
03066 
03067       bondedEnergy_ti_1 = reduction->item(REDUCTION_BONDED_ENERGY_TI_1);
03068       electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
03069       ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
03070       bondedEnergy_ti_2 = reduction->item(REDUCTION_BONDED_ENERGY_TI_2);
03071       electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
03072       ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
03073 //fepe
03074     }
03075 
03076     if ( minimize || ! ( step % slowFreq ) )
03077     {
03078       electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
03079 //fepb
03080       electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
03081 
03082       electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
03083       electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
03084       electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
03085       electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
03086 //fepe
03087     }
03088 
03089     if (simParameters->LJcorrection && volume) {
03090 #ifdef MEM_OPT_VERSION
03091       NAMD_bug("LJcorrection not supported in memory optimized build.");
03092 #else
03093       // Apply tail correction to energy.
03094       BigReal alchLambda = simParameters->getCurrentLambda(step);
03095       ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
03096 
03097       if (simParameters->alchOn) {
03098         if (simParameters->alchFepOn) {
03099           BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
03100           ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
03101         } else if (simParameters->alchThermIntOn) {
03102           ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
03103           ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
03104         }
03105       }
03106 #endif
03107     }
03108 
03109 //fepb BKR - Compute alchemical work if using dynamic lambda.  This is here so
03110 //           that the cumulative work can be given during a callback.
03111     if (simParameters->alchLambdaFreq > 0) {
03112       if (step <= 
03113           simParameters->firstTimestep + simParameters->alchEquilSteps) {
03114         cumAlchWork = 0.0;
03115       }
03116       alchWork = computeAlchWork(step);
03117       cumAlchWork += alchWork;
03118     }
03119 //fepe
03120 
03121     momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
03122     momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
03123     momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
03124     angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
03125     angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
03126     angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
03127 
03128     // Ported by JLai
03129     potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
03130   + improperEnergy + electEnergy + electEnergySlow + ljEnergy
03131   + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy 
03132         + groLJEnergy + groGaussEnergy);
03133     // End of port
03134     totalEnergy = potentialEnergy + kineticEnergy;
03135     flatEnergy = (totalEnergy
03136         + (1.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
03137     if ( !(step%slowFreq) ) {
03138       // only adjust based on most accurate energies
03139       BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
03140       if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
03141       if ( step != simParams->firstTimestep ) {
03142         smooth2_avg *= 0.9375;
03143         smooth2_avg += 0.0625 * s;
03144       }
03145     }
03146     smoothEnergy = (flatEnergy + smooth2_avg
03147         - (4.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
03148 
03149     // Reset values for accumulated heat and work.
03150     if (step <= simParameters->firstTimestep &&
03151         (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
03152       heat = 0.0;
03153       totalEnergy0 = totalEnergy;
03154     }
03155     if ( simParameters->outputMomenta && ! minimize &&
03156          ! ( step % simParameters->outputMomenta ) )
03157     {
03158       iout << "MOMENTUM: " << step 
03159            << " P: " << momentum
03160            << " L: " << angularMomentum
03161            << "\n" << endi;
03162     }
03163 
03164     if ( simParameters->outputPressure ) {
03165       pressure_tavg += pressure;
03166       groupPressure_tavg += groupPressure;
03167       tavg_count += 1;
03168       if ( minimize || ! ( step % simParameters->outputPressure ) ) {
03169         iout << "PRESSURE: " << step << " "
03170            << PRESSUREFACTOR * pressure << "\n"
03171            << "GPRESSURE: " << step << " "
03172            << PRESSUREFACTOR * groupPressure << "\n";
03173         if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
03174            << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
03175            << "GPRESSAVG: " << step << " "
03176            << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
03177         iout << endi;
03178         pressure_tavg = 0;
03179         groupPressure_tavg = 0;
03180         tavg_count = 0;
03181       }
03182     }
03183 
03184     // pressure profile reductions
03185     if (pressureProfileSlabs) {
03186       const int freq = simParams->pressureProfileFreq;
03187       const int arraysize = 3*pressureProfileSlabs;
03188       
03189       BigReal *total = new BigReal[arraysize];
03190       memset(total, 0, arraysize*sizeof(BigReal));
03191       const int first = simParams->firstTimestep;
03192 
03193       if (ppbonded)    ppbonded->getData(first, step, lattice, total);
03194       if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
03195       if (ppint)       ppint->getData(first, step, lattice, total);
03196       for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
03197       pressureProfileCount++;
03198 
03199       if (!(step % freq)) {
03200         // convert NAMD internal virial to pressure in units of bar 
03201         BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
03202    
03203         iout << "PRESSUREPROFILE: " << step << " ";
03204         if (step == first) {
03205           // output pressure profile for this step
03206           for (int i=0; i<arraysize; i++) {
03207             iout << total[i] * scalefac << " ";
03208           }
03209         } else {
03210           // output pressure profile averaged over the last count steps.
03211           scalefac /= pressureProfileCount;
03212           for (int i=0; i<arraysize; i++) 
03213             iout << pressureProfileAverage[i]*scalefac << " ";
03214         }
03215         iout << "\n" << endi; 
03216        
03217         // Clear the average for the next block
03218         memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
03219         pressureProfileCount = 0;
03220       }
03221       delete [] total;
03222     }
03223   
03224    if ( step != simParams->firstTimestep || stepInFullRun == 0 ) {  // skip repeated first step
03225     if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
03226      int benchPhase = stepInFullRun / simParams->firstLdbStep;
03227      if ( benchPhase > 0 && benchPhase < 10 ) {
03228 #ifdef NAMD_CUDA
03229       if ( simParams->outputEnergies < 60 ) {
03230         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
03231       }
03232 #endif
03233       iout << iINFO;
03234       if ( benchPhase < 4 ) iout << "Initial time: ";
03235       else iout << "Benchmark time: ";
03236       iout << CkNumPes() << " CPUs ";
03237       {
03238         BigReal wallPerStep =
03239     (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
03240   BigReal ns = simParams->dt / 1000000.0;
03241   BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
03242   BigReal daysPerNano = wallPerStep * days / ns;
03243   iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
03244         iout << memusage_MB() << " MB memory\n" << endi;
03245       }
03246      }
03247      startBenchTime = CmiWallTimer();
03248     }
03249     ++stepInFullRun;
03250    }
03251 
03252     printTiming(step);
03253 
03254     Vector pairVDWForce, pairElectForce;
03255     if ( simParameters->pairInteractionOn ) {
03256       GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
03257       GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
03258     }
03259 
03260     // Compute cumulative nonequilibrium work (including shadow work).
03261     if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03262       work = totalEnergy - totalEnergy0 - heat;
03263     }
03264 
03265     // callback to Tcl with whatever we can
03266 #ifdef NAMD_TCL
03267 #define CALLBACKDATA(LABEL,VALUE) \
03268                 labels << (LABEL) << " "; values << (VALUE) << " ";
03269 #define CALLBACKLIST(LABEL,VALUE) \
03270                 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
03271     if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
03272       std::ostringstream labels, values;
03273 #if CMK_BLUEGENEL
03274       // the normal version below gives a compiler error
03275       values.precision(16);
03276 #else
03277       values << std::setprecision(16);
03278 #endif
03279       CALLBACKDATA("TS",step);
03280       CALLBACKDATA("BOND",bondEnergy);
03281       CALLBACKDATA("ANGLE",angleEnergy);
03282       CALLBACKDATA("DIHED",dihedralEnergy);
03283       CALLBACKDATA("CROSS",crosstermEnergy);
03284       CALLBACKDATA("IMPRP",improperEnergy);
03285       CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
03286       CALLBACKDATA("VDW",ljEnergy);
03287       CALLBACKDATA("BOUNDARY",boundaryEnergy);
03288       CALLBACKDATA("MISC",miscEnergy);
03289       CALLBACKDATA("POTENTIAL",potentialEnergy);
03290       CALLBACKDATA("KINETIC",kineticEnergy);
03291       CALLBACKDATA("TOTAL",totalEnergy);
03292       CALLBACKDATA("TEMP",temperature);
03293       CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
03294       CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
03295       CALLBACKDATA("VOLUME",lattice.volume());
03296       CALLBACKLIST("CELL_A",lattice.a());
03297       CALLBACKLIST("CELL_B",lattice.b());
03298       CALLBACKLIST("CELL_C",lattice.c());
03299       CALLBACKLIST("CELL_O",lattice.origin());
03300       labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
03301                 << lattice.b_p() << " " << lattice.c_p() << "} ";
03302       if (simParameters->drudeOn) {
03303         CALLBACKDATA("DRUDEBOND",drudeBondTemp);
03304       }
03305       if ( simParameters->pairInteractionOn ) {
03306         CALLBACKLIST("VDW_FORCE",pairVDWForce);
03307         CALLBACKLIST("ELECT_FORCE",pairElectForce);
03308       }
03309       if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03310         CALLBACKLIST("HEAT",heat);
03311         CALLBACKLIST("WORK",work);
03312       }
03313       if (simParameters->alchOn) {
03314         if (simParameters->alchThermIntOn) {
03315           CALLBACKLIST("BOND1", bondedEnergy_ti_1);
03316           CALLBACKLIST("ELEC1", electEnergy_ti_1 + electEnergySlow_ti_1 +
03317                                 electEnergyPME_ti_1);
03318           CALLBACKLIST("VDW1", ljEnergy_ti_1);
03319           CALLBACKLIST("BOND2", bondedEnergy_ti_2);
03320           CALLBACKLIST("ELEC2", electEnergy_ti_2 + electEnergySlow_ti_2 +
03321                                 electEnergyPME_ti_2);
03322           CALLBACKLIST("VDW2", ljEnergy_ti_2);
03323           if (simParameters->alchLambdaFreq > 0) {
03324             CALLBACKLIST("CUMALCHWORK", cumAlchWork);
03325           }
03326         } else if (simParameters->alchFepOn) {
03327           CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
03328                                 improperEnergy + bondedEnergyDiff_f);
03329           CALLBACKLIST("ELEC2", electEnergy_f + electEnergySlow_f);
03330           CALLBACKLIST("VDW2", ljEnergy_f);
03331         } 
03332       }
03333 
03334       labels << '\0';  values << '\0';  // insane but makes Linux work
03335       state->callback_labelstring = labels.str();
03336       state->callback_valuestring = values.str();
03337       // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
03338     }
03339 #undef CALLBACKDATA
03340 #endif
03341 
03342     drudeBondTempAvg += drudeBondTemp;
03343 
03344     temp_avg += temperature;
03345     pressure_avg += trace(pressure)/3.;
03346     groupPressure_avg += trace(groupPressure)/3.;
03347     avg_count += 1;
03348 
03349     if ( simParams->outputPairlists && pairlistWarnings &&
03350                                 ! (step % simParams->outputPairlists) ) {
03351       iout << iINFO << pairlistWarnings <<
03352         " pairlist warnings in past " << simParams->outputPairlists <<
03353         " steps.\n" << endi;
03354       pairlistWarnings = 0;
03355     }
03356 
03357     BigReal enthalpy;
03358     if (simParameters->multigratorOn && ((step % simParameters->outputEnergies) == 0)) {
03359       enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
03360     }
03361 
03362     // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
03363     if ( ! minimize &&  step % simParameters->outputEnergies ) return;
03364     // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
03365 
03366     if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
03367       IMDEnergies energies;
03368       energies.tstep = step;
03369       energies.T = temp_avg/avg_count;
03370       energies.Etot = totalEnergy;
03371       energies.Epot = potentialEnergy;
03372       energies.Evdw = ljEnergy;
03373       energies.Eelec = electEnergy + electEnergySlow;
03374       energies.Ebond = bondEnergy;
03375       energies.Eangle = angleEnergy;
03376       energies.Edihe = dihedralEnergy + crosstermEnergy;
03377       energies.Eimpr = improperEnergy;
03378       Node::Object()->imd->gather_energies(&energies);
03379     }
03380 
03381     if ( marginViolations ) {
03382       iout << iERROR << marginViolations <<
03383         " margin violations detected since previous energy output.\n" << endi;
03384     }
03385     marginViolations = 0;
03386 
03387     if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
03388     {
03389         iout << "ETITLE:      TS";
03390         iout << FORMAT("BOND");
03391         iout << FORMAT("ANGLE");
03392         iout << FORMAT("DIHED");
03393         if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
03394         iout << FORMAT("IMPRP");
03395         iout << "     ";
03396         iout << FORMAT("ELECT");
03397         iout << FORMAT("VDW");
03398         iout << FORMAT("BOUNDARY");
03399         iout << FORMAT("MISC");
03400         iout << FORMAT("KINETIC");
03401         iout << "     ";
03402         iout << FORMAT("TOTAL");
03403         iout << FORMAT("TEMP");
03404         iout << FORMAT("POTENTIAL");
03405         // iout << FORMAT("TOTAL2");
03406         iout << FORMAT("TOTAL3");
03407         iout << FORMAT("TEMPAVG");
03408         if ( volume != 0. ) {
03409           iout << "     ";
03410           iout << FORMAT("PRESSURE");
03411           iout << FORMAT("GPRESSURE");
03412           iout << FORMAT("VOLUME");
03413           iout << FORMAT("PRESSAVG");
03414           iout << FORMAT("GPRESSAVG");
03415         }
03416         if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03417           iout << "     ";
03418           iout << FORMAT("HEAT");
03419           iout << FORMAT("WORK");
03420         }
03421         if (simParameters->drudeOn) {
03422           iout << "     ";
03423           iout << FORMAT("DRUDEBOND");
03424           iout << FORMAT("DRBONDAVG");
03425         }
03426         // Ported by JLai
03427         if (simParameters->goGroPair) {
03428           iout << "     ";
03429           iout << FORMAT("GRO_PAIR_LJ");
03430           iout << FORMAT("GRO_PAIR_GAUSS");
03431         }
03432 
03433         if (simParameters->goForcesOn) {
03434           iout << "     ";
03435           iout << FORMAT("NATIVE");
03436           iout << FORMAT("NONNATIVE");
03437           //iout << FORMAT("REL_NATIVE");
03438           //iout << FORMAT("REL_NONNATIVE");
03439           iout << FORMAT("GOTOTAL");
03440           //iout << FORMAT("GOAVG");
03441         }
03442         // End of port -- JLai
03443 
03444         if (simParameters->alchOn) {
03445           if (simParameters->alchThermIntOn) {
03446             iout << "\nTITITLE:     TS";
03447             iout << FORMAT("BOND1");
03448             iout << FORMAT("ELECT1");
03449             iout << FORMAT("VDW1");
03450             iout << FORMAT("BOND2");
03451             iout << "     ";
03452             iout << FORMAT("ELECT2");
03453             iout << FORMAT("VDW2");
03454             if (simParameters->alchLambdaFreq > 0) {
03455               iout << FORMAT("LAMBDA");
03456               iout << FORMAT("ALCHWORK");
03457               iout << FORMAT("CUMALCHWORK");
03458             }
03459           } else if (simParameters->alchFepOn) {
03460             iout << "\nFEPTITLE:    TS";
03461             iout << FORMAT("BOND2");
03462             iout << FORMAT("ELECT2");
03463             iout << FORMAT("VDW2");
03464             if (simParameters->alchLambdaFreq > 0) {
03465               iout << FORMAT("LAMBDA");
03466             }
03467           }
03468         }
03469 
03470         iout << "\n\n" << endi;
03471         
03472         if (simParameters->qmForcesOn) {
03473             iout << "QMETITLE:      TS";
03474             iout << FORMAT("QMID");
03475             iout << FORMAT("ENERGY");
03476             if (simParameters->PMEOn) iout << FORMAT("PMECORRENERGY");
03477             iout << "\n\n" << endi;
03478         }
03479         
03480     }
03481 
03482     // N.B.  HP's aCC compiler merges FORMAT calls in the same expression.
03483     //       Need separate statements because data returned in static array.
03484     iout << ETITLE(step);
03485     iout << FORMAT(bondEnergy);
03486     iout << FORMAT(angleEnergy);
03487     if ( simParameters->mergeCrossterms ) {
03488       iout << FORMAT(dihedralEnergy+crosstermEnergy);
03489     } else {
03490       iout << FORMAT(dihedralEnergy);
03491       iout << FORMAT(crosstermEnergy);
03492     }
03493     iout << FORMAT(improperEnergy);
03494     iout << "     ";
03495     iout << FORMAT(electEnergy+electEnergySlow);
03496     iout << FORMAT(ljEnergy);
03497     iout << FORMAT(boundaryEnergy);
03498     iout << FORMAT(miscEnergy);
03499     iout << FORMAT(kineticEnergy);
03500     iout << "     ";
03501     iout << FORMAT(totalEnergy);
03502     iout << FORMAT(temperature);
03503     iout << FORMAT(potentialEnergy);
03504     // iout << FORMAT(flatEnergy);
03505     iout << FORMAT(smoothEnergy);
03506     iout << FORMAT(temp_avg/avg_count);
03507     if ( volume != 0. )
03508     {
03509         iout << "     ";
03510         iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
03511         iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
03512         iout << FORMAT(volume);
03513         iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
03514         iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
03515     }
03516     if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
03517           iout << "     ";
03518           iout << FORMAT(heat);
03519           iout << FORMAT(work);
03520     }
03521     if (simParameters->drudeOn) {
03522         iout << "     ";
03523         iout << FORMAT(drudeBondTemp);
03524         iout << FORMAT(drudeBondTempAvg/avg_count);
03525     }
03526     // Ported by JLai
03527     if (simParameters->goGroPair) {
03528       iout << "     ";
03529       iout << FORMAT(groLJEnergy);
03530       iout << FORMAT(groGaussEnergy);
03531     }
03532 
03533     if (simParameters->goForcesOn) {
03534       iout << "     ";
03535       iout << FORMAT(goNativeEnergy);
03536       iout << FORMAT(goNonnativeEnergy);
03537       //iout << FORMAT(relgoNativeEnergy);
03538       //iout << FORMAT(relgoNonnativeEnergy);
03539       iout << FORMAT(goTotalEnergy);
03540       //iout << FORMAT("not implemented");
03541     } // End of port -- JLai
03542 
03543     if (simParameters->alchOn) { 
03544       if (simParameters->alchThermIntOn) {
03545         iout << "\n";
03546         iout << TITITLE(step);
03547         iout << FORMAT(bondedEnergy_ti_1);
03548         iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 + 
03549                        electEnergyPME_ti_1);
03550         iout << FORMAT(ljEnergy_ti_1);
03551         iout << FORMAT(bondedEnergy_ti_2);
03552         iout << "     ";
03553         iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 +
03554                        electEnergyPME_ti_2);
03555         iout << FORMAT(ljEnergy_ti_2);
03556         if (simParameters->alchLambdaFreq > 0) {
03557           iout << FORMAT(simParameters->getCurrentLambda(step));
03558           iout << FORMAT(alchWork);
03559           iout << FORMAT(cumAlchWork);
03560         }
03561       } else if (simParameters->alchFepOn) {
03562         iout << "\n";
03563         iout << FEPTITLE2(step);
03564         iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy 
03565                        + improperEnergy + bondedEnergyDiff_f);
03566         iout << FORMAT(electEnergy_f + electEnergySlow_f);
03567         iout << FORMAT(ljEnergy_f);
03568         if (simParameters->alchLambdaFreq > 0) {
03569           iout << FORMAT(simParameters->getCurrentLambda(step));
03570         }
03571       }
03572     }
03573 
03574     iout << "\n\n" << endi;
03575 
03576 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
03577      char webout[80];
03578      sprintf(webout,"%d %d %d %d",(int)totalEnergy,
03579              (int)(potentialEnergy),
03580              (int)kineticEnergy,(int)temperature);
03581      CApplicationDepositNode0Data(webout);
03582 #endif
03583 
03584     if (simParameters->pairInteractionOn) {
03585       iout << "PAIR INTERACTION:";
03586       iout << " STEP: " << step;
03587       iout << " VDW_FORCE: ";
03588       iout << FORMAT(pairVDWForce.x);
03589       iout << FORMAT(pairVDWForce.y);
03590       iout << FORMAT(pairVDWForce.z);
03591       iout << " ELECT_FORCE: ";
03592       iout << FORMAT(pairElectForce.x);
03593       iout << FORMAT(pairElectForce.y);
03594       iout << FORMAT(pairElectForce.z);
03595       iout << "\n" << endi;
03596     }
03597     drudeBondTempAvg = 0;
03598     temp_avg = 0;
03599     pressure_avg = 0;
03600     groupPressure_avg = 0;
03601     avg_count = 0;
03602 
03603     if ( fflush_count ) {
03604       --fflush_count;
03605       fflush(stdout);
03606     }
03607 }

void Controller::printFepMessage ( int  step  )  [protected]

Definition at line 1274 of file Controller.C.

References SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchLambda, SimParameters::alchLambda2, SimParameters::alchLambdaFreq, SimParameters::alchLambdaIDWS, SimParameters::alchOn, SimParameters::alchTemp, endi(), iout, and simParams.

Referenced by integrate().

01275 {
01276   if (simParams->alchOn && simParams->alchFepOn 
01277       && !simParams->alchLambdaFreq) {
01278     const BigReal alchLambda = simParams->alchLambda;
01279     const BigReal alchLambda2 = simParams->alchLambda2;
01280     const BigReal alchLambdaIDWS = simParams->alchLambdaIDWS;
01281     const BigReal alchTemp = simParams->alchTemp;
01282     const int alchEquilSteps = simParams->alchEquilSteps;
01283     iout << "FEP: RESETTING FOR NEW FEP WINDOW "
01284          << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2;
01285     if ( alchLambdaIDWS >= 0. ) {
01286       iout << " LAMBDA_IDWS " << alchLambdaIDWS;
01287     }
01288     iout << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
01289          << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
01290          << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp 
01291          << " K FOR FEP CALCULATION\n" << endi;
01292   }
01293 } 

void Controller::printMinimizeEnergies ( int  step  )  [protected]
void Controller::printTiMessage ( int  step  )  [protected]

Definition at line 1294 of file Controller.C.

References SimParameters::alchEquilSteps, SimParameters::alchLambda, SimParameters::alchLambdaFreq, SimParameters::alchOn, SimParameters::alchThermIntOn, endi(), iout, and simParams.

Referenced by integrate().

01295 {
01296   if (simParams->alchOn && simParams->alchThermIntOn 
01297       && !simParams->alchLambdaFreq) {
01298     const BigReal alchLambda = simParams->alchLambda;
01299     const int alchEquilSteps = simParams->alchEquilSteps;
01300     iout << "TI: RESETTING FOR NEW WINDOW "
01301          << "LAMBDA SET TO " << alchLambda 
01302          << "\nTI: WINDOW TO HAVE " << alchEquilSteps 
01303          << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
01304   }
01305 } 

void Controller::printTiming ( int  step  )  [protected]

Definition at line 2947 of file Controller.C.

References endi(), fflush_count, SimParameters::firstTimestep, iout, iWARN(), memusage_MB(), SimParameters::N, SimParameters::outputEnergies, SimParameters::outputTiming, and simParams.

Referenced by printEnergies().

02947                                      {
02948 
02949     if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
02950     {
02951       const double endWTime = CmiWallTimer() - firstWTime;
02952       const double endCTime = CmiTimer() - firstCTime;
02953 
02954       // fflush about once per minute
02955       if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
02956 
02957       const double elapsedW = 
02958         (endWTime - startWTime) / simParams->outputTiming;
02959       const double elapsedC = 
02960         (endCTime - startCTime) / simParams->outputTiming;
02961 
02962       const double remainingW = elapsedW * (simParams->N - step);
02963       const double remainingW_hours = remainingW / 3600;
02964 
02965       startWTime = endWTime;
02966       startCTime = endCTime;
02967 
02968 #ifdef NAMD_CUDA
02969       if ( simParams->outputEnergies < 60 &&
02970            step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
02971         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
02972       }
02973 #endif
02974       if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
02975         CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
02976                   ", %g hours remaining, %f MB of memory in use.\n",
02977                   step, endCTime, elapsedC, endWTime, elapsedW,
02978                   remainingW_hours, memusage_MB());
02979         if ( fflush_count ) { --fflush_count; fflush(stdout); }
02980       }
02981     }
02982 }

void Controller::reassignVelocities ( int  step  )  [protected]

Definition at line 1308 of file Controller.C.

References endi(), if(), iout, SimParameters::reassignFreq, SimParameters::reassignHold, SimParameters::reassignIncr, SimParameters::reassignTemp, and simParams.

Referenced by integrate().

01309 {
01310   const int reassignFreq = simParams->reassignFreq;
01311   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
01312     BigReal newTemp = simParams->reassignTemp;
01313     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
01314     if ( simParams->reassignIncr > 0.0 ) {
01315       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
01316         newTemp = simParams->reassignHold;
01317     } else {
01318       if ( newTemp < simParams->reassignHold )
01319         newTemp = simParams->reassignHold;
01320     }
01321     iout << "REASSIGNING VELOCITIES AT STEP " << step
01322          << " TO " << newTemp << " KELVIN.\n" << endi;
01323   }
01324 }

void Controller::rebalanceLoad ( int  step  )  [protected]

Definition at line 4105 of file Controller.C.

References fflush_count, LdbCoordinator::getNumStepsToRun(), ldbSteps, Node::Object(), LdbCoordinator::Object(), Node::outputPatchComputeMaps(), and LdbCoordinator::rebalance().

Referenced by integrate().

04106 {
04107   if ( ! ldbSteps ) { 
04108     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
04109   }
04110   if ( ! --ldbSteps ) {
04111     startBenchTime -= CmiWallTimer();
04112         Node::Object()->outputPatchComputeMaps("before_ldb", step);
04113     LdbCoordinator::Object()->rebalance(this);  
04114         startBenchTime += CmiWallTimer();
04115     fflush_count = 3;
04116   }
04117 }

void Controller::receivePressure ( int  step,
int  minimize = 0 
) [protected]

Definition at line 1404 of file Controller.C.

References SimParameters::accelMDDebugOn, SimParameters::accelMDOn, SimParameters::accelMDOutFreq, BOLTZMANN, calcPressure(), SimParameters::comMove, controlNumDegFreedom, controlPressure, controlPressure_nbond, controlPressure_normal, controlPressure_slow, drudeBondTemp, SimParameters::drudeOn, endi(), SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, SimParameters::fixedAtomsOn, GET_TENSOR, GET_VECTOR, groupPressure, groupPressure_nbond, groupPressure_normal, groupPressure_slow, iINFO(), iout, RequireReduction::item(), kineticEnergy, kineticEnergyCentered, kineticEnergyHalfstep, SimParameters::langevinOn, Node::molecule, momentumSqrSum, Molecule::num_deg_freedom(), Molecule::num_fixed_atoms(), Molecule::num_fixed_groups(), Molecule::num_group_deg_freedom(), Molecule::numAtoms, Molecule::numConstraints, numDegFreedom, Molecule::numDrudeAtoms, Molecule::numFepInitial, Molecule::numFixedAtoms, Molecule::numFixedGroups, Molecule::numFixedRigidBonds, Molecule::numHydrogenGroups, Molecule::numLonepairs, Molecule::numRigidBonds, Node::Object(), SimParameters::pairInteractionOn, pressure, pressure_amd, pressure_nbond, pressure_normal, pressure_slow, PRESSUREFACTOR, reduction, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, RequireReduction::require(), Node::simParameters, simParams, state, temperature, and SimParameters::useGroupPressure.

Referenced by integrate(), and printMinimizeEnergies().

01405 {
01406     Node *node = Node::Object();
01407     Molecule *molecule = node->molecule;
01408     SimParameters *simParameters = node->simParameters;
01409     Lattice &lattice = state->lattice;
01410 
01411     reduction->require();
01412 
01413     Tensor virial_normal;
01414     Tensor virial_nbond;
01415     Tensor virial_slow;
01416 #ifdef ALTVIRIAL
01417     Tensor altVirial_normal;
01418     Tensor altVirial_nbond;
01419     Tensor altVirial_slow;
01420 #endif
01421     Tensor intVirial_normal;
01422     Tensor intVirial_nbond;
01423     Tensor intVirial_slow;
01424     Vector extForce_normal;
01425     Vector extForce_nbond;
01426     Vector extForce_slow;
01427 
01428 #if 1
01429     numDegFreedom = molecule->num_deg_freedom();
01430     int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
01431     int numFixedGroups = molecule->num_fixed_groups();
01432     int numFixedAtoms = molecule->num_fixed_atoms();
01433 #endif
01434 #if 0
01435     int numAtoms = molecule->numAtoms;
01436     numDegFreedom = 3 * numAtoms;
01437     int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
01438     int numFixedAtoms =
01439         ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
01440     int numLonepairs = molecule->numLonepairs;
01441     int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
01442     if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
01443     if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01444     if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01445     if ( ! ( numFixedAtoms || molecule->numConstraints
01446         || simParameters->comMove || simParameters->langevinOn ) ) {
01447       numDegFreedom -= 3;
01448       numGroupDegFreedom -= 3;
01449     }
01450     if (simParameters->pairInteractionOn) {
01451       // this doesn't attempt to deal with fixed atoms or constraints
01452       numDegFreedom = 3 * molecule->numFepInitial;
01453     }
01454     int numRigidBonds = molecule->numRigidBonds;
01455     int numFixedRigidBonds =
01456         ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01457 
01458     // numLonepairs is subtracted here because all lonepairs have a rigid bond
01459     // to oxygen, but all of the LP degrees of freedom are dealt with above
01460     numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01461 #endif
01462 
01463     kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01464     kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01465 
01466     BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01467         reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01468     BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01469         reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01470 
01471     BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01472                                         / ( numDegFreedom * BOLTZMANN );
01473     BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01474                                         / ( numDegFreedom * BOLTZMANN );
01475     BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01476                                         / ( numGroupDegFreedom * BOLTZMANN );
01477     BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01478                                         / ( numGroupDegFreedom * BOLTZMANN );
01479 
01480     /*  test code for comparing different temperatures  
01481     iout << "TEMPTEST: " << step << " " << 
01482         atomTempHalfstep << " " <<
01483         atomTempCentered << " " <<
01484         groupTempHalfstep << " " <<
01485         groupTempCentered << "\n" << endi;
01486   iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
01487     numGroupDegFreedom << "\n" << endi;
01488      */
01489 
01490     GET_TENSOR(momentumSqrSum,reduction,REDUCTION_MOMENTUM_SQUARED);
01491 
01492     GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01493     GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01494     GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01495 
01496 #ifdef ALTVIRIAL
01497     GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01498     GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01499     GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01500 #endif
01501 
01502     GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01503     GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01504     GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01505 
01506     GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01507     GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01508     GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01509     // APH NOTE: These four lines are now done in calcPressure()
01510     // Vector extPosition = lattice.origin();
01511     // virial_normal -= outer(extForce_normal,extPosition);
01512     // virial_nbond -= outer(extForce_nbond,extPosition);
01513     // virial_slow -= outer(extForce_slow,extPosition);
01514 
01515     kineticEnergy = kineticEnergyCentered;
01516     temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );
01517 
01518     if (simParameters->drudeOn) {
01519       BigReal drudeComKE
01520         = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01521       BigReal drudeBondKE
01522         = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01523       int g_bond = 3 * molecule->numDrudeAtoms;
01524       int g_com = numDegFreedom - g_bond;
01525       temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
01526       drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
01527     }
01528 
01529     // Calculate number of degrees of freedom (controlNumDegFreedom)
01530     if ( simParameters->useGroupPressure )
01531     {
01532       controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01533       if ( ! ( numFixedAtoms || molecule->numConstraints
01534   || simParameters->comMove || simParameters->langevinOn ) ) {
01535         controlNumDegFreedom -= 1;
01536       }
01537     }
01538     else
01539     {
01540       controlNumDegFreedom = numDegFreedom / 3;
01541     }
01542     if (simParameters->fixCellDims) {
01543       if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01544       if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01545       if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01546     }
01547 
01548     // Calculate pressure tensors using virials
01549     calcPressure(step, minimize,
01550       virial_normal, virial_nbond, virial_slow,
01551       intVirial_normal, intVirial_nbond, intVirial_slow,
01552       extForce_normal, extForce_nbond, extForce_slow);
01553 
01554 #ifdef DEBUG_PRESSURE
01555     iout << iINFO << "Control pressure = " << controlPressure <<
01556       " = " << controlPressure_normal << " + " <<
01557       controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01558 #endif
01559     if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
01560         iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
01561              << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
01562              << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
01563              << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
01564              << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
01565              << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
01566              << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
01567              << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
01568              << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
01569              << endi;
01570    }
01571 }

void Controller::recvCheckpointAck ( checkpoint cp  )  [protected]

Definition at line 4096 of file Controller.C.

References checkpoint_task, Controller::checkpoint::lattice, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_SWAP, Controller::checkpoint::state, and state.

Referenced by Node::recvCheckpointAck().

04096                                                  {  // initiating replica
04097   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
04098     state->lattice = cp.lattice;
04099     *(ControllerState*)this = cp.state;
04100   }
04101   CkpvAccess(_qd)->process();
04102 }

void Controller::recvCheckpointReq ( const char *  key,
int  task,
checkpoint cp 
) [protected]

Definition at line 4066 of file Controller.C.

References checkpoints, NAMD_die(), SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_STORE, SCRIPT_CHECKPOINT_SWAP, and msm::swap().

Referenced by Node::recvCheckpointReq().

04066                                                                             {  // responding replica
04067   switch ( task ) {
04068     case SCRIPT_CHECKPOINT_STORE:
04069       if ( ! checkpoints.count(key) ) {
04070         checkpoints[key] = new checkpoint;
04071       }
04072       *checkpoints[key] = cp;
04073       break;
04074     case SCRIPT_CHECKPOINT_LOAD:
04075       if ( ! checkpoints.count(key) ) {
04076         NAMD_die("Unable to load checkpoint, requested key was never stored.");
04077       }
04078       cp = *checkpoints[key];
04079       break;
04080     case SCRIPT_CHECKPOINT_SWAP:
04081       if ( ! checkpoints.count(key) ) {
04082         NAMD_die("Unable to swap checkpoint, requested key was never stored.");
04083       }
04084       std::swap(cp,*checkpoints[key]);
04085       break;
04086     case SCRIPT_CHECKPOINT_FREE:
04087       if ( ! checkpoints.count(key) ) {
04088         NAMD_die("Unable to free checkpoint, requested key was never stored.");
04089       }
04090       delete checkpoints[key];
04091       checkpoints.erase(key);
04092       break;
04093   }
04094 }

void Controller::rescaleaccelMD ( int  step,
int  minimize = 0 
) [protected]

Definition at line 1823 of file Controller.C.

References SimParameters::accelMDalpha, SimParameters::accelMDDebugOn, SimParameters::accelMDdihe, SimParameters::accelMDdual, accelMDdVAverage, SimParameters::accelMDE, SimParameters::accelMDFirstStep, SimParameters::accelMDG, SimParameters::accelMDGcMDPrepSteps, SimParameters::accelMDGcMDSteps, SimParameters::accelMDGEquiPrepSteps, SimParameters::accelMDGEquiSteps, SimParameters::accelMDGiE, SimParameters::accelMDGresetVaftercmd, SimParameters::accelMDGRestart, SimParameters::accelMDGRestartFile, SimParameters::accelMDGSigma0D, SimParameters::accelMDGSigma0P, SimParameters::accelMDGStatWindow, SimParameters::accelMDLastStep, SimParameters::accelMDOn, SimParameters::accelMDOutFreq, ControllerBroadcasts::accelMDRescaleFactor, SimParameters::accelMDTalpha, SimParameters::accelMDTE, amd_reduction, broadcast, calc_accelMDG_E_k(), calc_accelMDG_force_factor(), calc_accelMDG_mean_std(), electEnergy, electEnergySlow, endi(), SimParameters::firstTimestep, GET_TENSOR, goNativeEnergy, goNonnativeEnergy, goTotalEnergy, groGaussEnergy, groLJEnergy, iout, RequireReduction::item(), SubmitReduction::item(), iWARN(), SimParameters::LJcorrection, ljEnergy, Node::molecule, SimParameters::N, NAMD_die(), nbondFreq, Node::Object(), PRESSUREFACTOR, SimpleBroadcastObject< T >::publish(), REDUCTION_ANGLE_ENERGY, REDUCTION_BC_ENERGY, REDUCTION_BOND_ENERGY, REDUCTION_CROSSTERM_ENERGY, REDUCTION_DIHEDRAL_ENERGY, REDUCTION_ELECT_ENERGY, REDUCTION_ELECT_ENERGY_SLOW, REDUCTION_GO_NATIVE_ENERGY, REDUCTION_GO_NONNATIVE_ENERGY, REDUCTION_GRO_GAUSS_ENERGY, REDUCTION_GRO_LJ_ENERGY, REDUCTION_IMPROPER_ENERGY, REDUCTION_LJ_ENERGY, REDUCTION_MAX_RESERVED, REDUCTION_MISC_ENERGY, RequireReduction::require(), SimParameters::restartFrequency, simParams, slowFreq, state, SubmitReduction::submit(), submit_reduction, Molecule::tail_corr_ener, virial_amd, Lattice::volume(), and write_accelMDG_rest_file().

Referenced by integrate(), and printMinimizeEnergies().

01824 {
01825     if ( !simParams->accelMDOn ) return;
01826 
01827     amd_reduction->require();
01828 
01829     // copy all to submit_reduction
01830     for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
01831       submit_reduction->item(i) += amd_reduction->item(i);
01832     }
01833     submit_reduction->submit();
01834 
01835     if (step == simParams->firstTimestep) accelMDdVAverage = 0;
01836 //    if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
01837     if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
01838 
01839     Node *node = Node::Object();
01840     Molecule *molecule = node->molecule;
01841     Lattice &lattice = state->lattice;
01842 
01843     const BigReal accelMDE = simParams->accelMDE;
01844     const BigReal accelMDalpha = simParams->accelMDalpha;
01845     const BigReal accelMDTE = simParams->accelMDTE;
01846     const BigReal accelMDTalpha = simParams->accelMDTalpha;
01847     const int accelMDOutFreq = simParams->accelMDOutFreq;
01848 
01849     //GaMD
01850     static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
01851     static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
01852     static BigReal k0P, kP, EP;
01853     static BigReal k0D, kD, ED;
01854     static int V_n = 1, iEusedD, iEusedP;
01855     static char warnD, warnP;
01856     static int restfreq;
01857 
01858     const int firststep = (simParams->accelMDFirstStep > simParams->firstTimestep) ? simParams->accelMDFirstStep : simParams->firstTimestep;
01859     const int ntcmdprep = (simParams->accelMDGcMDPrepSteps > 0) ? simParams->accelMDFirstStep + simParams->accelMDGcMDPrepSteps : 0;
01860     const int ntcmd     = simParams->accelMDFirstStep + simParams->accelMDGcMDSteps;
01861     const int ntebprep  = ntcmd + simParams->accelMDGEquiPrepSteps;
01862     const int nteb      = ntcmd + simParams->accelMDGEquiSteps;
01863     const int ntave     = simParams->accelMDGStatWindow;
01864 
01865     BigReal volume;
01866     BigReal bondEnergy;
01867     BigReal angleEnergy;
01868     BigReal dihedralEnergy;
01869     BigReal improperEnergy;
01870     BigReal crosstermEnergy;
01871     BigReal boundaryEnergy;
01872     BigReal miscEnergy;
01873     BigReal amd_electEnergy;
01874     BigReal amd_ljEnergy;
01875     BigReal amd_ljEnergy_Corr = 0.;
01876     BigReal amd_electEnergySlow;
01877     BigReal amd_groLJEnergy;
01878     BigReal amd_groGaussEnergy;
01879     BigReal amd_goTotalEnergy;
01880     BigReal potentialEnergy;
01881     BigReal factor_dihe = 1;
01882     BigReal factor_tot = 1;
01883     BigReal testV;
01884     BigReal dV = 0.;
01885     Vector  accelMDfactor;
01886     Tensor vir; //auto initialized to 0
01887     Tensor vir_dihe;
01888     Tensor vir_normal;
01889     Tensor vir_nbond;
01890     Tensor vir_slow;
01891 
01892     volume = lattice.volume();
01893 
01894     bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
01895     angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
01896     dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
01897     improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
01898     crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
01899     boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
01900     miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);
01901 
01902     GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
01903     GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
01904     GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
01905     GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
01906 
01907     if ( !( step % nbondFreq ) ) {
01908       amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
01909       amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
01910       amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
01911       amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
01912       BigReal goNativeEnergy = amd_reduction->item(REDUCTION_GO_NATIVE_ENERGY);
01913       BigReal goNonnativeEnergy = amd_reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
01914       amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
01915     } else {
01916       amd_electEnergy = electEnergy;
01917       amd_ljEnergy = ljEnergy;
01918       amd_groLJEnergy = groLJEnergy;
01919       amd_groGaussEnergy = groGaussEnergy;
01920       amd_goTotalEnergy = goTotalEnergy;
01921     }
01922 
01923     if( simParams->LJcorrection && volume ) {
01924         // Obtain tail correction (copied from printEnergies())
01925         // This value is only printed for sanity check
01926         // accelMD currently does not 'boost' tail correction
01927         amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
01928     }
01929 
01930     if ( !( step % slowFreq ) ) {
01931       amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
01932     } else {
01933       amd_electEnergySlow = electEnergySlow;
01934     }
01935 
01936     potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
01937         improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
01938         crosstermEnergy + boundaryEnergy + miscEnergy +
01939         amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
01940 
01941     //GaMD
01942     if(simParams->accelMDG){
01943        // if first time running accelMDG module
01944        if(step == firststep) {
01945           iEusedD = iEusedP = simParams->accelMDGiE;
01946           warnD   = warnP   = '\0';
01947 
01948           //restart file reading
01949           if(simParams->accelMDGRestart){
01950              FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
01951              char line[256];
01952              int dihe_n=0, tot_n=0;
01953              if(!rest){
01954                 sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
01955                 NAMD_die(line);
01956              }
01957 
01958              while(fgets(line, 256, rest)){
01959                 if(line[0] == 'D'){
01960                    dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
01961                                    &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
01962                 }
01963                 else if(line[0] == 'T'){
01964                    tot_n  = sscanf(line+1, " %d %la %la %la %la %la %la %la",
01965                                    &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
01966                 }
01967              }
01968 
01969              fclose(rest);
01970 
01971              V_n++;
01972 
01973              //check dihe settings
01974              if(simParams->accelMDdihe || simParams->accelMDdual){
01975                 if(dihe_n !=8)
01976                    NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
01977                 k0D = kD * (VmaxD - VminD);
01978                 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
01979                    << " Vmax " << VmaxD << " Vmin " << VminD 
01980                    << " Vavg " << VavgD << " sigmaV " << sigmaVD
01981                    << " E " << ED << " k " << kD
01982                    << "\n" << endi;
01983              }
01984 
01985              //check tot settings
01986              if(!simParams->accelMDdihe || simParams->accelMDdual){
01987                 if(tot_n !=8)
01988                    NAMD_die("Invalid total potential energy format in accelMDG restart file");
01989                 k0P = kP * (VmaxP - VminP);
01990                 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
01991                    << " Vmax " << VmaxP << " Vmin " << VminP 
01992                    << " Vavg " << VavgP << " sigmaV " << sigmaVP
01993                    << " E " << EP << " k " << kP
01994                    << "\n" << endi;
01995             }
01996 
01997             iEusedD = (ED == VmaxD) ? 1 : 2;
01998             iEusedP = (EP == VmaxP) ? 1 : 2;
01999           }
02000           //local restfreq follows NAMD restartfreq (default: 500)
02001           restfreq = simParams->restartFrequency ? simParams->restartFrequency : 500;
02002        }
02003 
02004        //for dihedral potential
02005        if(simParams->accelMDdihe || simParams->accelMDdual){
02006           testV = dihedralEnergy + crosstermEnergy;
02007 
02008           //write restart file every restartfreq steps
02009           if(step > firststep && step % restfreq == 0)
02010              write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
02011                    true, false);
02012           //write restart file at the end of the simulation
02013           if( simParams->accelMDLastStep > 0 ){
02014               if( step == simParams->accelMDLastStep )
02015                   write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
02016                           true, true);
02017           }
02018           else if(step == simParams->N)
02019               write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
02020                       true, true);
02021 
02022           //conventional MD
02023           if(step < ntcmd){
02024              //very first step
02025              if(step == firststep && !simParams->accelMDGRestart){
02026                 //initialize stat
02027                 VmaxD = VminD = VavgD = testV;
02028                 M2D = sigmaVD = 0.;
02029              }
02030              //first step after cmdprep
02031              else if(step == ntcmdprep && ntcmdprep != 0){
02032                 //reset stat
02033                 VmaxD = VminD = VavgD = testV;
02034                 M2D = sigmaVD = 0.;
02035                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
02036              }
02037              //every ntave steps
02038              else if(ntave > 0 && step % ntave == 0){
02039                 //update Vmax Vmin
02040                 if(testV > VmaxD) VmaxD = testV;
02041                 if(testV < VminD) VminD = testV;
02042                 //reset avg and std
02043                 VavgD = testV;
02044                 M2D = sigmaVD = 0.;
02045                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
02046              }
02047              //normal steps
02048              else
02049                 calc_accelMDG_mean_std(testV, V_n,
02050                       &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
02051 
02052              //last cmd step
02053              if(step == ntcmd - 1){
02054                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D, 
02055                       VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
02056              }
02057           }
02058           //equilibration
02059           else if(step < nteb){
02060              calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
02061                    &dV, &factor_dihe, &vir);
02062 
02063              //first step after cmd
02064              if(step == ntcmd && simParams->accelMDGresetVaftercmd){
02065                 //reset stat
02066                 VmaxD = VminD = VavgD = testV;
02067                 M2D = sigmaVD = 0.;
02068                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
02069              }
02070              //every ntave steps
02071              else if(ntave > 0 && step % ntave == 0){
02072                 //update Vmax Vmin
02073                 if(testV > VmaxD) VmaxD = testV;
02074                 if(testV < VminD) VminD = testV;
02075                 //reset avg and std
02076                 VavgD = testV;
02077                 M2D = sigmaVD = 0.;
02078                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
02079              }
02080              else
02081                 calc_accelMDG_mean_std(testV, V_n, 
02082                       &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
02083 
02084              //steps after ebprep
02085              if(step >= ntebprep)
02086                 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
02087                     calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D,
02088                           VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
02089           }
02090           //production
02091           else{
02092              calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
02093                    &dV, &factor_dihe, &vir);
02094           }
02095 
02096        }
02097        //for total potential
02098        if(!simParams->accelMDdihe || simParams->accelMDdual){
02099           Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
02100           testV = potentialEnergy;
02101           if(simParams->accelMDdual){
02102              testV -= dihedralEnergy + crosstermEnergy;
02103              vir_tot -= vir_dihe;
02104           }
02105 
02106           //write restart file every restartfreq steps
02107           if(step > firststep && step % restfreq == 0)
02108              write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
02109                    !simParams->accelMDdual, false);
02110           //write restart file at the end of simulation
02111           if( simParams->accelMDLastStep > 0 ){
02112               if( step == simParams->accelMDLastStep )
02113                   write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
02114                           !simParams->accelMDdual, true);
02115           }
02116           else if(step == simParams->N)
02117              write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
02118                    !simParams->accelMDdual, true);
02119 
02120           //conventional MD
02121           if(step < ntcmd){
02122              //very first step
02123              if(step == firststep && !simParams->accelMDGRestart){
02124                 //initialize stat
02125                 VmaxP = VminP = VavgP = testV;
02126                 M2P = sigmaVP = 0.;
02127              }
02128              //first step after cmdprep
02129              else if(step == ntcmdprep && ntcmdprep != 0){
02130                 //reset stat
02131                 VmaxP = VminP = VavgP = testV;
02132                 M2P = sigmaVP = 0.;
02133                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
02134              }
02135              //every ntave steps
02136              else if(ntave > 0 && step % ntave == 0){
02137                 //update Vmax Vmin
02138                 if(testV > VmaxP) VmaxP = testV;
02139                 if(testV < VminP) VminP = testV;
02140                 //reset avg and std
02141                 VavgP = testV;
02142                 M2P = sigmaVP = 0.;
02143                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
02144              }
02145              //normal steps
02146              else
02147                 calc_accelMDG_mean_std(testV, V_n,
02148                       &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
02149              //last cmd step
02150              if(step == ntcmd - 1){
02151                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P, 
02152                       VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
02153              }
02154           }
02155           //equilibration
02156           else if(step < nteb){
02157              calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
02158                    &dV, &factor_tot, &vir);
02159 
02160              //first step after cmd
02161              if(step == ntcmd && simParams->accelMDGresetVaftercmd){
02162                 //reset stat
02163                 VmaxP = VminP = VavgP = testV;
02164                 M2P = sigmaVP = 0.;
02165                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
02166              }
02167              //every ntave steps
02168              else if(ntave > 0 && step % ntave == 0){
02169                 //update Vmax Vmin
02170                 if(testV > VmaxP) VmaxP = testV;
02171                 if(testV < VminP) VminP = testV;
02172                 //reset avg and std
02173                 VavgP = testV;
02174                 M2P = sigmaVP = 0.;
02175                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
02176              }
02177              else
02178                 calc_accelMDG_mean_std(testV, V_n,
02179                       &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
02180 
02181              //steps after ebprep
02182              if(step >= ntebprep)
02183                 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
02184                     calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P,
02185                           VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
02186           }
02187           //production
02188           else{
02189              calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
02190                    &dV, &factor_tot, &vir);
02191           }
02192 
02193        }
02194        accelMDdVAverage += dV;
02195 
02196        //first step after ntcmdprep
02197        if((ntcmdprep > 0 && step == ntcmdprep) || 
02198           (step < nteb && ntave > 0 && step % ntave == 0) ||
02199           (simParams->accelMDGresetVaftercmd && step == ntcmd)){
02200           V_n = 1;
02201        }
02202 
02203        if(step < nteb)
02204            V_n++;
02205 
02206     }
02207     //aMD
02208     else{
02209         if (simParams->accelMDdihe) {
02210 
02211             testV = dihedralEnergy + crosstermEnergy;
02212             if ( testV < accelMDE ) {
02213                 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
02214                 factor_dihe *= factor_dihe;
02215                 vir = vir_dihe * (factor_dihe - 1.0);
02216                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02217                 accelMDdVAverage += dV;
02218             }  
02219 
02220         } else if (simParams->accelMDdual) {
02221 
02222             testV = dihedralEnergy + crosstermEnergy;
02223             if ( testV < accelMDE ) {
02224                 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
02225                 factor_dihe *= factor_dihe;
02226                 vir = vir_dihe * (factor_dihe - 1.0) ;
02227                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02228             }
02229 
02230             testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
02231             if ( testV < accelMDTE ) {
02232                 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
02233                 factor_tot *= factor_tot;
02234                 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
02235                 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
02236             }
02237             accelMDdVAverage += dV;
02238 
02239         } else {
02240 
02241             testV = potentialEnergy;
02242             if ( testV < accelMDE ) {
02243                 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
02244                 factor_tot *= factor_tot;
02245                 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
02246                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
02247                 accelMDdVAverage += dV;
02248             }
02249         } 
02250     }
02251 
02252     accelMDfactor[0]=factor_dihe;
02253     accelMDfactor[1]=factor_tot;
02254     accelMDfactor[2]=1;
02255     broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
02256     virial_amd = vir; 
02257 
02258     if ( factor_tot < 0.001 ) {
02259        iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
02260             << "\n" << endi;
02261     }
02262 
02263     if ( ! (step % accelMDOutFreq) ) {
02264         if ( !(step == simParams->firstTimestep) ) {
02265              accelMDdVAverage = accelMDdVAverage/accelMDOutFreq; 
02266         }
02267         iout << "ACCELERATED MD: STEP " << step
02268              << " dV "   << dV
02269              << " dVAVG " << accelMDdVAverage 
02270              << " BOND " << bondEnergy
02271              << " ANGLE " << angleEnergy
02272              << " DIHED " << dihedralEnergy+crosstermEnergy
02273              << " IMPRP " << improperEnergy
02274              << " ELECT " << amd_electEnergy+amd_electEnergySlow
02275              << " VDW "  << amd_ljEnergy
02276              << " POTENTIAL "  << potentialEnergy;
02277         if(amd_ljEnergy_Corr)
02278             iout << " LJcorr " << amd_ljEnergy_Corr;
02279         iout << "\n" << endi;
02280         if(simParams->accelMDG){
02281             if(simParams->accelMDdihe || simParams->accelMDdual)
02282                 iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD 
02283                     << " Vmax " << VmaxD << " Vmin " << VminD 
02284                     << " Vavg " << VavgD << " sigmaV " << sigmaVD
02285                     << " E " << ED << " k0 " << k0D << " k " << kD
02286                     << "\n" << endi;
02287             if(warnD & 1)
02288                 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
02289                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02290             if(warnD & 2)
02291                 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
02292                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02293             if(!simParams->accelMDdihe || simParams->accelMDdual)
02294                 iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP 
02295                     << " Vmax " << VmaxP << " Vmin " << VminP 
02296                     << " Vavg " << VavgP << " sigmaV " << sigmaVP
02297                     << " E " << EP << " k0 " << k0P << " k " << kP
02298                     << "\n" << endi;
02299             if(warnP & 1)
02300                 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
02301                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02302             if(warnP & 2)
02303                 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
02304                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
02305             warnD = warnP = '\0';
02306         }
02307 
02308         accelMDdVAverage = 0;
02309 
02310         if (simParams->accelMDDebugOn) {
02311            Tensor p_normal;
02312            Tensor p_nbond;
02313            Tensor p_slow;
02314            Tensor p;
02315            if ( volume != 0. )  {
02316                  p_normal = vir_normal/volume;
02317                  p_nbond  = vir_nbond/volume;
02318                  p_slow = vir_slow/volume;
02319                  p = vir/volume;
02320            }    
02321            iout << " accelMD Scaling Factor: " << accelMDfactor << "\n" 
02322                 << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
02323                 << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
02324                 << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
02325                 << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n" 
02326                 << endi;
02327         }
02328    }
02329 }

void Controller::rescaleVelocities ( int  step  )  [protected]

Definition at line 1233 of file Controller.C.

References SimParameters::adaptTempFirstStep, SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, SimParameters::adaptTempRescale, adaptTempT, broadcast, SimpleBroadcastObject< T >::publish(), SimParameters::rescaleFreq, SimParameters::rescaleTemp, rescaleVelocities_numTemps, rescaleVelocities_sumTemps, simParams, temperature, and ControllerBroadcasts::velocityRescaleFactor.

Referenced by integrate().

01234 {
01235   const int rescaleFreq = simParams->rescaleFreq;
01236   if ( rescaleFreq > 0 ) {
01237     rescaleVelocities_sumTemps += temperature;  ++rescaleVelocities_numTemps;
01238     if ( rescaleVelocities_numTemps == rescaleFreq ) {
01239       BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
01240       BigReal rescaleTemp = simParams->rescaleTemp;
01241       if ( simParams->adaptTempOn && simParams->adaptTempRescale && (step > simParams->adaptTempFirstStep ) && 
01242                 (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
01243         rescaleTemp = adaptTempT;
01244       }
01245       BigReal factor = sqrt(rescaleTemp/avgTemp);
01246       broadcast->velocityRescaleFactor.publish(step,factor);
01247       //iout << "RESCALING VELOCITIES AT STEP " << step
01248       //     << " FROM AVERAGE TEMPERATURE OF " << avgTemp
01249       //     << " TO " << rescaleTemp << " KELVIN.\n" << endi;
01250       rescaleVelocities_sumTemps = 0;  rescaleVelocities_numTemps = 0;
01251     }
01252   }
01253 }

void Controller::resumeAfterTraceBarrier ( int  step  ) 

Definition at line 4136 of file Controller.C.

References awaken(), broadcast, SimpleBroadcastObject< T >::publish(), and ControllerBroadcasts::traceBarrier.

Referenced by Node::resumeAfterTraceBarrier().

04136                                                 {
04137         broadcast->traceBarrier.publish(step,1);
04138         CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime); 
04139         awaken();
04140 }

void Controller::run ( void   ) 

Definition at line 268 of file Controller.C.

References awaken(), CTRL_STK_SZ, and DebugM.

Referenced by NamdState::runController().

00269 {
00270     // create a Thread and invoke it
00271     DebugM(4, "Starting thread in controller on this=" << this << "\n");
00272     thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
00273     CthSetStrategyDefault(thread);
00274 #if CMK_BLUEGENE_CHARM
00275     BgAttach(thread);
00276 #endif
00277     awaken();
00278 }

void Controller::stochRescaleVelocities ( int  step  )  [protected]

The Controller routine for stochastic velocity rescaling uses the most recent temperature reduction to calculate the velocity rescaling coefficient that is then broadcast to all patches.

Generate and broadcast the scale factor for stochastic velocity rescaling.

Stochastic velocity rescaling couples the system to a heat bath by globally scaling the velocites by a single factor. This factor is chosen based on the instantaneous and bath temperatures, a user-defined time scale, and a stochastic component linked to the number of degrees of freedom in the system. All of this information is combined here and sent to the Sequencer for the actual rescaling.

Parameters:
step the current timestep

Definition at line 1348 of file Controller.C.

References BOLTZMANN, broadcast, Random::gaussian(), heat, numDegFreedom, SimpleBroadcastObject< T >::publish(), random, simParams, stochRescale_count, ControllerBroadcasts::stochRescaleCoefficient, SimParameters::stochRescaleFreq, SimParameters::stochRescaleOn, SimParameters::stochRescaleTemp, stochRescaleTimefactor, Random::sum_of_squared_gaussians(), and temperature.

Referenced by integrate().

01349 {
01350   if ( simParams->stochRescaleOn )
01351   {
01352     ++stochRescale_count;
01353     if ( stochRescale_count == simParams->stochRescaleFreq )
01354     { 
01355       const BigReal stochRescaleTemp = simParams->stochRescaleTemp;
01356 
01357       BigReal coefficient = 1.;
01358       if ( temperature > 0. ) 
01359       {
01360         BigReal R1 = random->gaussian();
01361         // BigReal gammaShape = 0.5*(numDegFreedom - 1);
01362         // R2sum is the sum of (numDegFreedom - 1) squared normal variables, which is
01363         // chi-squared distributed. This is in turn a special case of the Gamma
01364         // distribution, which converges to a normal distribution in the limit of a
01365         // large shape parameter.
01366         // BigReal R2sum = 2*(gammaShape + sqrt(gammaShape)*random->gaussian()) + R1*R1;
01367         BigReal R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
01368         BigReal tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
01369 
01370         coefficient = sqrt(stochRescaleTimefactor + (1 - stochRescaleTimefactor)*tempfactor*R2sum
01371                   + 2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*stochRescaleTimefactor)); 
01372       }
01373       broadcast->stochRescaleCoefficient.publish(step,coefficient);
01374       heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient - 1.0);
01375       stochRescale_count = 0;
01376     }
01377   }
01378 }

void Controller::tcoupleVelocities ( int  step  )  [protected]

Definition at line 1326 of file Controller.C.

References broadcast, SimpleBroadcastObject< T >::publish(), simParams, ControllerBroadcasts::tcoupleCoefficient, SimParameters::tCoupleOn, SimParameters::tCoupleTemp, and temperature.

Referenced by integrate().

01327 {
01328   if ( simParams->tCoupleOn )
01329   {
01330     const BigReal tCoupleTemp = simParams->tCoupleTemp;
01331     BigReal coefficient = 1.;
01332     if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
01333     broadcast->tcoupleCoefficient.publish(step,coefficient);
01334   }
01335 }

void Controller::terminate ( void   )  [protected]

Definition at line 4157 of file Controller.C.

References awaken().

Referenced by algorithm().

04157                                {
04158   BackEnd::awaken();
04159   CthFree(thread);
04160   CthSuspend();
04161 }

void Controller::traceBarrier ( int  turnOnTrace,
int  step 
) [protected]

Definition at line 4129 of file Controller.C.

Referenced by integrate().

04129                                                        {
04130         CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);       
04131         CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04132         nd.traceBarrier(turnOnTrace, step);
04133         CthSuspend();
04134 }

void Controller::write_accelMDG_rest_file ( int  step_n,
char  type,
int  V_n,
BigReal  Vmax,
BigReal  Vmin,
BigReal  Vavg,
BigReal  sigmaV,
BigReal  M2,
BigReal  E,
BigReal  k,
bool  write_topic,
bool  lasttime 
) [protected]

Definition at line 1765 of file Controller.C.

References endi(), iout, iWARN(), NAMD_backup_file(), and simParams.

Referenced by rescaleaccelMD().

01765                                                                                                                                                              {
01766    FILE *rest;
01767    char timestepstr[20];
01768    char *filename;
01769    int baselen;
01770    char *restart_name;
01771    const char *bsuffix;
01772 
01773    if(lasttime || simParams->restartFrequency == 0){
01774            filename = simParams->outputFilename;
01775            bsuffix = ".BAK";
01776    }
01777    else{
01778            filename = simParams->restartFilename;
01779            bsuffix = ".old";
01780    }
01781 
01782    baselen = strlen(filename);
01783    restart_name = new char[baselen+26];
01784 
01785    strcpy(restart_name, filename);
01786    if ( simParams->restartSave && !lasttime) {
01787       sprintf(timestepstr,".%d",step_n);
01788       strcat(restart_name, timestepstr);
01789       bsuffix = ".BAK";
01790    }
01791    strcat(restart_name, ".gamd");
01792 
01793    if(write_topic){
01794       NAMD_backup_file(restart_name,bsuffix);
01795 
01796       rest = fopen(restart_name, "w");
01797       if(rest){
01798          fprintf(rest, "# NAMD accelMDG restart file\n"
01799                "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
01800                "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01801                type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01802          fclose(rest);
01803          iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
01804       }
01805       else
01806          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01807    }
01808    else{
01809       rest = fopen(restart_name, "a");
01810       if(rest){
01811          fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01812                           type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01813          fclose(rest);
01814       }
01815       else
01816          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01817    }
01818 
01819    delete[] restart_name;
01820 }

void Controller::writeExtendedSystemData ( int  step,
ofstream_namd file 
) [protected]

Definition at line 3622 of file Controller.C.

References Lattice::a(), Lattice::a_p(), Lattice::b(), Lattice::b_p(), Lattice::c(), Lattice::c_p(), ControllerState::langevinPiston_strainRate, SimParameters::langevinPistonOn, Lattice::origin(), simParams, state, Vector::x, Vector::y, and Vector::z.

Referenced by outputExtendedSystem().

03622                                                                       {
03623   Lattice &lattice = state->lattice;
03624   file.precision(12);
03625   file << step;
03626     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
03627     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
03628     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
03629     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
03630   if ( simParams->langevinPistonOn ) {
03631     Vector strainRate = diagonal(langevinPiston_strainRate);
03632     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
03633     file << " " << strainRate.x;
03634     file << " " << strainRate.y;
03635     file << " " << strainRate.z;
03636     file << " " << strainRate2.x;
03637     file << " " << strainRate2.y;
03638     file << " " << strainRate2.z;
03639   }
03640   file << std::endl;
03641 }

void Controller::writeExtendedSystemLabels ( ofstream_namd file  )  [protected]

Definition at line 3609 of file Controller.C.

References Lattice::a_p(), Lattice::b_p(), Lattice::c_p(), SimParameters::langevinPistonOn, simParams, and state.

Referenced by outputExtendedSystem().

03609                                                               {
03610   Lattice &lattice = state->lattice;
03611   file << "#$LABELS step";
03612   if ( lattice.a_p() ) file << " a_x a_y a_z";
03613   if ( lattice.b_p() ) file << " b_x b_y b_z";
03614   if ( lattice.c_p() ) file << " c_x c_y c_z";
03615   file << " o_x o_y o_z";
03616   if ( simParams->langevinPistonOn ) {
03617     file << " s_x s_y s_z s_u s_v s_w";
03618   }
03619   file << std::endl;
03620 }

void Controller::writeFepEnergyData ( int  step,
ofstream_namd file 
) [protected]

Definition at line 3887 of file Controller.C.

References SimParameters::alchEnsembleAvg, SimParameters::alchIDWSFreq, SimParameters::alchLambda, SimParameters::alchLambdaIDWS, SimParameters::alchTemp, BOLTZMANN, dE, dG, electEnergy, electEnergy_f, electEnergySlow, electEnergySlow_f, exp_dE_ByRT, fepFile, FepNo, FEPTITLE(), FEPTITLE_BACK(), SimParameters::firstTimestep, FORMAT(), ljEnergy, ljEnergy_f, net_dE, simParams, and temperature.

Referenced by outputFepEnergy().

03887                                                                  {
03888   BigReal eeng = electEnergy + electEnergySlow;
03889   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
03890   BigReal RT = BOLTZMANN * simParams->alchTemp;
03891 
03892   const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03893   const int stepInRun = step - simParams->firstTimestep;
03894   const BigReal alchLambda = simParams->alchLambda;
03895 
03896   if(stepInRun){
03897     if ( simParams->alchLambdaIDWS >= 0. &&
03898         (step / simParams->alchIDWSFreq) % 2 == 0 ) {
03899       // IDWS is active and we are on a "backward" timestep
03900       fepFile << FEPTITLE_BACK(step);
03901     } else {
03902       // "forward" timestep
03903       fepFile << FEPTITLE(step);
03904     }
03905     fepFile << FORMAT(eeng);
03906     fepFile << FORMAT(eeng_f);
03907     fepFile << FORMAT(ljEnergy);
03908     fepFile << FORMAT(ljEnergy_f);
03909     fepFile << FORMAT(dE);
03910     if(alchEnsembleAvg){
03911       BigReal dE_avg = net_dE / FepNo;
03912       fepFile << FORMAT(dE_avg);
03913     }
03914     fepFile << FORMAT(temperature);
03915     if(alchEnsembleAvg){
03916       dG = -(RT * log(exp_dE_ByRT / FepNo));
03917       fepFile << FORMAT(dG);
03918     }
03919     fepFile << std::endl;
03920   }
03921 }

void Controller::writeTiEnergyData ( int  step,
ofstream_namd file 
) [protected]

Friends And Related Function Documentation

friend class CheckpointMsg [friend]

Definition at line 54 of file Controller.h.

friend class Node [friend]

Definition at line 53 of file Controller.h.

friend class ScriptTcl [friend]

Definition at line 52 of file Controller.h.


Member Data Documentation

Definition at line 294 of file Controller.h.

Referenced by rescaleaccelMD().

Definition at line 317 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 311 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 310 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 306 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

int Controller::adaptTempBin [protected]

Definition at line 312 of file Controller.h.

Referenced by adaptTempUpdate().

int Controller::adaptTempBins [protected]

Definition at line 313 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 315 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 314 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 316 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 308 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 309 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 319 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 318 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 303 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 301 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 300 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 305 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 304 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 302 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart().

Definition at line 320 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempWriteRestart().

Definition at line 150 of file Controller.h.

Referenced by outputTiEnergy(), and printEnergies().

Definition at line 231 of file Controller.h.

Referenced by Controller(), rescaleaccelMD(), and ~Controller().

int Controller::avg_count [protected]

Definition at line 83 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 125 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 126 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 112 of file Controller.h.

Referenced by outputFepEnergy(), and printEnergies().

Definition at line 262 of file Controller.h.

Referenced by algorithm().

Definition at line 263 of file Controller.h.

Referenced by algorithm().

Definition at line 261 of file Controller.h.

Referenced by algorithm(), and Controller().

int Controller::checkpoint_task [protected]

Definition at line 270 of file Controller.h.

Referenced by recvCheckpointAck().

std::map<std::string,checkpoint*> Controller::checkpoints [protected]

Definition at line 269 of file Controller.h.

Referenced by algorithm(), and recvCheckpointReq().

Definition at line 242 of file Controller.h.

Referenced by enqueueCollections().

int Controller::computeChecksum [protected]

Definition at line 88 of file Controller.h.

Referenced by compareChecksums().

Definition at line 76 of file Controller.h.

Referenced by calcPressure(), langevinPiston1(), langevinPiston2(), and receivePressure().

Definition at line 75 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 77 of file Controller.h.

Referenced by calcPressure(), langevinPiston1(), langevinPiston2(), and receivePressure().

Definition at line 139 of file Controller.h.

Referenced by Controller(), printEnergies(), and writeTiEnergyData().

BigReal Controller::dE [protected]

Definition at line 118 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::dG [protected]

Definition at line 120 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 154 of file Controller.h.

Referenced by Controller(), printEnergies(), and receivePressure().

Definition at line 155 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 103 of file Controller.h.

Referenced by outputFepEnergy(), printEnergies(), rescaleaccelMD(), and writeFepEnergyData().

Definition at line 113 of file Controller.h.

Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData().

Definition at line 127 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 130 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 140 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 141 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 104 of file Controller.h.

Referenced by outputFepEnergy(), printEnergies(), rescaleaccelMD(), and writeFepEnergyData().

Definition at line 114 of file Controller.h.

Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData().

Definition at line 128 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 131 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 117 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 251 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

int Controller::FepNo [protected]

Definition at line 121 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 123 of file Controller.h.

Referenced by outputFepEnergy().

int Controller::fflush_count [protected]

Definition at line 215 of file Controller.h.

Referenced by printEnergies(), printTiming(), and rebalanceLoad().

Definition at line 108 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 109 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 110 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 107 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 106 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 167 of file Controller.h.

Referenced by calcPressure(), printEnergies(), and receivePressure().

Definition at line 82 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 73 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 72 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 74 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 85 of file Controller.h.

Referenced by Controller(), and printEnergies().

heat exchanged with the thermostat since firstTimestep

Definition at line 162 of file Controller.h.

Referenced by Controller(), printEnergies(), and stochRescaleVelocities().

Definition at line 159 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Definition at line 158 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Definition at line 197 of file Controller.h.

Referenced by Controller().

int Controller::ldbSteps [protected]

Definition at line 213 of file Controller.h.

Referenced by rebalanceLoad().

Definition at line 105 of file Controller.h.

Referenced by outputFepEnergy(), printEnergies(), rescaleaccelMD(), and writeFepEnergyData().

Definition at line 115 of file Controller.h.

Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData().

Definition at line 116 of file Controller.h.

Definition at line 129 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 132 of file Controller.h.

Referenced by computeAlchWork(), outputTiEnergy(), and printEnergies().

Definition at line 89 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Definition at line 93 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 94 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 95 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

int Controller::min_huge_count [protected]

Definition at line 97 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 59 of file Controller.h.

Referenced by Controller(), minimize(), and ~Controller().

Definition at line 96 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 204 of file Controller.h.

Referenced by multigratorPressure(), multigratorTemperature(), and receivePressure().

std::vector<BigReal> Controller::multigratorNu [protected]

Definition at line 206 of file Controller.h.

Referenced by Controller(), multigatorCalcEnthalpy(), and multigratorTemperature().

std::vector<BigReal> Controller::multigratorNuT [protected]

Definition at line 207 of file Controller.h.

Referenced by Controller(), and multigratorTemperature().

std::vector<BigReal> Controller::multigratorOmega [protected]

Definition at line 208 of file Controller.h.

Referenced by Controller(), multigatorCalcEnthalpy(), and multigratorTemperature().

Definition at line 210 of file Controller.h.

Referenced by Controller(), multigratorTemperature(), and ~Controller().

Definition at line 202 of file Controller.h.

Referenced by Controller(), multigatorCalcEnthalpy(), and multigratorPressure().

Definition at line 203 of file Controller.h.

Referenced by multigratorPressure().

std::vector<BigReal> Controller::multigratorZeta [protected]

Definition at line 209 of file Controller.h.

Referenced by Controller(), multigatorCalcEnthalpy(), and multigratorTemperature().

int Controller::nbondFreq [protected]

Definition at line 119 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 133 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 134 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 135 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 136 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 137 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 138 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int64_t Controller::numDegFreedom [protected]

Definition at line 274 of file Controller.h.

Referenced by Controller().

Definition at line 90 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Definition at line 199 of file Controller.h.

Referenced by langevinPiston1().

Definition at line 235 of file Controller.h.

Referenced by Controller(), printEnergies(), and ~Controller().

Definition at line 237 of file Controller.h.

Referenced by Controller(), printEnergies(), and ~Controller().

Definition at line 236 of file Controller.h.

Referenced by Controller(), printEnergies(), and ~Controller().

Definition at line 166 of file Controller.h.

Referenced by calcPressure(), printEnergies(), and receivePressure().

Definition at line 70 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 81 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 68 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 67 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 69 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 84 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 240 of file Controller.h.

Referenced by Controller(), printEnergies(), and ~Controller().

Definition at line 239 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 238 of file Controller.h.

Referenced by Controller(), and printEnergies().

Random* Controller::random [protected]

Definition at line 149 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 143 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 144 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 145 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 146 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 147 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 148 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::recent_TiNo [protected]

Definition at line 151 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 174 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

Definition at line 173 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

int Controller::slowFreq [protected]

Definition at line 165 of file Controller.h.

NamdState* const Controller::state [protected]
int Controller::stepInFullRun [protected]

Definition at line 101 of file Controller.h.

Referenced by printEnergies().

Count time steps until next stochastic velocity rescaling.

Definition at line 186 of file Controller.h.

Referenced by Controller(), and stochRescaleVelocities().

The timefactor for stochastic velocity rescaling depends on fixed configuration parameters, so can be precomputed.

Definition at line 189 of file Controller.h.

Referenced by Controller(), and stochRescaleVelocities().

Definition at line 198 of file Controller.h.

Referenced by langevinPiston1().

Definition at line 232 of file Controller.h.

Referenced by Controller(), rescaleaccelMD(), and ~Controller().

int Controller::tavg_count [protected]

Definition at line 86 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 80 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 255 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::TiNo [protected]

Definition at line 142 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 102 of file Controller.h.

Referenced by adaptTempUpdate(), printEnergies(), and printMinimizeEnergies().

totalEnergy at firstTimestep

Definition at line 163 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 71 of file Controller.h.

Referenced by calcPressure(), and rescaleaccelMD().

Definition at line 245 of file Controller.h.

Referenced by outputExtendedSystem().


The documentation for this class was generated from the following files:

Generated on 20 Oct 2019 for NAMD by  doxygen 1.6.1