Controller Class Reference

#include <Controller.h>

Inheritance diagram for Controller:

ControllerState List of all members.

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 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

Classes

struct  checkpoint

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

Controller::~Controller ( void   )  [virtual]

Definition at line 247 of file Controller.C.

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

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


Member Function Documentation

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

Definition at line 2284 of file Controller.C.

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

Referenced by adaptTempUpdate().

02284                                        {
02285     if (!simParams->adaptTempOn) return;
02286     iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
02287     adaptTempDtMin = 0;
02288     adaptTempDtMax = 0;
02289     adaptTempAutoDt = false;
02290     if (simParams->adaptTempBins == 0) {
02291       iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
02292       std::ifstream adaptTempRead(simParams->adaptTempInFile);
02293       if (adaptTempRead) {
02294         int readInt;
02295         BigReal readReal;
02296         bool readBool;
02297         // step
02298         adaptTempRead >> readInt;
02299         // Start with min and max temperatures
02300         adaptTempRead >> adaptTempT;     // KELVIN
02301         adaptTempRead >> adaptTempBetaMin;  // KELVIN
02302         adaptTempRead >> adaptTempBetaMax;  // KELVIN
02303         adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
02304         adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
02305         // In case file is manually edited
02306         if (adaptTempBetaMin > adaptTempBetaMax){
02307             readReal = adaptTempBetaMax;
02308             adaptTempBetaMax = adaptTempBetaMin;
02309             adaptTempBetaMin = adaptTempBetaMax;
02310         }
02311         adaptTempRead >> adaptTempBins;     
02312         adaptTempRead >> adaptTempCg;
02313         adaptTempRead >> adaptTempDt;
02314         adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02315         adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02316         adaptTempPotEnergySamples = new int[adaptTempBins];
02317         adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02318         adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02319         adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02320         adaptTempBetaN           = new BigReal[adaptTempBins];
02321         adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02322         for(int j = 0; j < adaptTempBins; ++j) {
02323           adaptTempRead >> adaptTempPotEnergyAve[j];
02324           adaptTempRead >> adaptTempPotEnergyVar[j];
02325           adaptTempRead >> adaptTempPotEnergySamples[j];
02326           adaptTempRead >> adaptTempPotEnergyAveNum[j];
02327           adaptTempRead >> adaptTempPotEnergyVarNum[j];
02328           adaptTempRead >> adaptTempPotEnergyAveDen[j];
02329           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02330         } 
02331         adaptTempRead.close();
02332       }
02333       else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
02334     } 
02335     else {
02336       adaptTempBins = simParams->adaptTempBins;
02337       adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02338       adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02339       adaptTempPotEnergySamples = new int[adaptTempBins];
02340       adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02341       adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02342       adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02343       adaptTempBetaN           = new BigReal[adaptTempBins];
02344       adaptTempBetaMax = 1./simParams->adaptTempTmin;
02345       adaptTempBetaMin = 1./simParams->adaptTempTmax;
02346       adaptTempCg = simParams->adaptTempCgamma;   
02347       adaptTempDt  = simParams->adaptTempDt;
02348       adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02349       adaptTempT = simParams->initialTemp; 
02350       if (simParams->langevinOn)
02351         adaptTempT = simParams->langevinTemp;
02352       else if (simParams->rescaleFreq > 0)
02353         adaptTempT = simParams->rescaleTemp;
02354       for(int j = 0; j < adaptTempBins; ++j){
02355           adaptTempPotEnergyAveNum[j] = 0.;
02356           adaptTempPotEnergyAveDen[j] = 0.;
02357           adaptTempPotEnergySamples[j] = 0;
02358           adaptTempPotEnergyVarNum[j] = 0.;
02359           adaptTempPotEnergyVar[j] = 0.;
02360           adaptTempPotEnergyAve[j] = 0.;
02361           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02362       }
02363     }
02364     if (simParams->adaptTempAutoDt > 0.0) {
02365        adaptTempAutoDt = true;
02366        adaptTempDtMin =  simParams->adaptTempAutoDt - 0.01;
02367        adaptTempDtMax =  simParams->adaptTempAutoDt + 0.01;
02368     }
02369     adaptTempDTave = 0;
02370     adaptTempDTavenum = 0;
02371     iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
02372     iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
02373     iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
02374     iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
02375     iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
02376     if (simParams->adaptTempRestartFreq > 0) {
02377       NAMD_backup_file(simParams->adaptTempRestartFile);
02378       adaptTempRestartFile.open(simParams->adaptTempRestartFile);
02379        if(!adaptTempRestartFile)
02380         NAMD_die("Error opening restart file");
02381     }
02382 }

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

Definition at line 2410 of file Controller.C.

References adaptTempBetaMax, adaptTempBetaMin, adaptTempBetaN, adaptTempBin, adaptTempDBeta, SimParameters::adaptTempDebug, SimParameters::adaptTempFreq, adaptTempInit(), SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, SimParameters::adaptTempOutFreq, adaptTempPotEnergyAve, adaptTempPotEnergyAveDen, adaptTempPotEnergyAveNum, adaptTempPotEnergySamples, adaptTempPotEnergyVar, adaptTempPotEnergyVarNum, endi(), SimParameters::firstTimestep, iout, iWARN(), j, kineticEnergy, simParams, and totalEnergy.

Referenced by integrate().

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

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

Definition at line 2384 of file Controller.C.

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

02384                                                {
02385     if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
02386         adaptTempRestartFile.seekp(std::ios::beg);        
02387         iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
02388         adaptTempRestartFile << step << " ";
02389         // Start with min and max temperatures
02390         adaptTempRestartFile << adaptTempT<< " ";     // KELVIN
02391         adaptTempRestartFile << 1./adaptTempBetaMin << " ";  // KELVIN
02392         adaptTempRestartFile << 1./adaptTempBetaMax << " ";  // KELVIN
02393         adaptTempRestartFile << adaptTempBins << " ";     
02394         adaptTempRestartFile << adaptTempCg << " ";
02395         adaptTempRestartFile << adaptTempDt ;
02396         adaptTempRestartFile << "\n" ;
02397         for(int j = 0; j < adaptTempBins; ++j) {
02398           adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
02399           adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
02400           adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
02401           adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
02402           adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
02403           adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
02404           adaptTempRestartFile << "\n";          
02405         }
02406         adaptTempRestartFile.flush(); 
02407     }
02408 }    

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

Definition at line 280 of file Controller.C.

References BackEnd::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(), SimParameters::initialTemp, integrate(), iout, Controller::checkpoint::lattice, NamdState::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().

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

void Controller::awaken ( void   )  [inline]

Definition at line 45 of file Controller.h.

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

00045 { CthAwaken(thread); };

void Controller::berendsenPressure ( int   )  [protected]

Definition at line 949 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, NamdState::lattice, LIMIT_SCALING, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), Lattice::rescale(), simParams, state, Tensor::xx, Tensor::yy, and Tensor::zz.

Referenced by integrate().

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

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 1714 of file Controller.C.

Referenced by rescaleaccelMD().

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

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 1743 of file Controller.C.

Referenced by rescaleaccelMD().

01744                                            {
01745    BigReal VE = testV - E;
01746    //if V < E, apply boost
01747    if(VE < 0){
01748       *factor = k * VE;
01749       *vir += vir_orig * (*factor);
01750       *dV += (*factor) * VE/2.;
01751       *factor += 1.;
01752    }
01753    else{
01754       *factor = 1.;
01755    }
01756 }

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 1694 of file Controller.C.

Referenced by rescaleaccelMD().

01695                                                                            {
01696    BigReal delta; 
01697 
01698    if(testV > *Vmax){
01699       *Vmax = testV;
01700    }
01701    else if(testV < *Vmin){
01702       *Vmin = testV;
01703    }
01704 
01705    //mean and std calculated by Online Algorithm
01706    delta = testV - *Vavg;
01707    *Vavg += delta / (BigReal)n;
01708    *M2 += delta * (testV - (*Vavg));
01709 
01710    *sigmaV = sqrt(*M2/n);
01711 }

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 1576 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(), NamdState::lattice, 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().

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

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

Definition at line 2708 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().

02708                                                          {
02709     Node *node = Node::Object();
02710     Molecule *molecule = node->molecule;
02711 
02712     // Some basic correctness checking
02713     BigReal checksum, checksum_b;
02714     char errmsg[256];
02715 
02716     checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
02717     if ( ((int)checksum) != molecule->numAtoms ) {
02718       sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
02719               (int)checksum, molecule->numAtoms);
02720       NAMD_bug(errmsg);
02721     }
02722 
02723     checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
02724     if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
02725       if ( ! computeChecksum ) {
02726         computesPartitioned = 0;
02727       } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
02728         computesPartitioned = 1;
02729       } else {
02730         NAMD_bug("Bad global compute count!\n");
02731       }
02732       computeChecksum = ((int)checksum);
02733     }
02734 
02735     checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
02736     if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
02737       sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
02738               (int)checksum, molecule->numCalcBonds);
02739       if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
02740         iout << iWARN << errmsg << endi;
02741       else NAMD_bug(errmsg);
02742     }
02743 
02744     checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
02745     if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
02746       sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
02747               (int)checksum, molecule->numCalcAngles);
02748       if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
02749         iout << iWARN << errmsg << endi;
02750       else NAMD_bug(errmsg);
02751     }
02752 
02753     checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
02754     if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
02755       sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
02756               (int)checksum, molecule->numCalcDihedrals);
02757       if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
02758         iout << iWARN << errmsg << endi;
02759       else NAMD_bug(errmsg);
02760     }
02761 
02762     checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
02763     if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
02764       sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
02765               (int)checksum, molecule->numCalcImpropers);
02766       if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
02767         iout << iWARN << errmsg << endi;
02768       else NAMD_bug(errmsg);
02769     }
02770 
02771     checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
02772     if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
02773       sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
02774               (int)checksum, molecule->numCalcTholes);
02775       if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
02776         iout << iWARN << errmsg << endi;
02777       else NAMD_bug(errmsg);
02778     }
02779 
02780     checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
02781     if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
02782       sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
02783               (int)checksum, molecule->numCalcAnisos);
02784       if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
02785         iout << iWARN << errmsg << endi;
02786       else NAMD_bug(errmsg);
02787     }
02788 
02789     checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
02790     if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
02791       sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
02792               (int)checksum, molecule->numCalcCrossterms);
02793       if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
02794         iout << iWARN << errmsg << endi;
02795       else NAMD_bug(errmsg);
02796     }
02797 
02798     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
02799     if ( ((int64)checksum) > molecule->numCalcExclusions &&
02800          ( ! simParams->mollyOn || step % slowFreq ) ) {
02801       if ( forgiving )
02802         iout << iWARN << "High global exclusion count ("
02803                       << ((int64)checksum) << " vs "
02804                       << molecule->numCalcExclusions << "), possible error!\n"
02805           << iWARN << "This warning is not unusual during minimization.\n"
02806           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02807       else {
02808         char errmsg[256];
02809         sprintf(errmsg, "High global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02810                 (long long)checksum, (long long)molecule->numCalcExclusions);
02811         NAMD_bug(errmsg);
02812       }
02813     }
02814 
02815     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
02816         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02817       if ( forgiving || ! simParams->fullElectFrequency ) {
02818         iout << iWARN << "Low global exclusion count!  ("
02819           << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
02820         if ( forgiving ) {
02821           iout << "\n"
02822             << iWARN << "This warning is not unusual during minimization.\n"
02823             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02824         } else {
02825           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02826         }
02827         iout << endi;
02828       } else {
02829         char errmsg[256];
02830         sprintf(errmsg, "Low global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
02831                 (long long)checksum, (long long)molecule->numCalcExclusions);
02832         NAMD_bug(errmsg);
02833       }
02834     }
02835 
02836 #ifdef NAMD_CUDA
02837     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
02838     if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
02839          ( ! simParams->mollyOn || step % slowFreq ) ) {
02840       if ( forgiving )
02841         iout << iWARN << "High global CUDA exclusion count ("
02842                       << ((int64)checksum) << " vs "
02843                       << molecule->numCalcFullExclusions << "), possible error!\n"
02844           << iWARN << "This warning is not unusual during minimization.\n"
02845           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02846       else {
02847         char errmsg[256];
02848         sprintf(errmsg, "High global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02849                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02850         NAMD_bug(errmsg);
02851       }
02852     }
02853 
02854     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
02855         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02856       if ( forgiving || ! simParams->fullElectFrequency ) {
02857         iout << iWARN << "Low global CUDA exclusion count!  ("
02858           << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
02859         if ( forgiving ) {
02860           iout << "\n"
02861             << iWARN << "This warning is not unusual during minimization.\n"
02862             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02863         } else {
02864           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02865         }
02866         iout << endi;
02867       } else {
02868         char errmsg[256];
02869         sprintf(errmsg, "Low global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
02870                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
02871         NAMD_bug(errmsg);
02872       }
02873     }
02874 #endif
02875 
02876     checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
02877     if ( ((int)checksum) && ! marginViolations ) {
02878       iout << iERROR << "Margin is too small for " << ((int)checksum) <<
02879         " atoms during timestep " << step << ".\n" << iERROR <<
02880       "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
02881     }
02882     marginViolations += (int)checksum;
02883 
02884     checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
02885     if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
02886       iout << iINFO <<
02887         "Pairlistdist is too small for " << ((int)checksum) <<
02888         " computes during timestep " << step << ".\n" << endi;
02889     }
02890     if ( simParams->outputPairlists )  pairlistWarnings += (int)checksum;
02891 
02892     checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
02893     if ( checksum ) {
02894       if ( forgiving )
02895         iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
02896       else NAMD_bug("Stray PME grid charges detected!\n");
02897     }
02898 }

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

Definition at line 3787 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().

03787                                                   {
03788   // alchemical scaling factors for groups 1/2 at the previous lambda
03789   const BigReal oldLambda = simParams->getCurrentLambda(step-1);
03790   const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
03791   const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
03792   const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
03793   const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
03794   const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
03795   const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
03796   // alchemical scaling factors for groups 1/2 at the new lambda
03797   const BigReal alchLambda = simParams->getCurrentLambda(step);
03798   const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
03799   const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
03800   const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
03801   const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
03802   const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
03803   const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda); 
03804 
03805   return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
03806           (elec_lambda_12 - elec_lambda_1)*
03807           (electEnergy_ti_1 + electEnergySlow_ti_1 +  electEnergyPME_ti_1) +
03808           (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
03809           (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
03810           (elec_lambda_22 - elec_lambda_2)*
03811           (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) + 
03812           (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
03813   );
03814 }

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

Definition at line 1254 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().

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

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

Definition at line 4087 of file Controller.C.

References broadcast.

Referenced by integrate().

04087                                                      {
04088 #if USE_BARRIER
04089         if (doBarrier) {
04090           broadcast->cycleBarrier.publish(step,1);
04091           CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
04092                   CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04093         }
04094 #endif
04095 }

void Controller::enqueueCollections ( int   )  [protected]

Definition at line 3572 of file Controller.C.

References collection, Output::coordinateNeeded(), CollectionMaster::enqueueForces(), CollectionMaster::enqueuePositions(), CollectionMaster::enqueueVelocities(), Output::forceNeeded(), NamdState::lattice, state, and Output::velocityNeeded().

Referenced by algorithm(), integrate(), and minimize().

03573 {
03574   if ( Output::coordinateNeeded(timestep) ){
03575     collection->enqueuePositions(timestep,state->lattice);
03576   }
03577   if ( Output::velocityNeeded(timestep) )
03578     collection->enqueueVelocities(timestep);
03579   if ( Output::forceNeeded(timestep) )
03580     collection->enqueueForces(timestep);
03581 }

void Controller::integrate ( int   )  [protected]

Definition at line 428 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().

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

void Controller::langevinPiston1 ( int   )  [protected]

Definition at line 986 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, NamdState::lattice, 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().

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

void Controller::langevinPiston2 ( int   )  [protected]

Definition at line 1141 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, NamdState::lattice, 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().

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

void Controller::minimize (  )  [protected]

Definition at line 593 of file Controller.C.

References broadcast, CALCULATE, minpoint::dudx, endi(), enqueueCollections(), SimParameters::firstTimestep, if(), 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, RequireReduction::require(), simParams, slowFreq, minpoint::u, minpoint::x, and x.

Referenced by algorithm().

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

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

Definition at line 918 of file Controller.C.

References BOLTZMANN, controlNumDegFreedom, kineticEnergy, NamdState::lattice, 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().

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

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

Definition at line 777 of file Controller.C.

References BOLTZMANN, broadcast, calcPressure(), controlNumDegFreedom, controlPressure, SimParameters::dt, GET_TENSOR, GET_VECTOR, Tensor::identity(), RequireReduction::item(), kineticEnergy, NamdState::lattice, 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, REDUCTION_EXT_FORCE_NBOND, REDUCTION_EXT_FORCE_NORMAL, REDUCTION_EXT_FORCE_SLOW, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_MOMENTUM_SQUARED, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, RequireReduction::require(), Lattice::rescale(), simParams, state, temperature, ControllerBroadcasts::velocityRescaleTensor, ControllerBroadcasts::velocityRescaleTensor2, and Lattice::volume().

Referenced by integrate().

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

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

Definition at line 852 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 3915 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().

03916 {
03917 
03918   if ( step >= 0 ) {
03919 
03920     // Write out eXtended System Trajectory (XST) file
03921     if ( simParams->xstFrequency &&
03922          ((step % simParams->xstFrequency) == 0) )
03923     {
03924       if ( ! xstFile.is_open() )
03925       {
03926         iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
03927         NAMD_backup_file(simParams->xstFilename);
03928         xstFile.open(simParams->xstFilename);
03929         while (!xstFile) {
03930           if ( errno == EINTR ) {
03931             CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
03932             xstFile.clear();
03933             xstFile.open(simParams->xstFilename);
03934             continue;
03935           }
03936           char err_msg[257];
03937           sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
03938           NAMD_err(err_msg);
03939         }
03940         xstFile << "# NAMD extended system trajectory file" << std::endl;
03941         writeExtendedSystemLabels(xstFile);
03942       }
03943       writeExtendedSystemData(step,xstFile);
03944       xstFile.flush();
03945     }
03946 
03947     // Write out eXtended System Configuration (XSC) files
03948     //  Output a restart file
03949     if ( simParams->restartFrequency &&
03950          ((step % simParams->restartFrequency) == 0) &&
03951          (step != simParams->firstTimestep) )
03952     {
03953       iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
03954                 << step << "\n" << endi;
03955       char fname[140];
03956       const char *bsuffix = ".old";
03957       strcpy(fname, simParams->restartFilename);
03958       if ( simParams->restartSave ) {
03959         char timestepstr[20];
03960         sprintf(timestepstr,".%d",step);
03961         strcat(fname, timestepstr);
03962         bsuffix = ".BAK";
03963       }
03964       strcat(fname, ".xsc");
03965       NAMD_backup_file(fname,bsuffix);
03966       ofstream_namd xscFile(fname);
03967       while (!xscFile) {
03968         if ( errno == EINTR ) {
03969           CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
03970           xscFile.clear();
03971           xscFile.open(fname);
03972           continue;
03973         }
03974         char err_msg[257];
03975         sprintf(err_msg, "Error opening XSC restart file %s",fname);
03976         NAMD_err(err_msg);
03977       } 
03978       xscFile << "# NAMD extended system configuration restart file" << std::endl;
03979       writeExtendedSystemLabels(xscFile);
03980       writeExtendedSystemData(step,xscFile);
03981       if (!xscFile) {
03982         char err_msg[257];
03983         sprintf(err_msg, "Error writing XSC restart file %s",fname);
03984         NAMD_err(err_msg);
03985       } 
03986     }
03987 
03988   }
03989 
03990   //  Output final coordinates
03991   if (step == FILE_OUTPUT || step == END_OF_RUN)
03992   {
03993     int realstep = ( step == FILE_OUTPUT ?
03994         simParams->firstTimestep : simParams->N );
03995     iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
03996                 << realstep << "\n" << endi;
03997     static char fname[140];
03998     strcpy(fname, simParams->outputFilename);
03999     strcat(fname, ".xsc");
04000     NAMD_backup_file(fname);
04001     ofstream_namd xscFile(fname);
04002     while (!xscFile) {
04003       if ( errno == EINTR ) {
04004         CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
04005         xscFile.clear();
04006         xscFile.open(fname);
04007         continue;
04008       }
04009       char err_msg[257];
04010       sprintf(err_msg, "Error opening XSC output file %s",fname);
04011       NAMD_err(err_msg);
04012     } 
04013     xscFile << "# NAMD extended system configuration output file" << std::endl;
04014     writeExtendedSystemLabels(xscFile);
04015     writeExtendedSystemData(realstep,xscFile);
04016     if (!xscFile) {
04017       char err_msg[257];
04018       sprintf(err_msg, "Error writing XSC output file %s",fname);
04019       NAMD_err(err_msg);
04020     } 
04021   }
04022 
04023   //  Close trajectory file
04024   if (step == END_OF_RUN) {
04025     if ( xstFile.is_open() ) {
04026       xstFile.close();
04027       iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
04028     }
04029   }
04030 
04031 }

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

Definition at line 3584 of file Controller.C.

References SimParameters::alchEnsembleAvg, SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchFepWhamOn, SimParameters::alchLambda, SimParameters::alchLambda2, 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().

03584                                          {
03585  if (simParams->alchOn && simParams->alchFepOn) {
03586   const int stepInRun = step - simParams->firstTimestep;
03587   const int alchEquilSteps = simParams->alchEquilSteps;
03588   const BigReal alchLambda = simParams->alchLambda;
03589   const BigReal alchLambda2 = simParams->alchLambda2;
03590   const bool alchEnsembleAvg = simParams->alchEnsembleAvg; 
03591   const bool FepWhamOn = simParams->alchFepWhamOn;
03592 
03593   if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
03594     FepNo = 0;
03595     exp_dE_ByRT = 0.0;
03596     net_dE = 0.0;
03597   }
03598   dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
03599        electEnergy - electEnergySlow - ljEnergy;
03600   BigReal RT = BOLTZMANN * simParams->alchTemp;
03601 
03602   if (alchEnsembleAvg){
03603     FepNo++;
03604     exp_dE_ByRT += exp(-dE/RT);
03605     net_dE += dE;
03606   }
03607 
03608   if (simParams->alchOutFreq) { 
03609     if (stepInRun == 0) {
03610       if (!fepFile.is_open()) {
03611         NAMD_backup_file(simParams->alchOutFile);
03612         fepFile.open(simParams->alchOutFile);
03613         iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
03614         if(alchEnsembleAvg){
03615           fepSum = 0.0;
03616           fepFile << "#            STEP                 Elec                            "
03617                   << "vdW                    dE           dE_avg         Temp             dG\n"
03618                   << "#                           l             l+dl      "
03619                   << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03620         }
03621         else{
03622           if(!FepWhamOn){ 
03623             fepFile << "#            STEP                 Elec                            "
03624                     << "vdW                    dE         Temp\n"
03625                     << "#                           l             l+dl      "
03626                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03627           } 
03628         }
03629       }
03630       if(!step){
03631         if(!FepWhamOn){ 
03632           fepFile << "#NEW FEP WINDOW: "
03633                   << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
03634                   << alchLambda2 << std::endl;
03635         }
03636       }
03637     }
03638     if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03639       fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03640               << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03641               << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03642     }
03643     if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
03644       writeFepEnergyData(step, fepFile);
03645       fepFile.flush();
03646     }
03647     if (alchEnsembleAvg && (step == simParams->N)) {
03648       fepSum = fepSum + dG;
03649       fepFile << "#Free energy change for lambda window [ " << alchLambda
03650               << " " << alchLambda2 << " ] is " << dG 
03651               << " ; net change until now is " << fepSum << std::endl;
03652       fepFile.flush();
03653     }
03654   }
03655  }
03656 }

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

Definition at line 3658 of file Controller.C.

References SimParameters::alchEquilSteps, SimParameters::alchLambda, SimParameters::alchLambda2, 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::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().

03658                                         {
03659  if (simParams->alchOn && simParams->alchThermIntOn) {
03660   const int stepInRun = step - simParams->firstTimestep;
03661   const int alchEquilSteps = simParams->alchEquilSteps;
03662   const int alchLambdaFreq = simParams->alchLambdaFreq;
03663 
03664   if (stepInRun == 0 || stepInRun == alchEquilSteps) {
03665     TiNo = 0;
03666     net_dEdl_bond_1 = 0;
03667     net_dEdl_bond_2 = 0;
03668     net_dEdl_elec_1 = 0;
03669     net_dEdl_elec_2 = 0;
03670     net_dEdl_lj_1 = 0;
03671     net_dEdl_lj_2 = 0;
03672   }
03673   TiNo++;
03674   net_dEdl_bond_1 += bondedEnergy_ti_1;
03675   net_dEdl_bond_2 += bondedEnergy_ti_2;
03676   net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
03677                       + electEnergyPME_ti_1);
03678   net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
03679                       + electEnergyPME_ti_2);
03680   net_dEdl_lj_1 += ljEnergy_ti_1;
03681   net_dEdl_lj_2 += ljEnergy_ti_2;
03682 
03683   // Don't accumulate block averages (you'll get a divide by zero!) or even 
03684   // open the TI output if alchOutFreq is zero.
03685   if (simParams->alchOutFreq) {
03686     if (stepInRun == 0 || stepInRun == alchEquilSteps 
03687         || (! ((step - 1) % simParams->alchOutFreq))) {
03688       // output of instantaneous dU/dl now replaced with running average
03689       // over last alchOutFreq steps (except for step 0)
03690       recent_TiNo = 0;
03691       recent_dEdl_bond_1 = 0;
03692       recent_dEdl_bond_2 = 0;
03693       recent_dEdl_elec_1 = 0;
03694       recent_dEdl_elec_2 = 0;
03695       recent_dEdl_lj_1 = 0;
03696       recent_dEdl_lj_2 = 0;
03697       recent_alchWork = 0;
03698     }
03699     recent_TiNo++;
03700     recent_dEdl_bond_1 += bondedEnergy_ti_1;
03701     recent_dEdl_bond_2 += bondedEnergy_ti_2;
03702     recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 
03703                            + electEnergyPME_ti_1); 
03704     recent_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 
03705                            + electEnergyPME_ti_2); 
03706     recent_dEdl_lj_1 += ljEnergy_ti_1;
03707     recent_dEdl_lj_2 += ljEnergy_ti_2;
03708     recent_alchWork += alchWork;
03709 
03710     if (stepInRun == 0) {
03711       if (!tiFile.is_open()) {
03712         NAMD_backup_file(simParams->alchOutFile);
03713         tiFile.open(simParams->alchOutFile);
03714         /* BKR - This has been rather drastically updated to better match 
03715            stdout. This was necessary for several reasons:
03716            1) PME global scaling is obsolete (now removed)
03717            2) scaling of bonded terms was added
03718            3) alchemical work is now accumulated when switching is active
03719          */
03720         iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
03721         tiFile << "#TITITLE:    TS";
03722         tiFile << FORMAT("BOND1");
03723         tiFile << FORMAT("AVGBOND1");
03724         tiFile << FORMAT("ELECT1");
03725         tiFile << FORMAT("AVGELECT1");
03726         tiFile << "     ";
03727         tiFile << FORMAT("VDW1");
03728         tiFile << FORMAT("AVGVDW1");
03729         tiFile << FORMAT("BOND2");
03730         tiFile << FORMAT("AVGBOND2");
03731         tiFile << FORMAT("ELECT2");
03732         tiFile << "     ";
03733         tiFile << FORMAT("AVGELECT2");
03734         tiFile << FORMAT("VDW2");
03735         tiFile << FORMAT("AVGVDW2");
03736         if (alchLambdaFreq > 0) {
03737           tiFile << FORMAT("ALCHWORK");
03738           tiFile << FORMAT("CUMALCHWORK");
03739         }
03740         tiFile << std::endl;
03741       }
03742 
03743       if (alchLambdaFreq > 0) {
03744         tiFile << "#ALCHEMICAL SWITCHING ACTIVE " 
03745                << simParams->alchLambda << " --> " << simParams->alchLambda2
03746                << "\n#LAMBDA SCHEDULE: " 
03747                << "dL: " << simParams->getLambdaDelta() 
03748                << " Freq: " << alchLambdaFreq;
03749       }
03750       else {
03751         const BigReal alchLambda = simParams->alchLambda;    
03752         const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
03753         const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
03754         const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
03755         const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
03756         const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
03757         const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
03758         tiFile << "#NEW TI WINDOW: "
03759                << "LAMBDA " << alchLambda 
03760                << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
03761                << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
03762                << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
03763                << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
03764       }
03765       tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
03766              << std::endl;
03767     }
03768 
03769     if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03770       tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03771              << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03772              << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03773     }
03774     if ((step%simParams->alchOutFreq) == 0) {
03775       writeTiEnergyData(step, tiFile);
03776       tiFile.flush();
03777     }
03778   }
03779  }
03780 }

void Controller::printDynamicsEnergies ( int   )  [protected]

Definition at line 2956 of file Controller.C.

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

Referenced by integrate().

02956                                                {
02957 
02958     Node *node = Node::Object();
02959     Molecule *molecule = node->molecule;
02960     SimParameters *simParameters = node->simParameters;
02961     Lattice &lattice = state->lattice;
02962 
02963     compareChecksums(step);
02964 
02965     printEnergies(step,0);
02966 }

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

Definition at line 2968 of file Controller.C.

References Lattice::a(), Lattice::a_p(), SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchLambda2, SimParameters::alchLambdaFreq, SimParameters::alchOn, SimParameters::alchThermIntOn, SimParameters::alchVdwLambdaEnd, alchWork, avg_count, Lattice::b(), Lattice::b_p(), bondedEnergy_ti_1, bondedEnergy_ti_2, bondedEnergyDiff_f, Lattice::c(), Lattice::c_p(), NamdState::callback_labelstring, NamdState::callback_valuestring, CALLBACKDATA, CALLBACKLIST, computeAlchWork(), cumAlchWork, ScriptTcl::doCallback(), 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(), PressureProfileReduction::getData(), Molecule::getEnergyTailCorr(), Node::getScript(), SimParameters::goForcesOn, SimParameters::goGroPair, goNativeEnergy, goNonnativeEnergy, goTotalEnergy, groGaussEnergy, groLJEnergy, groupPressure, groupPressure_avg, groupPressure_tavg, iERROR(), iINFO(), Node::imd, SimParameters::IMDfreq, SimParameters::IMDon, iout, RequireReduction::item(), iWARN(), kineticEnergy, kineticEnergyCentered, kineticEnergyHalfstep, NamdState::lattice, SimParameters::LJcorrection, ljEnergy, ljEnergy_f, ljEnergy_f_left, 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_F_LEFT, REDUCTION_LJ_ENERGY_TI_1, REDUCTION_LJ_ENERGY_TI_2, REDUCTION_MISC_ENERGY, REDUCTION_PAIR_ELECT_FORCE, REDUCTION_PAIR_VDW_FORCE, Node::simParameters, simParams, slowFreq, ControllerState::smooth2_avg, state, stepInFullRun, IMDEnergies::T, Molecule::tail_corr_dUdl_1, Molecule::tail_corr_dUdl_2, tavg_count, temp_avg, temperature, TITITLE(), totalEnergy, IMDEnergies::tstep, Lattice::volume(), Vector::x, XXXBIGREAL, Vector::y, and Vector::z.

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

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

void Controller::printFepMessage ( int   )  [protected]

Definition at line 1273 of file Controller.C.

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

Referenced by integrate().

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

void Controller::printMinimizeEnergies ( int   )  [protected]

Definition at line 2937 of file Controller.C.

References compareChecksums(), RequireReduction::item(), min_energy, min_f_dot_v, min_huge_count, min_v_dot_v, Node::molecule, Node::Object(), printEnergies(), receivePressure(), reduction, REDUCTION_MIN_F_DOT_F, REDUCTION_MIN_F_DOT_V, REDUCTION_MIN_HUGE_COUNT, REDUCTION_MIN_V_DOT_V, rescaleaccelMD(), and totalEnergy.

02937                                                {
02938 
02939     rescaleaccelMD(step,1);
02940     receivePressure(step,1);
02941 
02942     Node *node = Node::Object();
02943     Molecule *molecule = node->molecule;
02944     compareChecksums(step,1);
02945 
02946     printEnergies(step,1);
02947 
02948     min_energy = totalEnergy;
02949     min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
02950     min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
02951     min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
02952     min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
02953 
02954 }

void Controller::printTiMessage ( int   )  [protected]

Definition at line 1289 of file Controller.C.

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

Referenced by integrate().

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

void Controller::printTiming ( int   )  [protected]

Definition at line 2900 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().

02900                                      {
02901 
02902     if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
02903     {
02904       const double endWTime = CmiWallTimer() - firstWTime;
02905       const double endCTime = CmiTimer() - firstCTime;
02906 
02907       // fflush about once per minute
02908       if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
02909 
02910       const double elapsedW = 
02911         (endWTime - startWTime) / simParams->outputTiming;
02912       const double elapsedC = 
02913         (endCTime - startCTime) / simParams->outputTiming;
02914 
02915       const double remainingW = elapsedW * (simParams->N - step);
02916       const double remainingW_hours = remainingW / 3600;
02917 
02918       startWTime = endWTime;
02919       startCTime = endCTime;
02920 
02921 #ifdef NAMD_CUDA
02922       if ( simParams->outputEnergies < 60 &&
02923            step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
02924         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
02925       }
02926 #endif
02927       if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
02928         CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
02929                   ", %g hours remaining, %f MB of memory in use.\n",
02930                   step, endCTime, elapsedC, endWTime, elapsedW,
02931                   remainingW_hours, memusage_MB());
02932         if ( fflush_count ) { --fflush_count; fflush(stdout); }
02933       }
02934     }
02935 }

void Controller::reassignVelocities ( int   )  [protected]

Definition at line 1303 of file Controller.C.

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

Referenced by integrate().

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

void Controller::rebalanceLoad ( int   )  [protected]

Definition at line 4073 of file Controller.C.

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

Referenced by integrate().

04074 {
04075   if ( ! ldbSteps ) { 
04076     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
04077   }
04078   if ( ! --ldbSteps ) {
04079     startBenchTime -= CmiWallTimer();
04080         Node::Object()->outputPatchComputeMaps("before_ldb", step);
04081     LdbCoordinator::Object()->rebalance(this);  
04082         startBenchTime += CmiWallTimer();
04083     fflush_count = 3;
04084   }
04085 }

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

Definition at line 1398 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, NamdState::lattice, 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_EXT_FORCE_NBOND, REDUCTION_EXT_FORCE_NORMAL, REDUCTION_EXT_FORCE_SLOW, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_MOMENTUM_SQUARED, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, RequireReduction::require(), Node::simParameters, simParams, state, temperature, and SimParameters::useGroupPressure.

Referenced by integrate(), and printMinimizeEnergies().

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

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

Definition at line 4064 of file Controller.C.

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

Referenced by Node::recvCheckpointAck().

04064                                                  {  // initiating replica
04065   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
04066     state->lattice = cp.lattice;
04067     *(ControllerState*)this = cp.state;
04068   }
04069   CkpvAccess(_qd)->process();
04070 }

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

Definition at line 4034 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().

04034                                                                             {  // responding replica
04035   switch ( task ) {
04036     case SCRIPT_CHECKPOINT_STORE:
04037       if ( ! checkpoints.count(key) ) {
04038         checkpoints[key] = new checkpoint;
04039       }
04040       *checkpoints[key] = cp;
04041       break;
04042     case SCRIPT_CHECKPOINT_LOAD:
04043       if ( ! checkpoints.count(key) ) {
04044         NAMD_die("Unable to load checkpoint, requested key was never stored.");
04045       }
04046       cp = *checkpoints[key];
04047       break;
04048     case SCRIPT_CHECKPOINT_SWAP:
04049       if ( ! checkpoints.count(key) ) {
04050         NAMD_die("Unable to swap checkpoint, requested key was never stored.");
04051       }
04052       std::swap(cp,*checkpoints[key]);
04053       break;
04054     case SCRIPT_CHECKPOINT_FREE:
04055       if ( ! checkpoints.count(key) ) {
04056         NAMD_die("Unable to free checkpoint, requested key was never stored.");
04057       }
04058       delete checkpoints[key];
04059       checkpoints.erase(key);
04060       break;
04061   }
04062 }

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

Definition at line 1817 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::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(), NamdState::lattice, 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, REDUCTION_VIRIAL_AMD_DIHE, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, 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().

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

void Controller::rescaleVelocities ( int   )  [protected]

Definition at line 1232 of file Controller.C.

References SimParameters::adaptTempFirstStep, SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, SimParameters::adaptTempRescale, adaptTempT, broadcast, SimParameters::rescaleFreq, SimParameters::rescaleTemp, rescaleVelocities_numTemps, rescaleVelocities_sumTemps, simParams, temperature, and ControllerBroadcasts::velocityRescaleFactor.

Referenced by integrate().

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

void Controller::resumeAfterTraceBarrier ( int   ) 

Definition at line 4104 of file Controller.C.

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

Referenced by Node::resumeAfterTraceBarrier().

04104                                                 {
04105         broadcast->traceBarrier.publish(step,1);
04106         CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime); 
04107         awaken();
04108 }

void Controller::run ( void   ) 

Definition at line 267 of file Controller.C.

References awaken(), CTRL_STK_SZ, and DebugM.

Referenced by NamdState::runController().

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

void Controller::stochRescaleVelocities ( int   )  [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.

Definition at line 1343 of file Controller.C.

References broadcast, Random::gaussian(), numDegFreedom, random, simParams, stochRescale_count, ControllerBroadcasts::stochRescaleCoefficient, SimParameters::stochRescaleFreq, SimParameters::stochRescaleOn, SimParameters::stochRescaleTemp, stochRescaleTimefactor, Random::sum_of_squared_gaussians(), and temperature.

Referenced by integrate().

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

void Controller::tcoupleVelocities ( int   )  [protected]

Definition at line 1321 of file Controller.C.

References broadcast, simParams, ControllerBroadcasts::tcoupleCoefficient, SimParameters::tCoupleOn, SimParameters::tCoupleTemp, and temperature.

Referenced by integrate().

01322 {
01323   if ( simParams->tCoupleOn )
01324   {
01325     const BigReal tCoupleTemp = simParams->tCoupleTemp;
01326     BigReal coefficient = 1.;
01327     if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
01328     broadcast->tcoupleCoefficient.publish(step,coefficient);
01329   }
01330 }

void Controller::terminate ( void   )  [protected]

Definition at line 4125 of file Controller.C.

References BackEnd::awaken().

Referenced by algorithm().

04125                                {
04126   BackEnd::awaken();
04127   CthFree(thread);
04128   CthSuspend();
04129 }

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

Definition at line 4097 of file Controller.C.

Referenced by integrate().

04097                                                        {
04098         CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);       
04099         CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04100         nd.traceBarrier(turnOnTrace, step);
04101         CthSuspend();
04102 }

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 1759 of file Controller.C.

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

Referenced by rescaleaccelMD().

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

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

Definition at line 3551 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, NamdState::lattice, Lattice::origin(), simParams, state, Vector::x, Vector::y, and Vector::z.

Referenced by outputExtendedSystem().

03551                                                                       {
03552   Lattice &lattice = state->lattice;
03553   file.precision(12);
03554   file << step;
03555     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
03556     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
03557     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
03558     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
03559   if ( simParams->langevinPistonOn ) {
03560     Vector strainRate = diagonal(langevinPiston_strainRate);
03561     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
03562     file << " " << strainRate.x;
03563     file << " " << strainRate.y;
03564     file << " " << strainRate.z;
03565     file << " " << strainRate2.x;
03566     file << " " << strainRate2.y;
03567     file << " " << strainRate2.z;
03568   }
03569   file << std::endl;
03570 }

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

Definition at line 3538 of file Controller.C.

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

Referenced by outputExtendedSystem().

03538                                                               {
03539   Lattice &lattice = state->lattice;
03540   file << "#$LABELS step";
03541   if ( lattice.a_p() ) file << " a_x a_y a_z";
03542   if ( lattice.b_p() ) file << " b_x b_y b_z";
03543   if ( lattice.c_p() ) file << " c_x c_y c_z";
03544   file << " o_x o_y o_z";
03545   if ( simParams->langevinPistonOn ) {
03546     file << " s_x s_y s_z s_u s_v s_w";
03547   }
03548   file << std::endl;
03549 }

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

Definition at line 3816 of file Controller.C.

References SimParameters::alchDispLambda, SimParameters::alchElecLambda, SimParameters::alchEnsembleAvg, SimParameters::alchFepElecOn, SimParameters::alchFepWCADispOn, SimParameters::alchFepWCArcut1, SimParameters::alchFepWCArcut2, SimParameters::alchFepWCArcut3, SimParameters::alchFepWCARepuOn, SimParameters::alchFepWhamOn, SimParameters::alchLambda, SimParameters::alchRepLambda, SimParameters::alchTemp, BOLTZMANN, dE, dG, electEnergy, electEnergy_f, electEnergySlow, electEnergySlow_f, exp_dE_ByRT, fepFile, FepNo, FEPTITLE(), SimParameters::firstTimestep, FORMAT(), ljEnergy, ljEnergy_f, ljEnergy_f_left, net_dE, simParams, and temperature.

Referenced by outputFepEnergy().

03816                                                                  {
03817   BigReal eeng = electEnergy + electEnergySlow;
03818   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
03819   BigReal dE_Left = eeng_f + ljEnergy_f_left - eeng - ljEnergy;
03820   BigReal RT = BOLTZMANN * simParams->alchTemp;
03821 
03822   const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03823   const int stepInRun = step - simParams->firstTimestep;
03824   const bool FepWhamOn = simParams->alchFepWhamOn;
03825   const bool WCARepuOn = simParams->alchFepWCARepuOn;
03826   const BigReal WCArcut1 = simParams->alchFepWCArcut1;
03827   const BigReal WCArcut2 = simParams->alchFepWCArcut2;
03828   const BigReal WCArcut3 = simParams->alchFepWCArcut3;
03829   const BigReal alchLambda = simParams->alchLambda;
03830         
03831   const BigReal alchRepLambda = simParams->alchRepLambda;
03832   const BigReal alchDispLambda = simParams->alchDispLambda;
03833   const BigReal alchElecLambda = simParams->alchElecLambda;
03834 
03835   if(stepInRun){
03836     if(!FepWhamOn){
03837       fepFile << FEPTITLE(step);
03838       fepFile << FORMAT(eeng);
03839       fepFile << FORMAT(eeng_f);
03840       fepFile << FORMAT(ljEnergy);
03841       fepFile << FORMAT(ljEnergy_f);
03842     }
03843     else{ // FepWhamOn = ON
03844       if(WCARepuOn){
03845         if(WCArcut1<WCArcut2) { // [s1,s2]
03846           fepFile << "FEP_WCA_REP  ";
03847           fepFile << FORMAT(WCArcut1);
03848           fepFile << FORMAT(WCArcut2);
03849           fepFile << FORMAT(1.0);
03850           fepFile << FORMAT(dE_Left);
03851         }
03852         if(WCArcut2<WCArcut3) { // [s2,s3]
03853           if(WCArcut1<WCArcut2) fepFile << " BREAK ";
03854           fepFile << "FEP_WCA_REP  ";
03855           fepFile << FORMAT(WCArcut2);
03856           fepFile << FORMAT(WCArcut3);
03857           fepFile << FORMAT(0.0);
03858           fepFile << FORMAT(dE);
03859         }
03860         fepFile << std::endl;
03861       }
03862       else if(simParams->alchFepWCADispOn) {
03863         fepFile << "FEP_WCA_DISP ";
03864         fepFile << FORMAT(alchDispLambda);
03865       }
03866       else if(simParams->alchFepElecOn) {
03867         fepFile << "FEP_ELEC     ";
03868         fepFile << FORMAT(alchElecLambda);
03869       }
03870     }
03871     if( ! WCARepuOn ) {
03872       fepFile << FORMAT(dE);
03873     }
03874     if(alchEnsembleAvg){
03875       BigReal dE_avg = net_dE/FepNo;
03876       fepFile << FORMAT(dE_avg);
03877     }
03878     if(!FepWhamOn){
03879       fepFile << FORMAT(temperature);
03880     }
03881     if(alchEnsembleAvg){
03882       dG = -(RT * log(exp_dE_ByRT/FepNo));
03883       fepFile << FORMAT(dG);
03884     } 
03885     if( ! WCARepuOn ) {
03886       fepFile << std::endl;
03887     }
03888   }
03889 }

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

Definition at line 3891 of file Controller.C.

References SimParameters::alchLambdaFreq, cumAlchWork, FORMAT(), net_dEdl_bond_1, net_dEdl_bond_2, net_dEdl_elec_1, net_dEdl_elec_2, net_dEdl_lj_1, net_dEdl_lj_2, 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 TITITLE().

Referenced by outputTiEnergy().


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

BigReal Controller::accelMDdVAverage [protected]

Definition at line 291 of file Controller.h.

Referenced by rescaleaccelMD().

Bool Controller::adaptTempAutoDt [protected]

Definition at line 314 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempBetaMax [protected]

Definition at line 308 of file Controller.h.

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

BigReal Controller::adaptTempBetaMin [protected]

Definition at line 307 of file Controller.h.

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

BigReal* Controller::adaptTempBetaN [protected]

Definition at line 303 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

int Controller::adaptTempBin [protected]

Definition at line 309 of file Controller.h.

Referenced by adaptTempUpdate().

int Controller::adaptTempBins [protected]

Definition at line 310 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempCg [protected]

Definition at line 312 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempDBeta [protected]

Definition at line 311 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal Controller::adaptTempDt [protected]

Definition at line 313 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempDTave [protected]

Definition at line 305 of file Controller.h.

BigReal Controller::adaptTempDTavenum [protected]

Definition at line 306 of file Controller.h.

BigReal Controller::adaptTempDtMax [protected]

Definition at line 316 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempDtMin [protected]

Definition at line 315 of file Controller.h.

Referenced by adaptTempInit().

BigReal* Controller::adaptTempPotEnergyAve [protected]

Definition at line 300 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyAveDen [protected]

Definition at line 298 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyAveNum [protected]

Definition at line 297 of file Controller.h.

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

int* Controller::adaptTempPotEnergySamples [protected]

Definition at line 302 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyVar [protected]

Definition at line 301 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyVarNum [protected]

Definition at line 299 of file Controller.h.

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

ofstream_namd Controller::adaptTempRestartFile [protected]

Definition at line 317 of file Controller.h.

Referenced by adaptTempWriteRestart().

BigReal Controller::adaptTempT [protected]

Definition at line 304 of file Controller.h.

Referenced by adaptTempInit(), and rescaleVelocities().

BigReal Controller::alchWork [protected]

Definition at line 150 of file Controller.h.

Referenced by outputTiEnergy(), and printEnergies().

RequireReduction* Controller::amd_reduction [protected]

Definition at line 228 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().

BigReal Controller::bondedEnergy_ti_1 [protected]

Definition at line 125 of file Controller.h.

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

BigReal Controller::bondedEnergy_ti_2 [protected]

Definition at line 126 of file Controller.h.

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

BigReal Controller::bondedEnergyDiff_f [protected]

Definition at line 112 of file Controller.h.

Referenced by outputFepEnergy(), and printEnergies().

ControllerBroadcasts* Controller::broadcast [protected]

Definition at line 241 of file Controller.h.

Referenced by algorithm(), berendsenPressure(), Controller(), correctMomentum(), cycleBarrier(), langevinPiston1(), minimize(), multigratorPressure(), multigratorTemperature(), rescaleaccelMD(), rescaleVelocities(), resumeAfterTraceBarrier(), stochRescaleVelocities(), tcoupleVelocities(), and ~Controller().

Lattice Controller::checkpoint_lattice [protected]

Definition at line 259 of file Controller.h.

Referenced by algorithm().

ControllerState Controller::checkpoint_state [protected]

Definition at line 260 of file Controller.h.

Referenced by algorithm().

int Controller::checkpoint_stored [protected]

Definition at line 258 of file Controller.h.

Referenced by algorithm(), and Controller().

int Controller::checkpoint_task [protected]

Definition at line 267 of file Controller.h.

Referenced by recvCheckpointAck().

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

Definition at line 266 of file Controller.h.

Referenced by algorithm(), and recvCheckpointReq().

CollectionMaster* const Controller::collection [protected]

Definition at line 239 of file Controller.h.

Referenced by enqueueCollections().

int Controller::computeChecksum [protected]

Definition at line 88 of file Controller.h.

Referenced by compareChecksums().

int Controller::controlNumDegFreedom [protected]

Definition at line 165 of file Controller.h.

Referenced by langevinPiston1(), langevinPiston2(), multigatorCalcEnthalpy(), multigratorPressure(), and receivePressure().

Tensor Controller::controlPressure [protected]

Definition at line 166 of file Controller.h.

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

Tensor Controller::controlPressure_nbond [protected]

Definition at line 76 of file Controller.h.

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

Tensor Controller::controlPressure_normal [protected]

Definition at line 75 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::controlPressure_slow [protected]

Definition at line 77 of file Controller.h.

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

BigReal Controller::cumAlchWork [protected]

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().

BigReal Controller::drudeBondTemp [protected]

Definition at line 154 of file Controller.h.

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

BigReal Controller::drudeBondTempAvg [protected]

Definition at line 155 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::electEnergy [protected]

Definition at line 103 of file Controller.h.

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

BigReal Controller::electEnergy_f [protected]

Definition at line 113 of file Controller.h.

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

BigReal Controller::electEnergy_ti_1 [protected]

Definition at line 127 of file Controller.h.

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

BigReal Controller::electEnergy_ti_2 [protected]

Definition at line 130 of file Controller.h.

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

BigReal Controller::electEnergyPME_ti_1 [protected]

Definition at line 140 of file Controller.h.

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

BigReal Controller::electEnergyPME_ti_2 [protected]

Definition at line 141 of file Controller.h.

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

BigReal Controller::electEnergySlow [protected]

Definition at line 104 of file Controller.h.

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

BigReal Controller::electEnergySlow_f [protected]

Definition at line 114 of file Controller.h.

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

BigReal Controller::electEnergySlow_ti_1 [protected]

Definition at line 128 of file Controller.h.

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

BigReal Controller::electEnergySlow_ti_2 [protected]

Definition at line 131 of file Controller.h.

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

BigReal Controller::exp_dE_ByRT [protected]

Definition at line 117 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

ofstream_namd Controller::fepFile [protected]

Definition at line 248 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().

BigReal Controller::fepSum [protected]

Definition at line 123 of file Controller.h.

Referenced by outputFepEnergy().

int Controller::fflush_count [protected]

Definition at line 212 of file Controller.h.

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

BigReal Controller::goNativeEnergy [protected]

Definition at line 108 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::goNonnativeEnergy [protected]

Definition at line 109 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::goTotalEnergy [protected]

Definition at line 110 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::groGaussEnergy [protected]

Definition at line 107 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::groLJEnergy [protected]

Definition at line 106 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Tensor Controller::groupPressure [protected]

Definition at line 164 of file Controller.h.

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

BigReal Controller::groupPressure_avg [protected]

Definition at line 82 of file Controller.h.

Referenced by Controller(), and printEnergies().

Tensor Controller::groupPressure_nbond [protected]

Definition at line 73 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_normal [protected]

Definition at line 72 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_slow [protected]

Definition at line 74 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_tavg [protected]

Definition at line 85 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::kineticEnergy [protected]

Definition at line 157 of file Controller.h.

Referenced by adaptTempUpdate(), multigatorCalcEnthalpy(), multigratorPressure(), multigratorTemperature(), printEnergies(), and receivePressure().

BigReal Controller::kineticEnergyCentered [protected]

Definition at line 159 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

BigReal Controller::kineticEnergyHalfstep [protected]

Definition at line 158 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Tensor Controller::langevinPiston_origStrainRate [protected]

Definition at line 194 of file Controller.h.

Referenced by Controller().

int Controller::ldbSteps [protected]

Definition at line 210 of file Controller.h.

Referenced by rebalanceLoad().

BigReal Controller::ljEnergy [protected]

Definition at line 105 of file Controller.h.

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

BigReal Controller::ljEnergy_f [protected]

Definition at line 115 of file Controller.h.

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

BigReal Controller::ljEnergy_f_left [protected]

Definition at line 116 of file Controller.h.

Referenced by printEnergies(), and writeFepEnergyData().

BigReal Controller::ljEnergy_ti_1 [protected]

Definition at line 129 of file Controller.h.

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

BigReal Controller::ljEnergy_ti_2 [protected]

Definition at line 132 of file Controller.h.

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

int Controller::marginViolations [protected]

Definition at line 89 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

BigReal Controller::min_energy [protected]

Definition at line 93 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

BigReal Controller::min_f_dot_f [protected]

Definition at line 94 of file Controller.h.

Referenced by minimize().

BigReal Controller::min_f_dot_v [protected]

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().

RequireReduction* Controller::min_reduction [protected]

Definition at line 59 of file Controller.h.

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

BigReal Controller::min_v_dot_v [protected]

Definition at line 96 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Tensor Controller::momentumSqrSum [protected]

Definition at line 201 of file Controller.h.

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

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

Definition at line 203 of file Controller.h.

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

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

Definition at line 204 of file Controller.h.

Referenced by Controller(), and multigratorTemperature().

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

Definition at line 205 of file Controller.h.

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

RequireReduction* Controller::multigratorReduction [protected]

Definition at line 207 of file Controller.h.

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

BigReal Controller::multigratorXi [protected]

Definition at line 199 of file Controller.h.

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

BigReal Controller::multigratorXiT [protected]

Definition at line 200 of file Controller.h.

Referenced by multigratorPressure().

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

Definition at line 206 of file Controller.h.

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

int Controller::nbondFreq [protected]

Definition at line 78 of file Controller.h.

Referenced by calcPressure(), integrate(), langevinPiston1(), langevinPiston2(), minimize(), printEnergies(), and rescaleaccelMD().

BigReal Controller::net_dE [protected]

Definition at line 119 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::net_dEdl_bond_1 [protected]

Definition at line 133 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_bond_2 [protected]

Definition at line 134 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_elec_1 [protected]

Definition at line 135 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_elec_2 [protected]

Definition at line 136 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_lj_1 [protected]

Definition at line 137 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_lj_2 [protected]

Definition at line 138 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int64_t Controller::numDegFreedom [protected]

Definition at line 100 of file Controller.h.

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

Lattice Controller::origLattice [protected]

Definition at line 271 of file Controller.h.

Referenced by Controller().

int Controller::pairlistWarnings [protected]

Definition at line 90 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Tensor Controller::positionRescaleFactor [protected]

Definition at line 196 of file Controller.h.

Referenced by langevinPiston1().

PressureProfileReduction* Controller::ppbonded [protected]

Definition at line 232 of file Controller.h.

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

PressureProfileReduction* Controller::ppint [protected]

Definition at line 234 of file Controller.h.

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

PressureProfileReduction* Controller::ppnonbonded [protected]

Definition at line 233 of file Controller.h.

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

Tensor Controller::pressure [protected]

Definition at line 163 of file Controller.h.

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

Tensor Controller::pressure_amd [protected]

Definition at line 70 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

BigReal Controller::pressure_avg [protected]

Definition at line 81 of file Controller.h.

Referenced by Controller(), and printEnergies().

Tensor Controller::pressure_nbond [protected]

Definition at line 68 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_normal [protected]

Definition at line 67 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_slow [protected]

Definition at line 69 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_tavg [protected]

Definition at line 84 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal* Controller::pressureProfileAverage [protected]

Definition at line 237 of file Controller.h.

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

int Controller::pressureProfileCount [protected]

Definition at line 236 of file Controller.h.

Referenced by Controller(), and printEnergies().

int Controller::pressureProfileSlabs [protected]

Definition at line 235 of file Controller.h.

Referenced by Controller(), and printEnergies().

Random* Controller::random [protected]

Definition at line 224 of file Controller.h.

Referenced by Controller(), langevinPiston1(), langevinPiston2(), stochRescaleVelocities(), and ~Controller().

BigReal Controller::recent_alchWork [protected]

Definition at line 149 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_bond_1 [protected]

Definition at line 143 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_bond_2 [protected]

Definition at line 144 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_elec_1 [protected]

Definition at line 145 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_elec_2 [protected]

Definition at line 146 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_lj_1 [protected]

Definition at line 147 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_lj_2 [protected]

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().

RequireReduction* Controller::reduction [protected]

Definition at line 227 of file Controller.h.

Referenced by compareChecksums(), Controller(), correctMomentum(), multigratorPressure(), printEnergies(), printMinimizeEnergies(), receivePressure(), and ~Controller().

int Controller::rescaleVelocities_numTemps [protected]

Definition at line 171 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

BigReal Controller::rescaleVelocities_sumTemps [protected]

Definition at line 170 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

SimParameters* const Controller::simParams [protected]

Definition at line 225 of file Controller.h.

Referenced by adaptTempInit(), adaptTempUpdate(), adaptTempWriteRestart(), algorithm(), berendsenPressure(), compareChecksums(), computeAlchWork(), Controller(), integrate(), langevinPiston1(), langevinPiston2(), minimize(), multigatorCalcEnthalpy(), multigratorPressure(), multigratorTemperature(), outputExtendedSystem(), outputFepEnergy(), outputTiEnergy(), printEnergies(), printFepMessage(), printTiMessage(), printTiming(), reassignVelocities(), receivePressure(), rescaleaccelMD(), rescaleVelocities(), stochRescaleVelocities(), tcoupleVelocities(), writeExtendedSystemData(), writeExtendedSystemLabels(), writeFepEnergyData(), and writeTiEnergyData().

int Controller::slowFreq [protected]

Definition at line 79 of file Controller.h.

Referenced by calcPressure(), compareChecksums(), correctMomentum(), integrate(), langevinPiston1(), langevinPiston2(), minimize(), printEnergies(), and rescaleaccelMD().

BigReal Controller::smooth2_avg2 [protected]

Definition at line 162 of file Controller.h.

NamdState* const Controller::state [protected]

Definition at line 226 of file Controller.h.

Referenced by algorithm(), berendsenPressure(), calcPressure(), Controller(), enqueueCollections(), langevinPiston1(), langevinPiston2(), multigatorCalcEnthalpy(), multigratorPressure(), printDynamicsEnergies(), printEnergies(), receivePressure(), recvCheckpointAck(), rescaleaccelMD(), writeExtendedSystemData(), and writeExtendedSystemLabels().

int Controller::stepInFullRun [protected]

Definition at line 101 of file Controller.h.

Referenced by printEnergies().

int Controller::stochRescale_count [protected]

Count time steps until next stochastic velocity rescaling.

Definition at line 183 of file Controller.h.

Referenced by Controller(), and stochRescaleVelocities().

BigReal Controller::stochRescaleTimefactor [protected]

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

Definition at line 186 of file Controller.h.

Referenced by Controller(), and stochRescaleVelocities().

Tensor Controller::strainRate_old [protected]

Definition at line 195 of file Controller.h.

Referenced by langevinPiston1().

SubmitReduction* Controller::submit_reduction [protected]

Definition at line 229 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().

BigReal Controller::temp_avg [protected]

Definition at line 80 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::temperature [protected]

Definition at line 160 of file Controller.h.

Referenced by multigratorPressure(), multigratorTemperature(), printEnergies(), receivePressure(), rescaleVelocities(), stochRescaleVelocities(), tcoupleVelocities(), and writeFepEnergyData().

ofstream_namd Controller::tiFile [protected]

Definition at line 252 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().

BigReal Controller::totalEnergy [protected]

Definition at line 102 of file Controller.h.

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

Tensor Controller::virial_amd [protected]

Definition at line 71 of file Controller.h.

Referenced by calcPressure(), and rescaleaccelMD().

ofstream_namd Controller::xstFile [protected]

Definition at line 242 of file Controller.h.

Referenced by outputExtendedSystem().


The documentation for this class was generated from the following files:
Generated on Tue Jun 19 01:17:19 2018 for NAMD by  doxygen 1.4.7