NAMD
Classes | Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | Friends | List of all members
Controller Class Reference

#include <Controller.h>

Inheritance diagram for Controller:
ControllerState

Classes

struct  checkpoint
 

Public Member Functions

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

Public Attributes

BigReal accelMDdV
 

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)
 
double stochRescaleCoefficient ()
 
void berendsenPressure (int)
 
void langevinPiston1 (int)
 
void langevinPiston2 (int)
 
void multigratorPressure (int step, int callNumber)
 
void multigratorTemperature (int step, int callNumber)
 
BigReal multigatorCalcEnthalpy (BigReal potentialEnergy, int step, int minimize)
 
void rebalanceLoad (int)
 
void cycleBarrier (int, int)
 
void traceBarrier (int, int)
 
void terminate (void)
 
void outputExtendedSystem (int step)
 
void writeExtendedSystemLabels (ofstream_namd &file)
 
void writeExtendedSystemData (int step, ofstream_namd &file)
 
void outputFepEnergy (int step)
 
void writeFepEnergyData (int step, ofstream_namd &file)
 
void outputTiEnergy (int step)
 
BigReal computeAlchWork (const int step)
 
void writeTiEnergyData (int step, ofstream_namd &file)
 
void recvCheckpointReq (const char *key, int task, checkpoint &cp)
 
void recvCheckpointAck (checkpoint &cp)
 
void calc_accelMDG_mean_std (BigReal testV, int step_n, BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV)
 
void calc_accelMDG_E_k (int iE, int V_n, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal *k0, BigReal *k, BigReal *E, int *iEused, char *warn)
 
void calc_accelMDG_force_factor (BigReal k, BigReal E, BigReal testV, Tensor vir_orig, BigReal *dV, BigReal *factor, Tensor *vir)
 
void write_accelMDG_rest_file (int step_n, char type, int V_n, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal M2, BigReal E, BigReal k, bool write_topic, bool lasttime)
 
void rescaleaccelMD (int step, int minimize=0)
 
void adaptTempInit (int step)
 
void adaptTempUpdate (int step, int minimize=0)
 
void adaptTempWriteRestart (int step)
 

Protected Attributes

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

Friends

class ScriptTcl
 
class Node
 
class CheckpointMsg
 

Detailed Description

Definition at line 39 of file Controller.h.

Constructor & Destructor Documentation

Controller::Controller ( NamdState s)

Definition at line 128 of file Controller.C.

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

128  :
130  simParams(Node::Object()->simParameters),
131  state(s),
133  startCTime(0),
134  firstCTime(CmiTimer()),
135  startWTime(0),
136  firstWTime(CmiWallTimer()),
137  startBenchTime(0),
138  stepInFullRun(0),
139  ldbSteps(0),
140  fflush_count(3)
141 {
144  // for accelMD
145  if (simParams->accelMDOn) {
146  // REDUCTIONS_BASIC wil contain all potential energies and arrive first
148  // contents of amd_reduction be submitted to REDUCTIONS_AMD
150  // REDUCTIONS_AMD will contain Sequencer contributions and arrive second
152  } else {
153  amd_reduction = NULL;
154  submit_reduction = NULL;
156  }
157  // pressure profile reductions
160  ppbonded = ppnonbonded = ppint = NULL;
162  int ntypes = simParams->pressureProfileAtomTypes;
164  // number of partitions for pairwise interactions
165  int npairs = (ntypes * (ntypes+1))/2;
166  pressureProfileAverage = new BigReal[3*nslabs];
167  memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
168  int freq = simParams->pressureProfileFreq;
171  nslabs, npairs, "BONDED", freq);
173  nslabs, npairs, "NONBONDED", freq);
174  // internal partitions by atom type, but not in a pairwise fashion
176  nslabs, ntypes, "INTERNAL", freq);
177  } else {
178  // just doing Ewald, so only need nonbonded
180  nslabs, npairs, "NONBONDED", freq);
181  }
182  }
185 
186  heat = totalEnergy0 = 0.0;
189  stochRescale_count = 0;
190  if (simParams->stochRescaleOn) {
193  }
195  // strainRate tensor is symmetric to avoid rotation
198  if ( ! simParams->useFlexibleCell ) {
199  BigReal avg = trace(langevinPiston_strainRate) / 3.;
201  } else if ( simParams->useConstantRatio ) {
202 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
203  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
205 #undef AVGXY
206  }
208  if (simParams->multigratorOn) {
209  multigratorXi = 0.0;
212  Node *node = Node::Object();
213  Molecule *molecule = node->molecule;
214  BigReal Nf = molecule->num_deg_freedom();
216  multigratorNu.resize(n);
217  multigratorNuT.resize(n);
218  multigratorZeta.resize(n);
219  multigratorOmega.resize(n);
220  for (int i=0;i < n;i++) {
221  multigratorNu[i] = 0.0;
222  multigratorZeta[i] = 0.0;
223  if (i == 0) {
224  multigratorOmega[i] = Nf*kT0*tau*tau;
225  } else {
226  multigratorOmega[i] = kT0*tau*tau;
227  }
228  }
230  } else {
231  multigratorReduction = NULL;
232  }
233  origLattice = state->lattice;
235  temp_avg = 0;
236  pressure_avg = 0;
237  groupPressure_avg = 0;
238  avg_count = 0;
239  pressure_tavg = 0;
240  groupPressure_tavg = 0;
241  tavg_count = 0;
242  checkpoint_stored = 0;
243  drudeBondTemp = 0;
244  drudeBondTempAvg = 0;
245  cumAlchWork = 0;
246 }
static Node * Object()
Definition: Node.h:86
int checkpoint_stored
Definition: Controller.h:268
#define AVGXY(T)
#define XXXBIGREAL
Definition: Controller.C:60
int pressureProfileSlabs
Definition: Controller.h:245
BigReal smooth2_avg
Definition: Controller.h:36
#define BOLTZMANN
Definition: common.h:45
Definition: Node.h:78
BigReal temp_avg
Definition: Controller.h:81
int fflush_count
Definition: Controller.h:222
static PatchMap * Object()
Definition: PatchMap.h:27
BigReal totalEnergy0
Definition: Controller.h:164
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t cudaTextureObject_t const int const float const PatchPairRecord *__restrict__ const int *__restrict__ const int2 *__restrict__ const unsigned int *__restrict__ unsigned int *__restrict__ int *__restrict__ int *__restrict__ TileListStat *__restrict__ const BoundingBox *__restrict__ float4 *__restrict__ float4 *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ const int numPatches
BigReal pressure_avg
Definition: Controller.h:82
std::vector< BigReal > multigratorOmega
Definition: Controller.h:215
std::vector< BigReal > multigratorNuT
Definition: Controller.h:214
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
Tensor groupPressure_tavg
Definition: Controller.h:86
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
int pressureProfileSlabs
int stochRescale_count
Definition: Controller.h:192
void split(int iStream, int numStreams)
Definition: Random.h:77
RequireReduction * amd_reduction
Definition: Controller.h:238
int64_t num_deg_freedom(int isInitialReport=0) const
Definition: Molecule.h:523
BigReal heat
Definition: Controller.h:162
Tensor langevinPiston_strainRate
Definition: Controller.h:33
int pressureProfileCount
Definition: Controller.h:246
std::vector< BigReal > multigratorNu
Definition: Controller.h:213
NamdState *const state
Definition: Controller.h:236
PressureProfileReduction * ppint
Definition: Controller.h:244
BigReal stochRescalePeriod
Definition: Random.h:37
ControllerBroadcasts * broadcast
Definition: Controller.h:251
BigReal groupPressure_avg
Definition: Controller.h:83
BigReal drudeBondTempAvg
Definition: Controller.h:156
std::vector< BigReal > multigratorZeta
Definition: Controller.h:216
BigReal drudeBondTemp
Definition: Controller.h:155
zVector strainRate2
int computeChecksum
Definition: Controller.h:89
int marginViolations
Definition: Controller.h:90
SubmitReduction * submit_reduction
Definition: Controller.h:239
BigReal rescaleVelocities_sumTemps
Definition: Controller.h:174
int berendsenPressure_count
Definition: Controller.h:35
RequireReduction * multigratorReduction
Definition: Controller.h:217
int multigratorNoseHooverChainLength
BigReal multigratorTemperatureTarget
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
int ldbSteps
Definition: Controller.h:220
zVector strainRate
Bool pressureProfileEwaldOn
Random * random
Definition: Controller.h:234
BigReal multigratorTemperatureRelaxationTime
static CollectionMaster * Object()
int pairlistWarnings
Definition: Controller.h:91
unsigned int randomSeed
CollectionMaster *const collection
Definition: Controller.h:249
int pressureProfileAtomTypes
BigReal * pressureProfileAverage
Definition: Controller.h:247
RequireReduction * min_reduction
Definition: Controller.h:60
BigReal stochRescaleTimefactor
Definition: Controller.h:195
BigReal cumAlchWork
Definition: Controller.h:140
Tensor pressure_tavg
Definition: Controller.h:85
int rescaleVelocities_numTemps
Definition: Controller.h:175
Bool pressureProfileOn
PressureProfileReduction * ppnonbonded
Definition: Controller.h:243
RequireReduction * willRequire(int setID, int size=-1)
Definition: ReductionMgr.C:526
static Tensor symmetric(const Vector &v1, const Vector &v2)
Definition: Tensor.h:45
SimParameters *const simParams
Definition: Controller.h:235
int tavg_count
Definition: Controller.h:87
RequireReduction * reduction
Definition: Controller.h:237
PressureProfileReduction * ppbonded
Definition: Controller.h:242
Lattice origLattice
Definition: Controller.h:281
Tensor langevinPiston_origStrainRate
Definition: Controller.h:204
Molecule * molecule
Definition: Node.h:176
Tensor berendsenPressure_avg
Definition: Controller.h:34
int stepInFullRun
Definition: Controller.h:102
int avg_count
Definition: Controller.h:84
BigReal multigratorXi
Definition: Controller.h:209
double BigReal
Definition: common.h:112
Controller::~Controller ( void  )
virtual

Definition at line 248 of file Controller.C.

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

249 {
250  delete broadcast;
251  delete reduction;
252  delete min_reduction;
253  delete amd_reduction;
254  delete submit_reduction;
255  delete ppbonded;
256  delete ppnonbonded;
257  delete ppint;
258  delete [] pressureProfileAverage;
259  delete random;
261 }
RequireReduction * amd_reduction
Definition: Controller.h:238
PressureProfileReduction * ppint
Definition: Controller.h:244
ControllerBroadcasts * broadcast
Definition: Controller.h:251
SubmitReduction * submit_reduction
Definition: Controller.h:239
RequireReduction * multigratorReduction
Definition: Controller.h:217
Random * random
Definition: Controller.h:234
BigReal * pressureProfileAverage
Definition: Controller.h:247
RequireReduction * min_reduction
Definition: Controller.h:60
PressureProfileReduction * ppnonbonded
Definition: Controller.h:243
RequireReduction * reduction
Definition: Controller.h:237
PressureProfileReduction * ppbonded
Definition: Controller.h:242

Member Function Documentation

void Controller::adaptTempInit ( int  step)
protected

Definition at line 2343 of file Controller.C.

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

Referenced by adaptTempUpdate().

2343  {
2344  if (!simParams->adaptTempOn) return;
2345  iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
2346  adaptTempDtMin = 0;
2347  adaptTempDtMax = 0;
2348  adaptTempAutoDt = false;
2349  if (simParams->adaptTempBins == 0) {
2350  iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
2351  std::ifstream adaptTempRead(simParams->adaptTempInFile);
2352  if (adaptTempRead) {
2353  int readInt;
2354  BigReal readReal;
2355  bool readBool;
2356  // step
2357  adaptTempRead >> readInt;
2358  // Start with min and max temperatures
2359  adaptTempRead >> adaptTempT; // KELVIN
2360  adaptTempRead >> adaptTempBetaMin; // KELVIN
2361  adaptTempRead >> adaptTempBetaMax; // KELVIN
2362  adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
2363  adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
2364  // In case file is manually edited
2365  if (adaptTempBetaMin > adaptTempBetaMax){
2366  readReal = adaptTempBetaMax;
2367  adaptTempBetaMax = adaptTempBetaMin;
2368  adaptTempBetaMin = adaptTempBetaMax;
2369  }
2370  adaptTempRead >> adaptTempBins;
2371  adaptTempRead >> adaptTempCg;
2372  adaptTempRead >> adaptTempDt;
2380  adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
2381  for(int j = 0; j < adaptTempBins; ++j) {
2382  adaptTempRead >> adaptTempPotEnergyAve[j];
2383  adaptTempRead >> adaptTempPotEnergyVar[j];
2384  adaptTempRead >> adaptTempPotEnergySamples[j];
2385  adaptTempRead >> adaptTempPotEnergyAveNum[j];
2386  adaptTempRead >> adaptTempPotEnergyVarNum[j];
2387  adaptTempRead >> adaptTempPotEnergyAveDen[j];
2388  adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
2389  }
2390  adaptTempRead.close();
2391  }
2392  else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
2393  }
2394  else {
2395  adaptTempBins = simParams->adaptTempBins;
2403  adaptTempBetaMax = 1./simParams->adaptTempTmin;
2404  adaptTempBetaMin = 1./simParams->adaptTempTmax;
2405  adaptTempCg = simParams->adaptTempCgamma;
2406  adaptTempDt = simParams->adaptTempDt;
2407  adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
2408  adaptTempT = simParams->initialTemp;
2409  if (simParams->langevinOn)
2410  adaptTempT = simParams->langevinTemp;
2411  else if (simParams->rescaleFreq > 0)
2412  adaptTempT = simParams->rescaleTemp;
2413  for(int j = 0; j < adaptTempBins; ++j){
2414  adaptTempPotEnergyAveNum[j] = 0.;
2415  adaptTempPotEnergyAveDen[j] = 0.;
2417  adaptTempPotEnergyVarNum[j] = 0.;
2418  adaptTempPotEnergyVar[j] = 0.;
2419  adaptTempPotEnergyAve[j] = 0.;
2420  adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
2421  }
2422  }
2423  if (simParams->adaptTempAutoDt > 0.0) {
2424  adaptTempAutoDt = true;
2427  }
2428  adaptTempDTave = 0;
2429  adaptTempDTavenum = 0;
2430  iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
2431  iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
2432  iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
2433  iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
2434  iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
2435  if (simParams->adaptTempRestartFreq > 0) {
2439  NAMD_die("Error opening restart file");
2440  }
2441 }
BigReal adaptTempBetaMin
Definition: Controller.h:317
BigReal adaptTempTmax
int adaptTempRestartFreq
BigReal * adaptTempBetaN
Definition: Controller.h:313
int adaptTempBins
Definition: Controller.h:320
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
BigReal adaptTempCgamma
BigReal adaptTempTmin
BigReal adaptTempT
Definition: Controller.h:314
BigReal adaptTempDt
Definition: Controller.h:323
BigReal * adaptTempPotEnergyAveNum
Definition: Controller.h:307
#define iout
Definition: InfoStream.h:87
BigReal adaptTempDBeta
Definition: Controller.h:321
BigReal adaptTempDTavenum
Definition: Controller.h:316
BigReal adaptTempDTave
Definition: Controller.h:315
ofstream_namd adaptTempRestartFile
Definition: Controller.h:327
BigReal adaptTempDtMin
Definition: Controller.h:325
BigReal * adaptTempPotEnergyVar
Definition: Controller.h:311
BigReal rescaleTemp
char adaptTempRestartFile[128]
BigReal adaptTempDtMax
Definition: Controller.h:326
BigReal langevinTemp
char adaptTempInFile[128]
Bool adaptTempAutoDt
Definition: Controller.h:324
void NAMD_die(const char *err_msg)
Definition: common.C:83
BigReal adaptTempDt
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:159
BigReal initialTemp
BigReal * adaptTempPotEnergyAveDen
Definition: Controller.h:308
BigReal * adaptTempPotEnergyVarNum
Definition: Controller.h:309
BigReal adaptTempCg
Definition: Controller.h:322
int * adaptTempPotEnergySamples
Definition: Controller.h:312
BigReal adaptTempBetaMax
Definition: Controller.h:318
SimParameters *const simParams
Definition: Controller.h:235
BigReal adaptTempAutoDt
infostream & endi(infostream &s)
Definition: InfoStream.C:38
BigReal * adaptTempPotEnergyAve
Definition: Controller.h:310
void open(const char *_fname, std::ios_base::openmode _mode=std::ios_base::out)
Definition: fstream_namd.C:11
double BigReal
Definition: common.h:112
void Controller::adaptTempUpdate ( int  step,
int  minimize = 0 
)
protected

Definition at line 2469 of file Controller.C.

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

Referenced by integrate().

2470 {
2471  //Beta = 1./T
2472  if ( !simParams->adaptTempOn ) return;
2473  int j = 0;
2474  if (step == simParams->firstTimestep) {
2475  adaptTempInit(step);
2476  return;
2477  }
2478  if ( minimize || (step < simParams->adaptTempFirstStep ) ||
2479  ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
2480  const int adaptTempOutFreq = simParams->adaptTempOutFreq;
2481  const bool adaptTempDebug = simParams->adaptTempDebug;
2482  //Calculate Current inverse temperature and bin
2483  BigReal adaptTempBeta = 1./adaptTempT;
2484  adaptTempBin = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
2485 
2486  if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
2487  iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin
2488  << " adaptTempBeta: " << adaptTempBeta
2489  << " adaptTempDBeta: " << adaptTempDBeta
2490  << " betaMin:" << adaptTempBetaMin
2491  << " betaMax: " << adaptTempBetaMax << "\n";
2494 
2495  BigReal potentialEnergy;
2496  BigReal potEnergyAverage;
2497  BigReal potEnergyVariance;
2498  potentialEnergy = totalEnergy-kineticEnergy;
2499 
2500  //calculate new bin average and variance using adaptive averaging
2503  adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
2504 
2507  potEnergyVariance -= potEnergyAverage*potEnergyAverage;
2508 
2509  adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
2510  adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
2511 
2512  // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
2513  // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
2514  // is constant for each bin. This is to estimate <E(beta)> where beta \in
2515  // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
2516  if ( ! ( step % simParams->adaptTempFreq ) ) {
2517  // If adaptTempBin not at the edge, calculate improved average:
2518  if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
2519  // Get Averaging Limits:
2520  BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
2521  BigReal betaPlus;
2522  BigReal betaMinus;
2523  int nMinus =0;
2524  int nPlus = 0;
2525  if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
2526  if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
2527  betaMinus = adaptTempBeta - deltaBeta;
2528  betaPlus = adaptTempBeta + deltaBeta;
2529  BigReal betaMinus2 = betaMinus*betaMinus;
2530  BigReal betaPlus2 = betaPlus*betaPlus;
2531  nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
2532  nPlus = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
2533  // Variables for <E(beta)> estimate:
2534  BigReal potEnergyAve0 = 0.0;
2535  BigReal potEnergyAve1 = 0.0;
2536  // Integral terms
2537  BigReal A0 = 0;
2538  BigReal A1 = 0;
2539  BigReal A2 = 0;
2540  //A0 phi_s integral for beta_minus < beta < beta_{i+1}
2541  BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1];
2542  BigReal DeltaE2Ave;
2543  BigReal DeltaE2AveSum = 0;
2544  BigReal betaj;
2545  for (j = nMinus+1; j <= adaptTempBin; ++j) {
2546  potEnergyAve0 += adaptTempPotEnergyAve[j];
2547  DeltaE2Ave = adaptTempPotEnergyVar[j];
2548  DeltaE2AveSum += DeltaE2Ave;
2549  A0 += j*DeltaE2Ave;
2550  }
2551  A0 *= adaptTempDBeta;
2552  A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
2553  A0 *= adaptTempDBeta;
2554  betaj = adaptTempBetaN[nMinus+1]-betaMinus;
2555  betaj *= betaj;
2556  A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
2557  A0 /= (betaNp1-betaMinus);
2558 
2559  //A1 phi_s integral for beta_{i+1} < beta < beta_plus
2560  DeltaE2AveSum = 0;
2561  for (j = adaptTempBin+1; j < nPlus; j++) {
2562  potEnergyAve1 += adaptTempPotEnergyAve[j];
2563  DeltaE2Ave = adaptTempPotEnergyVar[j];
2564  DeltaE2AveSum += DeltaE2Ave;
2565  A1 += j*DeltaE2Ave;
2566  }
2567  A1 *= adaptTempDBeta;
2568  A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
2569  A1 *= adaptTempDBeta;
2570  if ( nPlus < adaptTempBins ) {
2571  betaj = betaPlus-adaptTempBetaN[nPlus];
2572  betaj *= betaj;
2573  A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
2574  }
2575  A1 /= (betaPlus-betaNp1);
2576 
2577  //A2 phi_t integral for beta_i
2578  A2 = 0.5*adaptTempDBeta*potEnergyVariance;
2579 
2580  // Now calculate a+ and a-
2581  DeltaE2Ave = A0-A1;
2582  BigReal aplus = 0;
2583  BigReal aminus = 0;
2584  if (DeltaE2Ave != 0) {
2585  aplus = (A0-A2)/(A0-A1);
2586  if (aplus < 0) {
2587  aplus = 0;
2588  }
2589  if (aplus > 1) {
2590  aplus = 1;
2591  }
2592  aminus = 1-aplus;
2593  potEnergyAve0 *= adaptTempDBeta;
2594  potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
2595  potEnergyAve0 /= (betaNp1-betaMinus);
2596  potEnergyAve1 *= adaptTempDBeta;
2597  if ( nPlus < adaptTempBins ) {
2598  potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
2599  }
2600  potEnergyAve1 /= (betaPlus-betaNp1);
2601  potEnergyAverage = aminus*potEnergyAve0;
2602  potEnergyAverage += aplus*potEnergyAve1;
2603  }
2604  if (simParams->adaptTempDebug) {
2605  iout << "ADAPTEMP DEBUG:" << "\n"
2606  << " adaptTempBin: " << adaptTempBin << "\n"
2607  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
2608  << " adaptTempBeta: " << adaptTempBeta << "\n"
2609  << " adaptTemp: " << adaptTempT<< "\n"
2610  << " betaMin: " << adaptTempBetaMin << "\n"
2611  << " betaMax: " << adaptTempBetaMax << "\n"
2612  << " gammaAve: " << gammaAve << "\n"
2613  << " deltaBeta: " << deltaBeta << "\n"
2614  << " betaMinus: " << betaMinus << "\n"
2615  << " betaPlus: " << betaPlus << "\n"
2616  << " nMinus: " << nMinus << "\n"
2617  << " nPlus: " << nPlus << "\n"
2618  << " A0: " << A0 << "\n"
2619  << " A1: " << A1 << "\n"
2620  << " A2: " << A2 << "\n"
2621  << " a+: " << aplus << "\n"
2622  << " a-: " << aminus << "\n"
2623  << endi;
2624  }
2625  }
2626  else {
2627  if (simParams->adaptTempDebug) {
2628  iout << "ADAPTEMP DEBUG:" << "\n"
2629  << " adaptTempBin: " << adaptTempBin << "\n"
2630  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
2631  << " adaptTempBeta: " << adaptTempBeta << "\n"
2632  << " adaptTemp: " << adaptTempT<< "\n"
2633  << " betaMin: " << adaptTempBetaMin << "\n"
2634  << " betaMax: " << adaptTempBetaMax << "\n"
2635  << " gammaAve: " << gammaAve << "\n"
2636  << " deltaBeta: N/A\n"
2637  << " betaMinus: N/A\n"
2638  << " betaPlus: N/A\n"
2639  << " nMinus: N/A\n"
2640  << " nPlus: N/A\n"
2641  << " A0: N/A\n"
2642  << " A1: N/A\n"
2643  << " A2: N/A\n"
2644  << " a+: N/A\n"
2645  << " a-: N/A\n"
2646  << endi;
2647  }
2648  }
2649 
2650  //dT is new temperature
2651  BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
2652  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
2653  dT += adaptTempT;
2654 
2655  // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
2656  // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
2657  if ( dT > 1./adaptTempBetaMin || dT < 1./adaptTempBetaMax ) {
2659  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
2660  dT += adaptTempT;
2661  // Check again, if not then keep original adaptTempTor assign random.
2662  if ( dT > 1./adaptTempBetaMin ) {
2663  if (!simParams->adaptTempRandom) {
2664  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
2665  // << " K higher than adaptTempTmax."
2666  // << " Keeping temperature at "
2667  // << adaptTempT<< "\n"<< endi;
2668  dT = adaptTempT;
2669  }
2670  else {
2671  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
2672  // << " K higher than adaptTempTmax."
2673  // << " Assigning random temperature in range\n" << endi;
2674  dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
2675  dT = 1./dT;
2676  }
2677  }
2678  else if ( dT < 1./adaptTempBetaMax ) {
2679  if (!simParams->adaptTempRandom) {
2680  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
2681  // << " K lower than adaptTempTmin."
2682  // << " Keeping temperature at " << adaptTempT<< "\n" << endi;
2683  dT = adaptTempT;
2684  }
2685  else {
2686  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
2687  // << " K lower than adaptTempTmin."
2688  // << " Assigning random temperature in range\n" << endi;
2689  dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
2690  dT = 1./dT;
2691  }
2692  }
2693  else if (adaptTempAutoDt) {
2694  //update temperature step size counter
2695  //FOR "TRUE" ADAPTIVE TEMPERING
2696  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
2697  if (adaptTempTdiff > 0) {
2698  adaptTempDTave += adaptTempTdiff;
2700 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
2701  }
2702  if(adaptTempDTavenum == 100){
2703  BigReal Frac;
2705  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
2706  Frac = adaptTempDTave/Frac;
2707  //if average temperature jump is > 3% of temperature range,
2708  //modify jump size to match 3%
2709  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
2710  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
2711  Frac = adaptTempDtMax/Frac;
2712  iout << "ADAPTEMP: Updating adaptTempDt to ";
2713  adaptTempDt *= Frac;
2714  iout << adaptTempDt << "\n" << endi;
2715  }
2716  adaptTempDTave = 0;
2717  adaptTempDTavenum = 0;
2718  }
2719  }
2720  }
2721  else if (adaptTempAutoDt) {
2722  //update temperature step size counter
2723  // FOR "TRUE" ADAPTIVE TEMPERING
2724  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
2725  if (adaptTempTdiff > 0) {
2726  adaptTempDTave += adaptTempTdiff;
2728 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
2729  }
2730  if(adaptTempDTavenum == 100){
2731  BigReal Frac;
2733  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
2734  Frac = adaptTempDTave/Frac;
2735  //if average temperature jump is > 3% of temperature range,
2736  //modify jump size to match 3%
2737  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
2738  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
2739  Frac = adaptTempDtMax/Frac;
2740  iout << "ADAPTEMP: Updating adaptTempDt to ";
2741  adaptTempDt *= Frac;
2742  iout << adaptTempDt << "\n" << endi;
2743  }
2744  adaptTempDTave = 0;
2745  adaptTempDTavenum = 0;
2746 
2747  }
2748 
2749  }
2750 
2751  adaptTempT = dT;
2753  }
2754  adaptTempWriteRestart(step);
2755  if ( ! (step % adaptTempOutFreq) ) {
2756  iout << "ADAPTEMP: STEP " << step
2757  << " TEMP " << adaptTempT
2758  << " ENERGY " << std::setprecision(10) << potentialEnergy
2759  << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
2760  << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
2761  iout << "\n" << endi;
2762  }
2763 
2764 }
BigReal adaptTempBetaMin
Definition: Controller.h:317
BigReal * adaptTempBetaN
Definition: Controller.h:313
int adaptTempBins
Definition: Controller.h:320
BigReal uniform(void)
Definition: Random.h:109
#define BOLTZMANN
Definition: common.h:45
BigReal gaussian(void)
Definition: Random.h:116
BigReal adaptTempT
Definition: Controller.h:314
void adaptTempWriteRestart(int step)
Definition: Controller.C:2443
BigReal adaptTempDt
Definition: Controller.h:323
void adaptTempInit(int step)
Definition: Controller.C:2343
void minimize()
Definition: Controller.C:594
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
SimpleBroadcastObject< BigReal > adaptTemperature
Definition: Broadcasts.h:86
BigReal * adaptTempPotEnergyAveNum
Definition: Controller.h:307
#define iout
Definition: InfoStream.h:87
BigReal adaptTempDBeta
Definition: Controller.h:321
BigReal adaptTempDTavenum
Definition: Controller.h:316
BigReal adaptTempDTave
Definition: Controller.h:315
int adaptTempBin
Definition: Controller.h:319
BigReal adaptTempDtMin
Definition: Controller.h:325
ControllerBroadcasts * broadcast
Definition: Controller.h:251
BigReal * adaptTempPotEnergyVar
Definition: Controller.h:311
BigReal adaptTempDtMax
Definition: Controller.h:326
Bool adaptTempAutoDt
Definition: Controller.h:324
void publish(int tag, const T &t)
Random * random
Definition: Controller.h:234
BigReal * adaptTempPotEnergyAveDen
Definition: Controller.h:308
BigReal * adaptTempPotEnergyVarNum
Definition: Controller.h:309
BigReal adaptTempCg
Definition: Controller.h:322
int * adaptTempPotEnergySamples
Definition: Controller.h:312
BigReal totalEnergy
Definition: Controller.h:103
BigReal adaptTempBetaMax
Definition: Controller.h:318
SimParameters *const simParams
Definition: Controller.h:235
infostream & endi(infostream &s)
Definition: InfoStream.C:38
BigReal * adaptTempPotEnergyAve
Definition: Controller.h:310
BigReal kineticEnergy
Definition: Controller.h:158
double BigReal
Definition: common.h:112
void Controller::adaptTempWriteRestart ( int  step)
protected

Definition at line 2443 of file Controller.C.

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

Referenced by adaptTempUpdate().

2443  {
2445  adaptTempRestartFile.seekp(std::ios::beg);
2446  iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
2447  adaptTempRestartFile << step << " ";
2448  // Start with min and max temperatures
2449  adaptTempRestartFile << adaptTempT<< " "; // KELVIN
2450  adaptTempRestartFile << 1./adaptTempBetaMin << " "; // KELVIN
2451  adaptTempRestartFile << 1./adaptTempBetaMax << " "; // KELVIN
2455  adaptTempRestartFile << "\n" ;
2456  for(int j = 0; j < adaptTempBins; ++j) {
2463  adaptTempRestartFile << "\n";
2464  }
2466  }
2467 }
BigReal adaptTempBetaMin
Definition: Controller.h:317
int adaptTempRestartFreq
int adaptTempBins
Definition: Controller.h:320
BigReal adaptTempT
Definition: Controller.h:314
BigReal adaptTempDt
Definition: Controller.h:323
ofstream_namd & flush()
Definition: fstream_namd.C:17
BigReal * adaptTempPotEnergyAveNum
Definition: Controller.h:307
#define iout
Definition: InfoStream.h:87
ofstream_namd adaptTempRestartFile
Definition: Controller.h:327
BigReal * adaptTempPotEnergyVar
Definition: Controller.h:311
BigReal * adaptTempPotEnergyAveDen
Definition: Controller.h:308
BigReal * adaptTempPotEnergyVarNum
Definition: Controller.h:309
BigReal adaptTempCg
Definition: Controller.h:322
int * adaptTempPotEnergySamples
Definition: Controller.h:312
BigReal adaptTempBetaMax
Definition: Controller.h:318
SimParameters *const simParams
Definition: Controller.h:235
infostream & endi(infostream &s)
Definition: InfoStream.C:38
BigReal * adaptTempPotEnergyAve
Definition: Controller.h:310
void Controller::algorithm ( void  )
protectedvirtual

Definition at line 281 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(), if(), SimParameters::initialTemp, integrate(), iout, Controller::checkpoint::lattice, minimize(), NAMD_bug(), NAMD_die(), Node::Object(), outputExtendedSystem(), SCRIPT_ATOMRECV, SCRIPT_ATOMSEND, SCRIPT_ATOMSENDRECV, SCRIPT_CHECKPOINT, SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_STORE, SCRIPT_CHECKPOINT_SWAP, SCRIPT_CONTINUE, SCRIPT_END, SCRIPT_FORCEOUTPUT, SCRIPT_MEASURE, SCRIPT_MINIMIZE, SCRIPT_OUTPUT, SCRIPT_REINITVELS, SCRIPT_RESCALESOLUTECHARGES, SCRIPT_RESCALEVELS, SCRIPT_REVERT, SCRIPT_RUN, SimParameters::scriptArg1, ControllerBroadcasts::scriptBarrier, SimParameters::scriptIntArg1, SimParameters::scriptStringArg1, Node::sendCheckpointReq(), simParams, state, Controller::checkpoint::state, msm::swap(), and terminate().

282 {
283  int scriptTask;
284  int scriptSeq = 0;
285  BackEnd::awaken();
286  while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
287  switch ( scriptTask ) {
288  case SCRIPT_OUTPUT:
291  break;
292  case SCRIPT_FORCEOUTPUT:
294  break;
295  case SCRIPT_MEASURE:
297  break;
298  case SCRIPT_REINITVELS:
299  iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
300  << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
301  break;
302  case SCRIPT_RESCALEVELS:
303  iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
304  << " BY " << simParams->scriptArg1 << "\n" << endi;
305  break;
307  // Parameter setting already reported in ScriptTcl
308  // Nothing to do!
309  break;
310  case SCRIPT_CHECKPOINT:
311  iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
312  << "\n" << endi;
313  checkpoint_stored = 1;
314  checkpoint_lattice = state->lattice;
316  break;
317  case SCRIPT_REVERT:
318  iout << "REVERTING AT STEP " << simParams->firstTimestep
319  << "\n" << endi;
320  if ( ! checkpoint_stored )
321  NAMD_die("Unable to revert, checkpoint was never called!");
322  state->lattice = checkpoint_lattice;
324  break;
326  iout << "STORING CHECKPOINT AT STEP " << simParams->firstTimestep
327  << " TO KEY " << simParams->scriptStringArg1;
328  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
329  iout << "\n" << endi;
330  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
331  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
332  checkpoints[simParams->scriptStringArg1] = new checkpoint;
333  }
334  checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
335  cp.lattice = state->lattice;
336  cp.state = *(ControllerState*)this;
337  } else {
339  scriptTask, state->lattice, *(ControllerState*)this);
340  }
341  break;
343  iout << "LOADING CHECKPOINT AT STEP " << simParams->firstTimestep
344  << " FROM KEY " << simParams->scriptStringArg1;
345  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
346  iout << "\n" << endi;
347  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
348  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
349  NAMD_die("Unable to load checkpoint, requested key was never stored.");
350  }
351  checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
352  state->lattice = cp.lattice;
353  *(ControllerState*)this = cp.state;
354  } else {
356  scriptTask, state->lattice, *(ControllerState*)this);
357  }
358  break;
360  iout << "SWAPPING CHECKPOINT AT STEP " << simParams->firstTimestep
361  << " FROM KEY " << simParams->scriptStringArg1;
362  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
363  iout << "\n" << endi;
364  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
365  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
366  NAMD_die("Unable to swap checkpoint, requested key was never stored.");
367  }
368  checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
369  std::swap(state->lattice,cp.lattice);
370  std::swap(*(ControllerState*)this,cp.state);
371  } else {
373  scriptTask, state->lattice, *(ControllerState*)this);
374  }
375  break;
377  iout << "FREEING CHECKPOINT AT STEP " << simParams->firstTimestep
378  << " FROM KEY " << simParams->scriptStringArg1;
379  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
380  iout << "\n" << endi;
381  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
382  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
383  NAMD_die("Unable to free checkpoint, requested key was never stored.");
384  }
387  } else {
389  scriptTask, state->lattice, *(ControllerState*)this);
390  }
391  break;
392  case SCRIPT_ATOMSENDRECV:
393  case SCRIPT_ATOMSEND:
394  case SCRIPT_ATOMRECV:
395  break;
396  case SCRIPT_MINIMIZE:
397  minimize();
398  break;
399  case SCRIPT_RUN:
400  case SCRIPT_CONTINUE:
401  integrate(scriptTask);
402  break;
403  default:
404  NAMD_bug("Unknown task in Controller::algorithm");
405  }
406  BackEnd::awaken();
407  }
410  terminate();
411 }
static Node * Object()
Definition: Node.h:86
int checkpoint_stored
Definition: Controller.h:268
void enqueueCollections(int)
Definition: Controller.C:3655
char scriptStringArg1[128]
static void awaken(void)
Definition: BackEnd.C:316
std::map< std::string, checkpoint * > checkpoints
Definition: Controller.h:276
#define FILE_OUTPUT
Definition: Output.h:25
void minimize()
Definition: Controller.C:594
#define EVAL_MEASURE
Definition: Output.h:27
void integrate(int)
Definition: Controller.C:429
if(ComputeNonbondedUtil::goMethod==2)
#define iout
Definition: InfoStream.h:87
void swap(Array< T > &s, Array< T > &t)
Definition: MsmMap.h:319
NamdState *const state
Definition: Controller.h:236
Lattice checkpoint_lattice
Definition: Controller.h:269
ControllerState checkpoint_state
Definition: Controller.h:270
ControllerBroadcasts * broadcast
Definition: Controller.h:251
void NAMD_bug(const char *err_msg)
Definition: common.C:123
SimpleBroadcastObject< int > scriptBarrier
Definition: Broadcasts.h:83
BigReal scriptArg1
#define END_OF_RUN
Definition: Output.h:26
void NAMD_die(const char *err_msg)
Definition: common.C:83
void terminate(void)
Definition: Controller.C:4169
BigReal initialTemp
SimParameters *const simParams
Definition: Controller.h:235
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void sendCheckpointReq(int remote, const char *key, int task, Lattice &lat, ControllerState &cs)
Definition: Node.C:1271
void outputExtendedSystem(int step)
Definition: Controller.C:3959
#define FORCE_OUTPUT
Definition: Output.h:28
void Controller::awaken ( void  )
inline

Definition at line 45 of file Controller.h.

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

45 { CthAwaken(thread); };
void Controller::berendsenPressure ( int  step)
protected

Definition at line 950 of file Controller.C.

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

Referenced by integrate().

951 {
955  const int freq = simParams->berendsenPressureFreq;
956  if ( ! (berendsenPressure_count % freq) ) {
960  // We only use on-diagonal terms (for now)
961  factor = Tensor::diagonal(diagonal(factor));
963  factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
965  factor += Tensor::identity(1.0);
966 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
967  if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
968  if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
969  int limited = 0;
970  LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
971  LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
972  LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
973 #undef LIMIT_SCALING
974  if ( limited ) {
975  iout << iERROR << "Step " << step <<
976  " cell rescaling factor limited.\n" << endi;
977  }
979  state->lattice.rescale(factor);
980  }
981  } else {
984  }
985 }
BigReal berendsenPressureCompressibility
Bool berendsenPressureOn
BigReal berendsenPressureRelaxationTime
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
Tensor controlPressure
Definition: Controller.h:170
int berendsenPressureFreq
#define iout
Definition: InfoStream.h:87
virial zz
NamdState *const state
Definition: Controller.h:236
void rescale(Tensor factor)
Definition: Lattice.h:60
BigReal berendsenPressureTarget
ControllerBroadcasts * broadcast
Definition: Controller.h:251
int berendsenPressure_count
Definition: Controller.h:35
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
void publish(int tag, const T &t)
virial yy
BigReal xx
Definition: Tensor.h:17
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
Definition: Tensor.h:15
#define LIMIT_SCALING(VAR, MIN, MAX, FLAG)
SimParameters *const simParams
Definition: Controller.h:235
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Tensor berendsenPressure_avg
Definition: Controller.h:34
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 
)
inlineprotected

Definition at line 1728 of file Controller.C.

Referenced by rescaleaccelMD().

1729  {
1730  switch(iE){
1731  case 2:
1732  *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
1733  // if k0 <=0 OR k0 > 1, switch to iE=1 mode
1734  if(*k0 > 1.)
1735  *warn |= 1;
1736  else if(*k0 <= 0.)
1737  *warn |= 2;
1738  //else stay in iE=2 mode
1739  else{
1740  *E = Vmin + (Vmax-Vmin)/(*k0);
1741  *iEused = 2;
1742  break;
1743  }
1744  case 1:
1745  *E = Vmax;
1746  *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
1747  if(*k0 > 1.)
1748  *k0 = 1.;
1749  *iEused = 1;
1750  break;
1751  }
1752 
1753  *k = *k0 / (Vmax - Vmin);
1754 }
void Controller::calc_accelMDG_force_factor ( BigReal  k,
BigReal  E,
BigReal  testV,
Tensor  vir_orig,
BigReal dV,
BigReal factor,
Tensor vir 
)
inlineprotected

Definition at line 1757 of file Controller.C.

Referenced by rescaleaccelMD().

1758  {
1759  BigReal VE = testV - E;
1760  //if V < E, apply boost
1761  if(VE < 0){
1762  *factor = k * VE;
1763  *vir += vir_orig * (*factor);
1764  *dV += (*factor) * VE/2.;
1765  *factor += 1.;
1766  }
1767  else{
1768  *factor = 1.;
1769  }
1770 }
double BigReal
Definition: common.h:112
void Controller::calc_accelMDG_mean_std ( BigReal  testV,
int  step_n,
BigReal Vmax,
BigReal Vmin,
BigReal Vavg,
BigReal M2,
BigReal sigmaV 
)
inlineprotected

Definition at line 1708 of file Controller.C.

Referenced by rescaleaccelMD().

1709  {
1710  BigReal delta;
1711 
1712  if(testV > *Vmax){
1713  *Vmax = testV;
1714  }
1715  else if(testV < *Vmin){
1716  *Vmin = testV;
1717  }
1718 
1719  //mean and std calculated by Online Algorithm
1720  delta = testV - *Vavg;
1721  *Vavg += delta / (BigReal)n;
1722  *M2 += delta * (testV - (*Vavg));
1723 
1724  *sigmaV = sqrt(*M2/n);
1725 }
double BigReal
Definition: common.h:112
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 1590 of file Controller.C.

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

Referenced by multigratorPressure(), and receivePressure().

1593  {
1594 
1595  Tensor virial_normal = virial_normal_in;
1596  Tensor virial_nbond = virial_nbond_in;
1597  Tensor virial_slow = virial_slow_in;
1598 
1599  // Tensor tmp = virial_normal;
1600  // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
1601  // tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
1602 
1603  Node *node = Node::Object();
1604  Molecule *molecule = node->molecule;
1605  SimParameters *simParameters = node->simParameters;
1606  Lattice &lattice = state->lattice;
1607 
1608  BigReal volume;
1609 
1610  Vector extPosition = lattice.origin();
1611  virial_normal -= outer(extForce_normal,extPosition);
1612  virial_nbond -= outer(extForce_nbond,extPosition);
1613  virial_slow -= outer(extForce_slow,extPosition);
1614 
1615  if ( (volume=lattice.volume()) != 0. )
1616  {
1617 
1618  if (simParameters->LJcorrection && volume) {
1619 #ifdef MEM_OPT_VERSION
1620  NAMD_bug("LJcorrection not supported in memory optimized build.");
1621 #else
1622  // Apply tail correction to pressure
1623  BigReal alchLambda = simParameters->getCurrentLambda(step);
1624  virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
1625 #endif
1626  }
1627 
1628  // kinetic energy component included in virials
1629  pressure_normal = virial_normal / volume;
1630  groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
1631 
1632  if (simParameters->accelMDOn) {
1633  pressure_amd = virial_amd / volume;
1636  }
1637 
1638  if ( minimize || ! ( step % nbondFreq ) )
1639  {
1640  pressure_nbond = virial_nbond / volume;
1641  groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
1642  }
1643 
1644  if ( minimize || ! ( step % slowFreq ) )
1645  {
1646  pressure_slow = virial_slow / volume;
1647  groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
1648  }
1649 
1653  }
1654  else
1655  {
1656  pressure = Tensor();
1657  groupPressure = Tensor();
1658  }
1659 
1660  if ( simParameters->useGroupPressure )
1661  {
1666  }
1667  else
1668  {
1673  }
1674 
1675  if ( simParameters->useFlexibleCell ) {
1676  // use symmetric pressure to control rotation
1677  // controlPressure_normal = symmetric(controlPressure_normal);
1678  // controlPressure_nbond = symmetric(controlPressure_nbond);
1679  // controlPressure_slow = symmetric(controlPressure_slow);
1680  // controlPressure = symmetric(controlPressure);
1681  // only use on-diagonal components for now
1686  if ( simParameters->useConstantRatio ) {
1687 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
1688  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
1693 #undef AVGXY
1694  }
1695  } else {
1702  controlPressure =
1703  Tensor::identity(trace(controlPressure)/3.);
1704  }
1705 }
static Node * Object()
Definition: Node.h:86
#define AVGXY(T)
Tensor controlPressure_slow
Definition: Controller.h:78
int nbondFreq
Definition: Controller.h:79
Tensor groupPressure_nbond
Definition: Controller.h:74
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
Definition: Node.h:78
Tensor controlPressure
Definition: Controller.h:170
void minimize()
Definition: Controller.C:594
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
Tensor pressure_amd
Definition: Controller.h:71
Vector origin() const
Definition: Lattice.h:262
Tensor pressure_slow
Definition: Controller.h:70
NamdState *const state
Definition: Controller.h:236
Tensor pressure
Definition: Controller.h:167
Tensor groupPressure_slow
Definition: Controller.h:75
void NAMD_bug(const char *err_msg)
Definition: common.C:123
Tensor groupPressure_normal
Definition: Controller.h:73
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
BigReal getVirialTailCorr(const BigReal)
Tensor virial_amd
Definition: Controller.h:72
BigReal volume(void) const
Definition: Lattice.h:277
Tensor outer(const Vector &v1, const Vector &v2)
Definition: Tensor.h:241
Definition: Tensor.h:15
Tensor pressure_normal
Definition: Controller.h:68
Tensor groupPressure
Definition: Controller.h:168
BigReal getCurrentLambda(const int)
int slowFreq
Definition: Controller.h:80
Tensor controlPressure_normal
Definition: Controller.h:76
Molecule * molecule
Definition: Node.h:176
Tensor pressure_nbond
Definition: Controller.h:69
Tensor controlPressure_nbond
Definition: Controller.h:77
double BigReal
Definition: common.h:112
void Controller::compareChecksums ( int  step,
int  forgiving = 0 
)
protected

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

2767  {
2768  Node *node = Node::Object();
2769  Molecule *molecule = node->molecule;
2770 
2771  // Some basic correctness checking
2772  BigReal checksum, checksum_b;
2773  char errmsg[256];
2774 
2775  checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
2776  if ( ((int)checksum) != molecule->numAtoms ) {
2777  sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
2778  (int)checksum, molecule->numAtoms);
2779  NAMD_bug(errmsg);
2780  }
2781 
2783  if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
2784  if ( ! computeChecksum ) {
2785  computesPartitioned = 0;
2786  } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
2787  computesPartitioned = 1;
2788  } else {
2789  NAMD_bug("Bad global compute count!\n");
2790  }
2791  computeChecksum = ((int)checksum);
2792  }
2793 
2794  checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
2795  if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
2796  sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
2797  (int)checksum, molecule->numCalcBonds);
2798  if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
2799  iout << iWARN << errmsg << endi;
2800  else NAMD_bug(errmsg);
2801  }
2802 
2804  if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
2805  sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
2806  (int)checksum, molecule->numCalcAngles);
2807  if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
2808  iout << iWARN << errmsg << endi;
2809  else NAMD_bug(errmsg);
2810  }
2811 
2813  if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
2814  sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
2815  (int)checksum, molecule->numCalcDihedrals);
2816  if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
2817  iout << iWARN << errmsg << endi;
2818  else NAMD_bug(errmsg);
2819  }
2820 
2822  if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
2823  sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
2824  (int)checksum, molecule->numCalcImpropers);
2825  if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
2826  iout << iWARN << errmsg << endi;
2827  else NAMD_bug(errmsg);
2828  }
2829 
2831  if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
2832  sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
2833  (int)checksum, molecule->numCalcTholes);
2834  if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
2835  iout << iWARN << errmsg << endi;
2836  else NAMD_bug(errmsg);
2837  }
2838 
2840  if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
2841  sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
2842  (int)checksum, molecule->numCalcAnisos);
2843  if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
2844  iout << iWARN << errmsg << endi;
2845  else NAMD_bug(errmsg);
2846  }
2847 
2849  if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
2850  sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
2851  (int)checksum, molecule->numCalcCrossterms);
2852  if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
2853  iout << iWARN << errmsg << endi;
2854  else NAMD_bug(errmsg);
2855  }
2856 
2858  if ( ((int64)checksum) > molecule->numCalcExclusions &&
2859  ( ! simParams->mollyOn || step % slowFreq ) ) {
2860  if ( forgiving )
2861  iout << iWARN << "High global exclusion count ("
2862  << ((int64)checksum) << " vs "
2863  << molecule->numCalcExclusions << "), possible error!\n"
2864  << iWARN << "This warning is not unusual during minimization.\n"
2865  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
2866  else {
2867  char errmsg[256];
2868  sprintf(errmsg, "High global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2869  (long long)checksum, (long long)molecule->numCalcExclusions);
2870  NAMD_bug(errmsg);
2871  }
2872  }
2873 
2874  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
2876  if ( forgiving || ! simParams->fullElectFrequency ) {
2877  iout << iWARN << "Low global exclusion count! ("
2878  << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
2879  if ( forgiving ) {
2880  iout << "\n"
2881  << iWARN << "This warning is not unusual during minimization.\n"
2882  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
2883  } else {
2884  iout << " System unstable or pairlistdist or cutoff too small.\n";
2885  }
2886  iout << endi;
2887  } else {
2888  char errmsg[256];
2889  sprintf(errmsg, "Low global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2890  (long long)checksum, (long long)molecule->numCalcExclusions);
2891  NAMD_bug(errmsg);
2892  }
2893  }
2894 
2895 #ifdef NAMD_CUDA
2897  if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
2898  ( ! simParams->mollyOn || step % slowFreq ) ) {
2899  if ( forgiving )
2900  iout << iWARN << "High global CUDA exclusion count ("
2901  << ((int64)checksum) << " vs "
2902  << molecule->numCalcFullExclusions << "), possible error!\n"
2903  << iWARN << "This warning is not unusual during minimization.\n"
2904  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
2905  else {
2906  char errmsg[256];
2907  sprintf(errmsg, "High global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2908  (long long)checksum, (long long)molecule->numCalcFullExclusions);
2909  NAMD_bug(errmsg);
2910  }
2911  }
2912 
2913  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
2915  if ( forgiving || ! simParams->fullElectFrequency ) {
2916  iout << iWARN << "Low global CUDA exclusion count! ("
2917  << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
2918  if ( forgiving ) {
2919  iout << "\n"
2920  << iWARN << "This warning is not unusual during minimization.\n"
2921  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
2922  } else {
2923  iout << " System unstable or pairlistdist or cutoff too small.\n";
2924  }
2925  iout << endi;
2926  } else {
2927  char errmsg[256];
2928  sprintf(errmsg, "Low global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2929  (long long)checksum, (long long)molecule->numCalcFullExclusions);
2930  NAMD_bug(errmsg);
2931  }
2932  }
2933 #endif
2934 
2936  if ( ((int)checksum) && ! marginViolations ) {
2937  iout << iERROR << "Margin is too small for " << ((int)checksum) <<
2938  " atoms during timestep " << step << ".\n" << iERROR <<
2939  "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
2940  }
2941  marginViolations += (int)checksum;
2942 
2944  if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
2945  iout << iINFO <<
2946  "Pairlistdist is too small for " << ((int)checksum) <<
2947  " computes during timestep " << step << ".\n" << endi;
2948  }
2949  if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
2950 
2952  if ( checksum ) {
2953  if ( forgiving )
2954  iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
2955  else NAMD_bug("Stray PME grid charges detected!\n");
2956  }
2957 }
static Node * Object()
Definition: Node.h:86
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
int numCalcBonds
Definition: Molecule.h:621
int numCalcAnisos
Definition: Molecule.h:631
Definition: Node.h:78
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
#define iout
Definition: InfoStream.h:87
int numCalcCrossterms
Definition: Molecule.h:625
int64 numCalcFullExclusions
Definition: Molecule.h:627
void NAMD_bug(const char *err_msg)
Definition: common.C:123
int computeChecksum
Definition: Controller.h:89
int marginViolations
Definition: Controller.h:90
int numCalcDihedrals
Definition: Molecule.h:623
int numCalcImpropers
Definition: Molecule.h:624
int numAtoms
Definition: Molecule.h:556
BigReal item(int i) const
Definition: ReductionMgr.h:340
int64 numCalcExclusions
Definition: Molecule.h:626
int pairlistWarnings
Definition: Controller.h:91
long long int64
Definition: common.h:34
int slowFreq
Definition: Controller.h:80
SimParameters *const simParams
Definition: Controller.h:235
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
RequireReduction * reduction
Definition: Controller.h:237
int numCalcTholes
Definition: Molecule.h:630
Molecule * molecule
Definition: Node.h:176
int numCalcAngles
Definition: Molecule.h:622
double BigReal
Definition: common.h:112
BigReal Controller::computeAlchWork ( const int  step)
protected

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

3870  {
3871  // alchemical scaling factors for groups 1/2 at the previous lambda
3872  const BigReal oldLambda = simParams->getCurrentLambda(step-1);
3873  const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
3874  const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
3875  const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
3876  const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
3877  const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
3878  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
3879  // alchemical scaling factors for groups 1/2 at the new lambda
3880  const BigReal alchLambda = simParams->getCurrentLambda(step);
3881  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
3882  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
3883  const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
3884  const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
3885  const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
3886  const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda);
3887 
3888  return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
3889  (elec_lambda_12 - elec_lambda_1)*
3891  (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
3892  (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
3893  (elec_lambda_22 - elec_lambda_2)*
3895  (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
3896  );
3897 }
BigReal electEnergy_ti_1
Definition: Controller.h:128
BigReal ljEnergy_ti_1
Definition: Controller.h:130
BigReal getBondLambda(const BigReal)
BigReal electEnergySlow_ti_1
Definition: Controller.h:129
BigReal electEnergyPME_ti_1
Definition: Controller.h:141
BigReal electEnergy_ti_2
Definition: Controller.h:131
BigReal electEnergyPME_ti_2
Definition: Controller.h:142
BigReal getCurrentLambda(const int)
BigReal electEnergySlow_ti_2
Definition: Controller.h:132
BigReal getVdwLambda(const BigReal)
SimParameters *const simParams
Definition: Controller.h:235
BigReal bondedEnergy_ti_2
Definition: Controller.h:127
BigReal bondedEnergy_ti_1
Definition: Controller.h:126
BigReal getElecLambda(const BigReal)
double BigReal
Definition: common.h:112
BigReal ljEnergy_ti_2
Definition: Controller.h:133
void Controller::correctMomentum ( int  step)
protected

Definition at line 1255 of file Controller.C.

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

Referenced by integrate().

1255  {
1256 
1257  Vector momentum;
1258  momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
1259  momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
1260  momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
1262 
1263  if ( momentum.length2() == 0. )
1264  iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
1265  if ( mass == 0. )
1266  NAMD_die("Total mass is zero in Controller::correctMomentum");
1267 
1268  momentum *= (-1./mass);
1269 
1270  broadcast->momentumCorrection.publish(step+slowFreq,momentum);
1271 }
SimpleBroadcastObject< Vector > momentumCorrection
Definition: Broadcasts.h:79
Definition: Vector.h:64
BigReal z
Definition: Vector.h:66
#define iout
Definition: InfoStream.h:87
ControllerBroadcasts * broadcast
Definition: Controller.h:251
BigReal x
Definition: Vector.h:66
void NAMD_die(const char *err_msg)
Definition: common.C:83
void publish(int tag, const T &t)
BigReal item(int i) const
Definition: ReductionMgr.h:340
BigReal length2(void) const
Definition: Vector.h:173
BigReal y
Definition: Vector.h:66
int slowFreq
Definition: Controller.h:80
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
RequireReduction * reduction
Definition: Controller.h:237
double BigReal
Definition: common.h:112
void Controller::cycleBarrier ( int  doBarrier,
int  step 
)
protected

Definition at line 4131 of file Controller.C.

References broadcast.

Referenced by integrate().

4131  {
4132 #if USE_BARRIER
4133  if (doBarrier) {
4134  broadcast->cycleBarrier.publish(step,1);
4135  CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
4136  CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4137  }
4138 #endif
4139 }
ControllerBroadcasts * broadcast
Definition: Controller.h:251
void Controller::enqueueCollections ( int  timestep)
protected

Definition at line 3655 of file Controller.C.

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

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

3656 {
3657  if ( Output::coordinateNeeded(timestep) ){
3658  collection->enqueuePositions(timestep,state->lattice);
3659  }
3660  if ( Output::velocityNeeded(timestep) )
3661  collection->enqueueVelocities(timestep);
3662  if ( Output::forceNeeded(timestep) )
3663  collection->enqueueForces(timestep);
3664 }
static int velocityNeeded(int)
Definition: Output.C:365
static int coordinateNeeded(int)
Definition: Output.C:198
void enqueueVelocities(int seq)
void enqueuePositions(int seq, Lattice &lattice)
static int forceNeeded(int)
Definition: Output.C:460
NamdState *const state
Definition: Controller.h:236
CollectionMaster *const collection
Definition: Controller.h:249
void enqueueForces(int seq)
void Controller::integrate ( int  scriptTask)
protected

Definition at line 429 of file Controller.C.

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

Referenced by algorithm().

429  {
430  char traceNote[24];
431 
432  int step = simParams->firstTimestep;
433 
434  const int numberOfSteps = simParams->N;
435  const int stepsPerCycle = simParams->stepsPerCycle;
436 
437  const int zeroMomentum = simParams->zeroMomentum;
438 
440  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
441  if (dofull)
443  else
445  if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
446 
447  if ( scriptTask == SCRIPT_RUN ) {
448 
449  reassignVelocities(step); // only for full-step velecities
450  rescaleaccelMD(step);
451  adaptTempUpdate(step); // Init adaptive tempering;
452 
453  receivePressure(step);
454  if ( zeroMomentum && dofull && ! (step % slowFreq) )
455  correctMomentum(step);
456  printFepMessage(step);
457  printTiMessage(step);
458  printDynamicsEnergies(step);
459  outputFepEnergy(step);
460  outputTiEnergy(step);
461  if(traceIsOn()){
462  traceUserEvent(eventEndOfTimeStep);
463  sprintf(traceNote, "s:%d", step);
464  traceUserSuppliedNote(traceNote);
465  }
466  outputExtendedSystem(step);
467  rebalanceLoad(step);
468 
469  }
470 
471  // Handling SIGINT doesn't seem to be working on Lemieux, and it
472  // sometimes causes the net-xxx versions of NAMD to segfault on exit,
473  // so disable it for now.
474  // namd_sighandler_t oldhandler = signal(SIGINT,
475  // (namd_sighandler_t)my_sigint_handler);
476  for ( ++step ; step <= numberOfSteps; ++step )
477  {
478 
479  adaptTempUpdate(step);
480  rescaleVelocities(step);
481  tcoupleVelocities(step);
483  berendsenPressure(step);
484  langevinPiston1(step);
485  rescaleaccelMD(step);
486  enqueueCollections(step); // after lattice scaling!
487  receivePressure(step);
488  if ( zeroMomentum && dofull && ! (step % slowFreq) )
489  correctMomentum(step);
490  langevinPiston2(step);
491  reassignVelocities(step);
492 
493  multigratorTemperature(step, 1);
494  multigratorPressure(step, 1);
495  multigratorPressure(step, 2);
496  multigratorTemperature(step, 2);
497 
498  printDynamicsEnergies(step);
499  outputFepEnergy(step);
500  outputTiEnergy(step);
501  if(traceIsOn()){
502  traceUserEvent(eventEndOfTimeStep);
503  sprintf(traceNote, "s:%d", step);
504  traceUserSuppliedNote(traceNote);
505  }
506  // if (gotsigint) {
507  // iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
508  // NAMD_quit();
509  // }
510  outputExtendedSystem(step);
511 #if CYCLE_BARRIER
512  cycleBarrier(!((step+1) % stepsPerCycle),step);
513 #elif PME_BARRIER
514  cycleBarrier(dofull && !(step%slowFreq),step);
515 #elif STEP_BARRIER
516  cycleBarrier(1, step);
517 #endif
518 
519  if(Node::Object()->specialTracing || simParams->statsOn){
520  int bstep = simParams->traceStartStep;
521  int estep = bstep + simParams->numTraceSteps;
522  if(step == bstep){
523  traceBarrier(1, step);
524  }
525  if(step == estep){
526  traceBarrier(0, step);
527  }
528  }
529 
530 #ifdef MEASURE_NAMD_WITH_PAPI
531  if(simParams->papiMeasure) {
532  int bstep = simParams->papiMeasureStartStep;
533  int estep = bstep + simParams->numPapiMeasureSteps;
534  if(step == bstep) {
535  papiMeasureBarrier(1, step);
536  }
537  if(step == estep) {
538  papiMeasureBarrier(0, step);
539  }
540  }
541 #endif
542 
543  rebalanceLoad(step);
544 
545 #if PME_BARRIER
546  cycleBarrier(dofull && !((step+1)%slowFreq),step); // step before PME
547 #endif
548  }
549  // signal(SIGINT, oldhandler);
550 }
static Node * Object()
Definition: Node.h:86
void enqueueCollections(int)
Definition: Controller.C:3655
void cycleBarrier(int, int)
Definition: Controller.C:4131
void rescaleVelocities(int)
Definition: Controller.C:1233
void rescaleaccelMD(int step, int minimize=0)
Definition: Controller.C:1831
int nbondFreq
Definition: Controller.h:79
int eventEndOfTimeStep
Definition: Node.C:285
void printTiMessage(int)
Definition: Controller.C:1294
void multigratorTemperature(int step, int callNumber)
Definition: Controller.C:853
void langevinPiston1(int)
Definition: Controller.C:987
void outputTiEnergy(int step)
Definition: Controller.C:3741
void stochRescaleVelocities(int)
Definition: Controller.C:1348
void correctMomentum(int step)
Definition: Controller.C:1255
void berendsenPressure(int)
Definition: Controller.C:950
void langevinPiston2(int)
Definition: Controller.C:1142
void receivePressure(int step, int minimize=0)
Definition: Controller.C:1412
void reassignVelocities(int)
Definition: Controller.C:1308
void outputFepEnergy(int step)
Definition: Controller.C:3667
void printFepMessage(int)
Definition: Controller.C:1274
int slowFreq
Definition: Controller.h:80
void traceBarrier(int, int)
Definition: Controller.C:4141
void printDynamicsEnergies(int)
Definition: Controller.C:3015
void multigratorPressure(int step, int callNumber)
Definition: Controller.C:778
void tcoupleVelocities(int)
Definition: Controller.C:1326
SimParameters *const simParams
Definition: Controller.h:235
void rebalanceLoad(int)
Definition: Controller.C:4117
void outputExtendedSystem(int step)
Definition: Controller.C:3959
void adaptTempUpdate(int step, int minimize=0)
Definition: Controller.C:2469
void Controller::langevinPiston1 ( int  step)
protected

Definition at line 987 of file Controller.C.

References BOLTZMANN, broadcast, Lattice::c(), controlNumDegFreedom, controlPressure, controlPressure_nbond, controlPressure_slow, Tensor::diagonal(), SimParameters::dt, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, Random::gaussian(), Random::gaussian_vector(), Tensor::identity(), iINFO(), iout, ControllerState::langevinPiston_strainRate, SimParameters::langevinPistonBarrier, SimParameters::langevinPistonDecay, SimParameters::langevinPistonOn, SimParameters::langevinPistonPeriod, SimParameters::langevinPistonTarget, SimParameters::langevinPistonTemp, nbondFreq, ControllerBroadcasts::positionRescaleFactor, 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().

988 {
990  {
991  Tensor &strainRate = langevinPiston_strainRate;
992  int cellDims = simParams->useFlexibleCell ? 1 : 3;
993  BigReal dt = simParams->dt;
994  BigReal dt_long = slowFreq * dt;
997  BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
998 
999 #ifdef DEBUG_PRESSURE
1000  iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
1001 #endif
1002 
1003  if ( ! ( (step-1) % slowFreq ) )
1004  {
1005  BigReal gamma = 1 / simParams->langevinPistonDecay;
1006  BigReal f1 = exp( -0.5 * dt_long * gamma );
1007  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1008  strainRate *= f1;
1009  if ( simParams->useFlexibleCell ) {
1010  // We only use on-diagonal terms (for now)
1011  if ( simParams->useConstantRatio ) {
1012  BigReal r = f2 * random->gaussian();
1013  strainRate.xx += r;
1014  strainRate.yy += r;
1015  strainRate.zz += f2 * random->gaussian();
1016  } else {
1017  strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
1018  }
1019  } else {
1020  strainRate += f2 * Tensor::identity(random->gaussian());
1021  }
1022 
1023 #ifdef DEBUG_PRESSURE
1024  iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
1025 #endif
1026  }
1027 
1028  // Apply surface tension. If surfaceTensionTarget is zero, we get
1029  // the default (isotropic pressure) case.
1030 
1031  Tensor ptarget;
1032  ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
1033  if ( simParams->surfaceTensionTarget != 0. ) {
1034  ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
1035  simParams->surfaceTensionTarget / state->lattice.c().z;
1036  }
1037 
1038  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1039  ( controlPressure - ptarget );
1040 
1041  if (simParams->fixCellDims) {
1042  if (simParams->fixCellDimX) strainRate.xx = 0;
1043  if (simParams->fixCellDimY) strainRate.yy = 0;
1044  if (simParams->fixCellDimZ) strainRate.zz = 0;
1045  }
1046 
1047 
1048 #ifdef DEBUG_PRESSURE
1049  iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
1050 #endif
1051 
1053  if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
1054  {
1055  // We only use on-diagonal terms (for now)
1056  Tensor factor;
1057  if ( !simParams->useConstantArea ) {
1058  factor.xx = exp( dt_long * strainRate.xx );
1059  factor.yy = exp( dt_long * strainRate.yy );
1060  } else {
1061  factor.xx = factor.yy = 1;
1062  }
1063  factor.zz = exp( dt_long * strainRate.zz );
1064  broadcast->positionRescaleFactor.publish(step,factor);
1065  state->lattice.rescale(factor);
1066 #ifdef DEBUG_PRESSURE
1067  iout << iINFO << "rescaling by: " << factor << "\n";
1068 #endif
1069  }
1070  } else { // langevinPistonBarrier
1071  if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
1072  {
1073  if ( ! positionRescaleFactor.xx ) { // first pass without MTS
1074  // We only use on-diagonal terms (for now)
1075  Tensor factor;
1076  if ( !simParams->useConstantArea ) {
1077  factor.xx = exp( dt_long * strainRate.xx );
1078  factor.yy = exp( dt_long * strainRate.yy );
1079  } else {
1080  factor.xx = factor.yy = 1;
1081  }
1082  factor.zz = exp( dt_long * strainRate.zz );
1083  broadcast->positionRescaleFactor.publish(step,factor);
1084  positionRescaleFactor = factor;
1085  strainRate_old = strainRate;
1086  }
1088 #ifdef DEBUG_PRESSURE
1089  iout << iINFO << "rescaling by: " << positionRescaleFactor << "\n";
1090 #endif
1091  }
1092  if ( ! ( (step-slowFreq/2) % slowFreq ) )
1093  {
1094  // We only use on-diagonal terms (for now)
1095  Tensor factor;
1096  if ( !simParams->useConstantArea ) {
1097  factor.xx = exp( (dt+dt_long) * strainRate.xx - dt * strainRate_old.xx );
1098  factor.yy = exp( (dt+dt_long) * strainRate.yy - dt * strainRate_old.yy );
1099  } else {
1100  factor.xx = factor.yy = 1;
1101  }
1102  factor.zz = exp( (dt+dt_long) * strainRate.zz - dt * strainRate_old.zz );
1103  broadcast->positionRescaleFactor.publish(step+1,factor);
1104  positionRescaleFactor = factor;
1105  strainRate_old = strainRate;
1106  }
1107  }
1108 
1109  // corrections to integrator
1110  if ( ! ( step % nbondFreq ) )
1111  {
1112 #ifdef DEBUG_PRESSURE
1113  iout << iINFO << "correcting strain rate for nbond, ";
1114 #endif
1115  strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1116  ( (nbondFreq - 1) * controlPressure_nbond );
1117 #ifdef DEBUG_PRESSURE
1118  iout << "strain rate: " << strainRate << "\n";
1119 #endif
1120  }
1121  if ( ! ( step % slowFreq ) )
1122  {
1123 #ifdef DEBUG_PRESSURE
1124  iout << iINFO << "correcting strain rate for slow, ";
1125 #endif
1126  strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1127  ( (slowFreq - 1) * controlPressure_slow );
1128 #ifdef DEBUG_PRESSURE
1129  iout << "strain rate: " << strainRate << "\n";
1130 #endif
1131  }
1132  if (simParams->fixCellDims) {
1133  if (simParams->fixCellDimX) strainRate.xx = 0;
1134  if (simParams->fixCellDimY) strainRate.yy = 0;
1135  if (simParams->fixCellDimZ) strainRate.zz = 0;
1136  }
1137 
1138  }
1139 
1140 }
Vector gaussian_vector(void)
Definition: Random.h:208
Tensor controlPressure_slow
Definition: Controller.h:78
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
int nbondFreq
Definition: Controller.h:79
BigReal langevinPistonTemp
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
#define BOLTZMANN
Definition: common.h:45
BigReal gaussian(void)
Definition: Random.h:116
Tensor controlPressure
Definition: Controller.h:170
BigReal surfaceTensionTarget
BigReal z
Definition: Vector.h:66
Bool langevinPistonBarrier
#define iout
Definition: InfoStream.h:87
Tensor langevinPiston_strainRate
Definition: Controller.h:33
BigReal langevinPistonDecay
NamdState *const state
Definition: Controller.h:236
void rescale(Tensor factor)
Definition: Lattice.h:60
ControllerBroadcasts * broadcast
Definition: Controller.h:251
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
void publish(int tag, const T &t)
BigReal volume(void) const
Definition: Lattice.h:277
Random * random
Definition: Controller.h:234
BigReal xx
Definition: Tensor.h:17
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
BigReal langevinPistonPeriod
BigReal zz
Definition: Tensor.h:19
Definition: Tensor.h:15
Tensor strainRate_old
Definition: Controller.h:205
int slowFreq
Definition: Controller.h:80
BigReal yy
Definition: Tensor.h:18
SimParameters *const simParams
Definition: Controller.h:235
Tensor positionRescaleFactor
Definition: Controller.h:206
BigReal langevinPistonTarget
int controlNumDegFreedom
Definition: Controller.h:169
Vector c() const
Definition: Lattice.h:254
Tensor controlPressure_nbond
Definition: Controller.h:77
double BigReal
Definition: common.h:112
void Controller::langevinPiston2 ( int  step)
protected

Definition at line 1142 of file Controller.C.

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

Referenced by integrate().

1143 {
1144  if ( simParams->langevinPistonOn )
1145  {
1146  Tensor &strainRate = langevinPiston_strainRate;
1147  int cellDims = simParams->useFlexibleCell ? 1 : 3;
1148  BigReal dt = simParams->dt;
1149  BigReal dt_long = slowFreq * dt;
1152  BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
1153 
1154  // corrections to integrator
1155  if ( ! ( step % nbondFreq ) )
1156  {
1157 #ifdef DEBUG_PRESSURE
1158  iout << iINFO << "correcting strain rate for nbond, ";
1159 #endif
1160  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1161  ( (nbondFreq - 1) * controlPressure_nbond );
1162 #ifdef DEBUG_PRESSURE
1163  iout << "strain rate: " << strainRate << "\n";
1164 #endif
1165  }
1166  if ( ! ( step % slowFreq ) )
1167  {
1168 #ifdef DEBUG_PRESSURE
1169  iout << iINFO << "correcting strain rate for slow, ";
1170 #endif
1171  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1172  ( (slowFreq - 1) * controlPressure_slow );
1173 #ifdef DEBUG_PRESSURE
1174  iout << "strain rate: " << strainRate << "\n";
1175 #endif
1176  }
1177 
1178  // Apply surface tension. If surfaceTensionTarget is zero, we get
1179  // the default (isotropic pressure) case.
1180 
1181  Tensor ptarget;
1182  ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
1183  if ( simParams->surfaceTensionTarget != 0. ) {
1184  ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
1185  simParams->surfaceTensionTarget / state->lattice.c().z;
1186  }
1187 
1188  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1189  ( controlPressure - ptarget );
1190 
1191 
1192 #ifdef DEBUG_PRESSURE
1193  iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
1194 #endif
1195 
1196  if ( ! ( step % slowFreq ) )
1197  {
1198  BigReal gamma = 1 / simParams->langevinPistonDecay;
1199  BigReal f1 = exp( -0.5 * dt_long * gamma );
1200  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1201  strainRate *= f1;
1202  if ( simParams->useFlexibleCell ) {
1203  // We only use on-diagonal terms (for now)
1204  if ( simParams->useConstantRatio ) {
1205  BigReal r = f2 * random->gaussian();
1206  strainRate.xx += r;
1207  strainRate.yy += r;
1208  strainRate.zz += f2 * random->gaussian();
1209  } else {
1210  strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
1211  }
1212  } else {
1213  strainRate += f2 * Tensor::identity(random->gaussian());
1214  }
1215 #ifdef DEBUG_PRESSURE
1216  iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
1217 #endif
1218  }
1219 
1220 #ifdef DEBUG_PRESSURE
1221  iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
1222 #endif
1223  if (simParams->fixCellDims) {
1224  if (simParams->fixCellDimX) strainRate.xx = 0;
1225  if (simParams->fixCellDimY) strainRate.yy = 0;
1226  if (simParams->fixCellDimZ) strainRate.zz = 0;
1227  }
1228  }
1229 
1230 
1231 }
Vector gaussian_vector(void)
Definition: Random.h:208
Tensor controlPressure_slow
Definition: Controller.h:78
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
int nbondFreq
Definition: Controller.h:79
BigReal langevinPistonTemp
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
#define BOLTZMANN
Definition: common.h:45
BigReal gaussian(void)
Definition: Random.h:116
Tensor controlPressure
Definition: Controller.h:170
BigReal surfaceTensionTarget
BigReal z
Definition: Vector.h:66
#define iout
Definition: InfoStream.h:87
Tensor langevinPiston_strainRate
Definition: Controller.h:33
BigReal langevinPistonDecay
NamdState *const state
Definition: Controller.h:236
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
BigReal volume(void) const
Definition: Lattice.h:277
Random * random
Definition: Controller.h:234
BigReal xx
Definition: Tensor.h:17
BigReal langevinPistonPeriod
BigReal zz
Definition: Tensor.h:19
Definition: Tensor.h:15
int slowFreq
Definition: Controller.h:80
BigReal yy
Definition: Tensor.h:18
SimParameters *const simParams
Definition: Controller.h:235
BigReal langevinPistonTarget
int controlNumDegFreedom
Definition: Controller.h:169
Vector c() const
Definition: Lattice.h:254
Tensor controlPressure_nbond
Definition: Controller.h:77
double BigReal
Definition: common.h:112
void Controller::minimize ( )
protected

Definition at line 594 of file Controller.C.

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

Referenced by algorithm().

594  {
595  // iout << "Controller::minimize() called.\n" << endi;
596 
597  const int minVerbose = simParams->minVerbose;
598  const int numberOfSteps = simParams->N;
599  int step = simParams->firstTimestep;
600  slowFreq = nbondFreq = 1;
601  BigReal linegoal = simParams->minLineGoal; // 1.0e-3
602  const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
603 
604  CALCULATE
605 
606  int minSeq = 0;
607 
608  // just move downhill until initial bad contacts go away or 100 steps
609  int old_min_huge_count = 2000000000;
610  int steps_at_same_min_huge_count = 0;
611  for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
612  if ( min_huge_count >= old_min_huge_count ) {
613  if ( ++steps_at_same_min_huge_count > 10 ) break;
614  } else {
615  old_min_huge_count = min_huge_count;
616  steps_at_same_min_huge_count = 0;
617  }
618  iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
619  broadcast->minimizeCoefficient.publish(minSeq++,1.);
620  enqueueCollections(step);
621  CALCULATE
622  }
623  if ( min_huge_count ) {
624  iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
625  }
626  iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
627 
628  int atStart = 2;
629  int errorFactor = 10;
630  BigReal old_f_dot_f = min_f_dot_f;
631  broadcast->minimizeCoefficient.publish(minSeq++,0.);
632  broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
633  int newDir = 1;
636  while ( 1 ) {
637  // line minimization
638  // bracket minimum on line
639  minpoint lo,hi,mid,last;
640  BigReal x = 0;
641  lo.x = x;
642  lo.u = min_energy;
643  lo.dudx = -1. * min_f_dot_v;
645  mid = lo;
646  last = mid;
647  if ( minVerbose ) {
648  iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
649  if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
650  iout << "\n" << endi;
651  }
652  BigReal tol = fabs( linegoal * min_f_dot_v );
653  iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
654  fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
655  int start_with_huge = last.noGradient;
657  BigReal maxstep = 0.1 / sqrt(min_reduction->item(0));
658  x = maxstep; MOVETO(x);
659  // bracket minimum on line
660  while ( last.u < mid.u ) {
661  if ( last.dudx < mid.dudx * (5.5 - x/maxstep)/5. ) {
662  x = 6 * maxstep; break;
663  }
664  lo = mid; mid = last;
665  x += maxstep;
666  if ( x > 5.5 * maxstep ) break;
667  MOVETO(x)
668  }
669  if ( x > 5.5 * maxstep ) {
670  iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" << endi;
671  broadcast->minimizeCoefficient.publish(minSeq++,0.);
672  broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
673  newDir = 1;
674  old_f_dot_f = min_f_dot_f;
677  continue;
678  }
679  hi = last;
680 #define PRINT_BRACKET \
681  iout << "LINE MINIMIZER BRACKET: DX " \
682  << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
683  " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
684  lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
686  // converge on minimum on line
687  int itcnt;
688  for ( itcnt = 10; itcnt > 0; --itcnt ) {
689  int progress = 1;
690  // select new position
691  if ( mid.noGradient ) {
692  if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) { // subdivide left side
693  x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
694  MOVETO(x)
695  if ( last.u <= mid.u ) {
696  hi = mid; mid = last;
697  } else {
698  lo = last;
699  }
700  } else { // subdivide right side
701  x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
702  MOVETO(x)
703  if ( last.u <= mid.u ) {
704  lo = mid; mid = last;
705  } else {
706  hi = last;
707  }
708  }
709  } else {
710  if ( mid.dudx > 0. ) { // subdivide left side
711  BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
712  BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
713  x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
714  BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
715  if ( xdiv ) x /= xdiv; else x = last.x;
716  if ( x > altxhi ) x = altxhi;
717  if ( x < altxlo ) x = altxlo;
718  if ( x-last.x == 0 ) break;
719  MOVETO(x)
720  if ( last.u <= mid.u ) {
721  hi = mid; mid = last;
722  } else {
723  if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
724  lo = last;
725  }
726  } else { // subdivide right side
727  BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
728  BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
729  x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
730  BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
731  if ( xdiv ) x /= xdiv; else x = last.x;
732  if ( x < altxlo ) x = altxlo;
733  if ( x > altxhi ) x = altxhi;
734  if ( x-last.x == 0 ) break;
735  MOVETO(x)
736  if ( last.u <= mid.u ) {
737  lo = mid; mid = last;
738  } else {
739  if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
740  hi = last;
741  }
742  }
743  }
745  if ( ! progress ) {
746  MOVETO(mid.x);
747  break;
748  }
749  if ( fabs(last.dudx) < tol ) break;
750  if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
751  if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
752  }
753  // new direction
754  broadcast->minimizeCoefficient.publish(minSeq++,0.);
755  BigReal c = min_f_dot_f / old_f_dot_f;
756  c = ( c > 1.5 ? 1.5 : c );
757  if ( atStart ) { c = 0; --atStart; }
758  if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
759  c = 0;
760  if ( errorFactor < 100 ) errorFactor += 10;
761  }
762  if ( c == 0 ) {
763  iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
764  }
765  broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
766  newDir = 1;
767  old_f_dot_f = min_f_dot_f;
770  }
771 
772 }
void enqueueCollections(int)
Definition: Controller.C:3655
int nbondFreq
Definition: Controller.h:79
BigReal min_v_dot_v
Definition: Controller.h:97
BigReal noGradient
Definition: Controller.C:590
int min_huge_count
Definition: Controller.h:98
BigReal minLineGoal
if(ComputeNonbondedUtil::goMethod==2)
BigReal u
Definition: Controller.C:590
#define iout
Definition: InfoStream.h:87
BigReal min_energy
Definition: Controller.h:94
#define MOVETO(X)
Definition: Controller.C:560
void require(void)
Definition: ReductionMgr.h:341
ControllerBroadcasts * broadcast
Definition: Controller.h:251
BigReal x
Definition: Controller.C:590
BigReal min_f_dot_f
Definition: Controller.h:95
void publish(int tag, const T &t)
SimpleBroadcastObject< BigReal > minimizeCoefficient
Definition: Broadcasts.h:78
BigReal item(int i) const
Definition: ReductionMgr.h:340
BigReal dudx
Definition: Controller.C:590
RequireReduction * min_reduction
Definition: Controller.h:60
#define CALCULATE
Definition: Controller.C:553
#define PRINT_BRACKET
int slowFreq
Definition: Controller.h:80
SimParameters *const simParams
Definition: Controller.h:235
infostream & endi(infostream &s)
Definition: InfoStream.C:38
gridSize x
double BigReal
Definition: common.h:112
BigReal min_f_dot_v
Definition: Controller.h:96
BigReal Controller::multigatorCalcEnthalpy ( BigReal  potentialEnergy,
int  step,
int  minimize 
)
protected

Definition at line 919 of file Controller.C.

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

Referenced by printEnergies().

919  {
920  Node *node = Node::Object();
921  Molecule *molecule = node->molecule;
922  SimParameters *simParameters = node->simParameters;
923 
924  BigReal V = state->lattice.volume();
928  BigReal Nf = numDegFreedom;
930  BigReal sumZeta = 0.0;
931  for (int i=1;i < simParams->multigratorNoseHooverChainLength;i++) {
932  sumZeta += multigratorZeta[i];
933  }
934  BigReal nuOmegaSum = 0.0;
935  for (int i=0;i < simParams->multigratorNoseHooverChainLength;i++) {
936  nuOmegaSum += multigratorNu[i]*multigratorNu[i]/(2.0*multigratorOmega[i]);
937  }
938  BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
939  BigReal eta = sqrt(kT0*W)*multigratorXi;
940 
941  BigReal enthalpy = kineticEnergy + potentialEnergy + eta*eta/(2.0*W) + P0*V + nuOmegaSum + kT0*(Nf*multigratorZeta[0] + sumZeta);
942 
943 // if (!(step % 100))
944  // fprintf(stderr, "enthalpy %lf %lf %lf %lf %lf %lf %lf\n", enthalpy,
945  // kineticEnergy, potentialEnergy, eta*eta/(2.0*W), P0*V, nuOmegaSum, kT0*(Nf*multigratorZeta[0] + sumZeta));
946 
947  return enthalpy;
948 }
static Node * Object()
Definition: Node.h:86
#define BOLTZMANN
Definition: common.h:45
Definition: Node.h:78
SimParameters * simParameters
Definition: Node.h:178
BigReal multigratorPressureTarget
std::vector< BigReal > multigratorOmega
Definition: Controller.h:215
std::vector< BigReal > multigratorNu
Definition: Controller.h:213
NamdState *const state
Definition: Controller.h:236
int64_t numDegFreedom
Definition: Controller.h:101
std::vector< BigReal > multigratorZeta
Definition: Controller.h:216
int multigratorNoseHooverChainLength
BigReal multigratorTemperatureTarget
BigReal multigratorPressureRelaxationTime
BigReal volume(void) const
Definition: Lattice.h:277
SimParameters *const simParams
Definition: Controller.h:235
Molecule * molecule
Definition: Node.h:176
int controlNumDegFreedom
Definition: Controller.h:169
BigReal multigratorXi
Definition: Controller.h:209
BigReal kineticEnergy
Definition: Controller.h:158
double BigReal
Definition: common.h:112
void Controller::multigratorPressure ( int  step,
int  callNumber 
)
protected

Definition at line 778 of file Controller.C.

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

Referenced by integrate().

778  {
783  BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
786  {
787  // Compute new scaling factors and send them to Sequencer
788  BigReal V = state->lattice.volume();
789  BigReal Pinst = trace(controlPressure)/3.0;
790  BigReal PGsum = trace(momentumSqrSum);
791  //
792  multigratorXiT = multigratorXi + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
793  BigReal scale = exp(s*NGfac*multigratorXiT);
794  BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
795  // fprintf(stderr, "%d | T %lf P %lf V %1.3lf\n", step, temperature, Pinst*PRESSUREFACTOR, V);
796  Tensor scaleTensor = Tensor::identity(scale);
797  Tensor volScaleTensor = Tensor::identity(scale);
798  Tensor velScaleTensor = Tensor::identity(velScale);
799  state->lattice.rescale(volScaleTensor);
800  if (callNumber == 1) {
801  broadcast->positionRescaleFactor.publish(step,scaleTensor);
802  broadcast->velocityRescaleTensor.publish(step,velScaleTensor);
803  } else {
804  broadcast->positionRescaleFactor2.publish(step,scaleTensor);
805  broadcast->velocityRescaleTensor2.publish(step,velScaleTensor);
806  }
807  }
808 
809  {
810  // Wait here for Sequencer to finish scaling and force computation
811  reduction->require();
812  Tensor virial_normal;
813  Tensor virial_nbond;
814  Tensor virial_slow;
815  Tensor intVirial_normal;
816  Tensor intVirial_nbond;
817  Tensor intVirial_slow;
818  Vector extForce_normal;
819  Vector extForce_nbond;
820  Vector extForce_slow;
821  GET_TENSOR(momentumSqrSum, reduction, REDUCTION_MOMENTUM_SQUARED);
822  GET_TENSOR(virial_normal, reduction, REDUCTION_VIRIAL_NORMAL);
823  GET_TENSOR(virial_nbond, reduction, REDUCTION_VIRIAL_NBOND);
824  GET_TENSOR(virial_slow, reduction, REDUCTION_VIRIAL_SLOW);
825  GET_TENSOR(intVirial_normal, reduction, REDUCTION_INT_VIRIAL_NORMAL);
826  GET_TENSOR(intVirial_nbond, reduction, REDUCTION_INT_VIRIAL_NBOND);
827  GET_TENSOR(intVirial_slow, reduction, REDUCTION_INT_VIRIAL_SLOW);
828  GET_VECTOR(extForce_normal, reduction, REDUCTION_EXT_FORCE_NORMAL);
829  GET_VECTOR(extForce_nbond, reduction, REDUCTION_EXT_FORCE_NBOND);
830  GET_VECTOR(extForce_slow, reduction, REDUCTION_EXT_FORCE_SLOW);
831  calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
832  intVirial_normal, intVirial_nbond, intVirial_slow,
833  extForce_normal, extForce_nbond, extForce_slow);
834  if (callNumber == 2) {
835  // Update temperature for the Temperature Cycle that is coming next
837  temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
838  }
839  }
840 
841  {
842  // Update pressure integrator
843  BigReal V = state->lattice.volume();
844  BigReal Pinst = trace(controlPressure)/3.0;
845  BigReal PGsum = trace(momentumSqrSum);
846  //
847  multigratorXi = multigratorXiT + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
848  }
849 
850  }
851 }
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)
Definition: Controller.C:1590
#define BOLTZMANN
Definition: common.h:45
Tensor controlPressure
Definition: Controller.h:170
Definition: Vector.h:64
BigReal multigratorPressureTarget
#define GET_TENSOR(O, R, A)
Definition: ReductionMgr.h:59
BigReal temperature
Definition: Controller.h:161
NamdState *const state
Definition: Controller.h:236
void rescale(Tensor factor)
Definition: Lattice.h:60
void require(void)
Definition: ReductionMgr.h:341
ControllerBroadcasts * broadcast
Definition: Controller.h:251
int64_t numDegFreedom
Definition: Controller.h:101
SimpleBroadcastObject< Tensor > velocityRescaleTensor2
Definition: Broadcasts.h:72
BigReal multigratorTemperatureTarget
#define GET_VECTOR(O, R, A)
Definition: ReductionMgr.h:54
BigReal multigratorPressureRelaxationTime
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
void publish(int tag, const T &t)
BigReal volume(void) const
Definition: Lattice.h:277
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
BigReal item(int i) const
Definition: ReductionMgr.h:340
Definition: Tensor.h:15
int multigratorPressureFreq
SimParameters *const simParams
Definition: Controller.h:235
RequireReduction * reduction
Definition: Controller.h:237
SimpleBroadcastObject< Tensor > positionRescaleFactor2
Definition: Broadcasts.h:74
SimpleBroadcastObject< Tensor > velocityRescaleTensor
Definition: Broadcasts.h:71
int controlNumDegFreedom
Definition: Controller.h:169
BigReal multigratorXiT
Definition: Controller.h:210
BigReal multigratorXi
Definition: Controller.h:209
BigReal kineticEnergy
Definition: Controller.h:158
Tensor momentumSqrSum
Definition: Controller.h:211
double BigReal
Definition: common.h:112
void Controller::multigratorTemperature ( int  step,
int  callNumber 
)
protected

Definition at line 853 of file Controller.C.

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

Referenced by integrate().

853  {
858  BigReal Nf = numDegFreedom;
860  BigReal T1, T2, v;
861  {
862  T1 = temperature;
863  BigReal kTinst = BOLTZMANN * temperature;
864  for (int i=n-1;i >= 0;i--) {
865  if (i == 0) {
866  BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
867  multigratorNuT[0] = exp(-0.5*t*NuOmega)*multigratorNu[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
868  } else if (i == n-1) {
869  multigratorNuT[n-1] = multigratorNu[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
870  } else {
871  BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
872  multigratorNuT[i] = exp(-0.5*t*NuOmega)*multigratorNu[i] +
873  0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
874  }
875  }
876  BigReal velScale = exp(-t*multigratorNuT[0]/multigratorOmega[0]);
877  v = velScale;
878  if (callNumber == 1)
879  broadcast->velocityRescaleFactor.publish(step,velScale);
880  else
881  broadcast->velocityRescaleFactor2.publish(step,velScale);
882  }
883 
884  {
885  // Wait here for Sequencer to finish scaling and re-calculating kinetic energy
888  temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
889  T2 = temperature;
890  if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
891  // If this is pressure cycle, receive new momentum product
892  GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
893  }
894  }
895 
896  // fprintf(stderr, "%d | T %lf scale %lf T' %lf\n", step, T1, v, T2);
897 
898  {
899  BigReal kTinst = BOLTZMANN * temperature;
900  for (int i=0;i < n;i++) {
901  if (i == 0) {
902  BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
903  multigratorNu[0] = exp(-0.5*t*NuOmega)*multigratorNuT[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
904  } else if (i == n-1) {
905  multigratorNu[n-1] = multigratorNuT[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
906  } else {
907  BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
908  multigratorNu[i] = exp(-0.5*t*NuOmega)*multigratorNuT[i] +
909  0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
910  }
912  }
913  }
914 
915  }
916 }
#define BOLTZMANN
Definition: common.h:45
std::vector< BigReal > multigratorOmega
Definition: Controller.h:215
std::vector< BigReal > multigratorNuT
Definition: Controller.h:214
#define GET_TENSOR(O, R, A)
Definition: ReductionMgr.h:59
BigReal temperature
Definition: Controller.h:161
std::vector< BigReal > multigratorNu
Definition: Controller.h:213
void require(void)
Definition: ReductionMgr.h:341
ControllerBroadcasts * broadcast
Definition: Controller.h:251
int64_t numDegFreedom
Definition: Controller.h:101
std::vector< BigReal > multigratorZeta
Definition: Controller.h:216
SimpleBroadcastObject< BigReal > velocityRescaleFactor2
Definition: Broadcasts.h:73
RequireReduction * multigratorReduction
Definition: Controller.h:217
int multigratorNoseHooverChainLength
BigReal multigratorTemperatureTarget
SimpleBroadcastObject< BigReal > velocityRescaleFactor
Definition: Broadcasts.h:68
void publish(int tag, const T &t)
BigReal multigratorTemperatureRelaxationTime
BigReal item(int i) const
Definition: ReductionMgr.h:340
int multigratorPressureFreq
SimParameters *const simParams
Definition: Controller.h:235
int multigratorTemperatureFreq
BigReal kineticEnergy
Definition: Controller.h:158
Tensor momentumSqrSum
Definition: Controller.h:211
double BigReal
Definition: common.h:112
void Controller::outputExtendedSystem ( int  step)
protected

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

3960 {
3961 
3962  if ( step >= 0 ) {
3963 
3964  // Write out eXtended System Trajectory (XST) file
3965  if ( simParams->xstFrequency &&
3966  ((step % simParams->xstFrequency) == 0) )
3967  {
3968  if ( ! xstFile.is_open() )
3969  {
3970  iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
3973  while (!xstFile) {
3974  if ( errno == EINTR ) {
3975  CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
3976  xstFile.clear();
3978  continue;
3979  }
3980  char err_msg[257];
3981  sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
3982  NAMD_err(err_msg);
3983  }
3984  xstFile << "# NAMD extended system trajectory file" << std::endl;
3986  }
3988  xstFile.flush();
3989  }
3990 
3991  // Write out eXtended System Configuration (XSC) files
3992  // Output a restart file
3993  if ( simParams->restartFrequency &&
3994  ((step % simParams->restartFrequency) == 0) &&
3995  (step != simParams->firstTimestep) )
3996  {
3997  iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
3998  << step << "\n" << endi;
3999  char fname[140];
4000  const char *bsuffix = ".old";
4001  strcpy(fname, simParams->restartFilename);
4002  if ( simParams->restartSave ) {
4003  char timestepstr[20];
4004  sprintf(timestepstr,".%d",step);
4005  strcat(fname, timestepstr);
4006  bsuffix = ".BAK";
4007  }
4008  strcat(fname, ".xsc");
4009  NAMD_backup_file(fname,bsuffix);
4010  ofstream_namd xscFile(fname);
4011  while (!xscFile) {
4012  if ( errno == EINTR ) {
4013  CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
4014  xscFile.clear();
4015  xscFile.open(fname);
4016  continue;
4017  }
4018  char err_msg[257];
4019  sprintf(err_msg, "Error opening XSC restart file %s",fname);
4020  NAMD_err(err_msg);
4021  }
4022  xscFile << "# NAMD extended system configuration restart file" << std::endl;
4023  writeExtendedSystemLabels(xscFile);
4024  writeExtendedSystemData(step,xscFile);
4025  if (!xscFile) {
4026  char err_msg[257];
4027  sprintf(err_msg, "Error writing XSC restart file %s",fname);
4028  NAMD_err(err_msg);
4029  }
4030  }
4031 
4032  }
4033 
4034  // Output final coordinates
4035  if (step == FILE_OUTPUT || step == END_OF_RUN)
4036  {
4037  int realstep = ( step == FILE_OUTPUT ?
4039  iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
4040  << realstep << "\n" << endi;
4041  static char fname[140];
4042  strcpy(fname, simParams->outputFilename);
4043  strcat(fname, ".xsc");
4044  NAMD_backup_file(fname);
4045  ofstream_namd xscFile(fname);
4046  while (!xscFile) {
4047  if ( errno == EINTR ) {
4048  CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
4049  xscFile.clear();
4050  xscFile.open(fname);
4051  continue;
4052  }
4053  char err_msg[257];
4054  sprintf(err_msg, "Error opening XSC output file %s",fname);
4055  NAMD_err(err_msg);
4056  }
4057  xscFile << "# NAMD extended system configuration output file" << std::endl;
4058  writeExtendedSystemLabels(xscFile);
4059  writeExtendedSystemData(realstep,xscFile);
4060  if (!xscFile) {
4061  char err_msg[257];
4062  sprintf(err_msg, "Error writing XSC output file %s",fname);
4063  NAMD_err(err_msg);
4064  }
4065  }
4066 
4067  // Close trajectory file
4068  if (step == END_OF_RUN) {
4069  if ( xstFile.is_open() ) {
4070  xstFile.close();
4071  iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
4072  }
4073  }
4074 
4075 }
char restartFilename[128]
void NAMD_err(const char *err_msg)
Definition: common.C:102
#define FILE_OUTPUT
Definition: Output.h:25
void writeExtendedSystemLabels(ofstream_namd &file)
Definition: Controller.C:3621
ofstream_namd & flush()
Definition: fstream_namd.C:17
#define iout
Definition: InfoStream.h:87
ofstream_namd xstFile
Definition: Controller.h:252
char outputFilename[128]
void writeExtendedSystemData(int step, ofstream_namd &file)
Definition: Controller.C:3634
bool is_open() const
Definition: fstream_namd.h:30
char xstFilename[128]
#define END_OF_RUN
Definition: Output.h:26
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:159
SimParameters *const simParams
Definition: Controller.h:235
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void open(const char *_fname, std::ios_base::openmode _mode=std::ios_base::out)
Definition: fstream_namd.C:11
void Controller::outputFepEnergy ( int  step)
protected

Definition at line 3667 of file Controller.C.

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

Referenced by integrate().

3667  {
3668  if (simParams->alchOn && simParams->alchFepOn) {
3669  const int stepInRun = step - simParams->firstTimestep;
3670  const int alchEquilSteps = simParams->alchEquilSteps;
3671  const BigReal alchLambda = simParams->alchLambda;
3672  const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
3673 
3674  if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
3675  FepNo = 0;
3676  exp_dE_ByRT = 0.0;
3677  net_dE = 0.0;
3678  }
3682 
3683  if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
3684  (step / simParams->alchIDWSFreq) % 2 == 1 )) {
3685  // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
3686  FepNo++;
3687  exp_dE_ByRT += exp(-dE/RT);
3688  net_dE += dE;
3689  }
3690 
3691  if (simParams->alchOutFreq) {
3692  if (stepInRun == 0) {
3693  if (!fepFile.is_open()) {
3696  iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
3697  if(alchEnsembleAvg){
3698  fepSum = 0.0;
3699  fepFile << "# STEP Elec "
3700  << "vdW dE dE_avg Temp dG\n"
3701  << "# l l+dl "
3702  << " l l+dl E(l+dl)-E(l)" << std::endl;
3703  }
3704  else{
3705  fepFile << "# STEP Elec "
3706  << "vdW dE Temp\n"
3707  << "# l l+dl "
3708  << " l l+dl E(l+dl)-E(l)" << std::endl;
3709  }
3710  }
3711  if(!step){
3712  fepFile << "#NEW FEP WINDOW: "
3713  << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
3714  << simParams->alchLambda2;
3715  if ( simParams->alchLambdaIDWS >= 0. ) {
3716  fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
3717  }
3718  fepFile << std::endl;
3719  }
3720  }
3721  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3722  fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3723  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3724  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3725  }
3726  if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
3727  writeFepEnergyData(step, fepFile);
3728  fepFile.flush();
3729  }
3730  if (alchEnsembleAvg && (step == simParams->N)) {
3731  fepSum = fepSum + dG;
3732  fepFile << "#Free energy change for lambda window [ " << alchLambda
3733  << " " << simParams->alchLambda2 << " ] is " << dG
3734  << " ; net change until now is " << fepSum << std::endl;
3735  fepFile.flush();
3736  }
3737  }
3738  }
3739 }
BigReal net_dE
Definition: Controller.h:120
BigReal dG
Definition: Controller.h:121
BigReal fepSum
Definition: Controller.h:124
#define BOLTZMANN
Definition: common.h:45
BigReal ljEnergy
Definition: Controller.h:106
BigReal alchLambda2
BigReal electEnergySlow_f
Definition: Controller.h:115
ofstream_namd & flush()
Definition: fstream_namd.C:17
#define iout
Definition: InfoStream.h:87
BigReal electEnergySlow
Definition: Controller.h:105
BigReal alchLambda
BigReal electEnergy_f
Definition: Controller.h:114
BigReal alchLambdaIDWS
BigReal bondedEnergyDiff_f
Definition: Controller.h:113
bool is_open() const
Definition: fstream_namd.h:30
BigReal dE
Definition: Controller.h:119
BigReal exp_dE_ByRT
Definition: Controller.h:118
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:159
BigReal alchTemp
char alchOutFile[128]
BigReal electEnergy
Definition: Controller.h:104
BigReal ljEnergy_f
Definition: Controller.h:116
SimParameters *const simParams
Definition: Controller.h:235
infostream & endi(infostream &s)
Definition: InfoStream.C:38
ofstream_namd fepFile
Definition: Controller.h:258
void writeFepEnergyData(int step, ofstream_namd &file)
Definition: Controller.C:3899
void open(const char *_fname, std::ios_base::openmode _mode=std::ios_base::out)
Definition: fstream_namd.C:11
double BigReal
Definition: common.h:112
void Controller::outputTiEnergy ( int  step)
protected

Definition at line 3741 of file Controller.C.

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

Referenced by integrate().

3741  {
3743  const int stepInRun = step - simParams->firstTimestep;
3744  const int alchEquilSteps = simParams->alchEquilSteps;
3745  const int alchLambdaFreq = simParams->alchLambdaFreq;
3746 
3747  if (stepInRun == 0 || stepInRun == alchEquilSteps) {
3748  TiNo = 0;
3749  net_dEdl_bond_1 = 0;
3750  net_dEdl_bond_2 = 0;
3751  net_dEdl_elec_1 = 0;
3752  net_dEdl_elec_2 = 0;
3753  net_dEdl_lj_1 = 0;
3754  net_dEdl_lj_2 = 0;
3755  }
3756  TiNo++;
3765 
3766  // Don't accumulate block averages (you'll get a divide by zero!) or even
3767  // open the TI output if alchOutFreq is zero.
3768  if (simParams->alchOutFreq) {
3769  if (stepInRun == 0 || stepInRun == alchEquilSteps
3770  || (! ((step - 1) % simParams->alchOutFreq))) {
3771  // output of instantaneous dU/dl now replaced with running average
3772  // over last alchOutFreq steps (except for step 0)
3773  recent_TiNo = 0;
3774  recent_dEdl_bond_1 = 0;
3775  recent_dEdl_bond_2 = 0;
3776  recent_dEdl_elec_1 = 0;
3777  recent_dEdl_elec_2 = 0;
3778  recent_dEdl_lj_1 = 0;
3779  recent_dEdl_lj_2 = 0;
3780  recent_alchWork = 0;
3781  }
3782  recent_TiNo++;
3786  + electEnergyPME_ti_1);
3788  + electEnergyPME_ti_2);
3792 
3793  if (stepInRun == 0) {
3794  if (!tiFile.is_open()) {
3797  /* BKR - This has been rather drastically updated to better match
3798  stdout. This was necessary for several reasons:
3799  1) PME global scaling is obsolete (now removed)
3800  2) scaling of bonded terms was added
3801  3) alchemical work is now accumulated when switching is active
3802  */
3803  iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
3804  tiFile << "#TITITLE: TS";
3805  tiFile << FORMAT("BOND1");
3806  tiFile << FORMAT("AVGBOND1");
3807  tiFile << FORMAT("ELECT1");
3808  tiFile << FORMAT("AVGELECT1");
3809  tiFile << " ";
3810  tiFile << FORMAT("VDW1");
3811  tiFile << FORMAT("AVGVDW1");
3812  tiFile << FORMAT("BOND2");
3813  tiFile << FORMAT("AVGBOND2");
3814  tiFile << FORMAT("ELECT2");
3815  tiFile << " ";
3816  tiFile << FORMAT("AVGELECT2");
3817  tiFile << FORMAT("VDW2");
3818  tiFile << FORMAT("AVGVDW2");
3819  if (alchLambdaFreq > 0) {
3820  tiFile << FORMAT("ALCHWORK");
3821  tiFile << FORMAT("CUMALCHWORK");
3822  }
3823  tiFile << std::endl;
3824  }
3825 
3826  if (alchLambdaFreq > 0) {
3827  tiFile << "#ALCHEMICAL SWITCHING ACTIVE "
3828  << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
3829  << "\n#LAMBDA SCHEDULE: "
3830  << "dL: " << simParams->getLambdaDelta()
3831  << " Freq: " << alchLambdaFreq;
3832  }
3833  else {
3834  const BigReal alchLambda = simParams->alchLambda;
3835  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
3836  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
3837  const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
3838  const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
3839  const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
3840  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
3841  tiFile << "#NEW TI WINDOW: "
3842  << "LAMBDA " << alchLambda
3843  << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
3844  << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
3845  << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
3846  << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
3847  }
3848  tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
3849  << std::endl;
3850  }
3851 
3852  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3853  tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3854  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3855  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3856  }
3857  if ((step%simParams->alchOutFreq) == 0) {
3858  writeTiEnergyData(step, tiFile);
3859  tiFile.flush();
3860  }
3861  }
3862  }
3863 }
ofstream_namd tiFile
Definition: Controller.h:262
BigReal electEnergy_ti_1
Definition: Controller.h:128
BigReal net_dEdl_lj_2
Definition: Controller.h:139
BigReal getLambdaDelta(void)
BigReal net_dEdl_lj_1
Definition: Controller.h:138
BigReal net_dEdl_bond_2
Definition: Controller.h:135
ofstream_namd & flush()
Definition: fstream_namd.C:17
BigReal recent_alchWork
Definition: Controller.h:150
#define iout
Definition: InfoStream.h:87
BigReal recent_dEdl_bond_2
Definition: Controller.h:145
BigReal alchLambda
BigReal alchWork
Definition: Controller.h:151
BigReal ljEnergy_ti_1
Definition: Controller.h:130
BigReal recent_dEdl_elec_2
Definition: Controller.h:147
BigReal recent_dEdl_elec_1
Definition: Controller.h:146
BigReal getBondLambda(const BigReal)
int recent_TiNo
Definition: Controller.h:152
bool is_open() const
Definition: fstream_namd.h:30
BigReal getCurrentLambda2(const int)
BigReal electEnergySlow_ti_1
Definition: Controller.h:129
BigReal electEnergyPME_ti_1
Definition: Controller.h:141
static char * FORMAT(BigReal X)
Definition: ComputeQM.C:367
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:159
BigReal electEnergy_ti_2
Definition: Controller.h:131
BigReal net_dEdl_elec_1
Definition: Controller.h:136
BigReal alchTemp
char alchOutFile[128]
BigReal recent_dEdl_lj_1
Definition: Controller.h:148
BigReal electEnergyPME_ti_2
Definition: Controller.h:142
BigReal electEnergySlow_ti_2
Definition: Controller.h:132
BigReal getVdwLambda(const BigReal)
SimParameters *const simParams
Definition: Controller.h:235
BigReal net_dEdl_elec_2
Definition: Controller.h:137
infostream & endi(infostream &s)
Definition: InfoStream.C:38
BigReal recent_dEdl_lj_2
Definition: Controller.h:149
BigReal bondedEnergy_ti_2
Definition: Controller.h:127
void writeTiEnergyData(int step, ofstream_namd &file)
Definition: Controller.C:3935
void open(const char *_fname, std::ios_base::openmode _mode=std::ios_base::out)
Definition: fstream_namd.C:11
BigReal net_dEdl_bond_1
Definition: Controller.h:134
BigReal bondedEnergy_ti_1
Definition: Controller.h:126
BigReal getElecLambda(const BigReal)
double BigReal
Definition: common.h:112
BigReal ljEnergy_ti_2
Definition: Controller.h:133
BigReal recent_dEdl_bond_1
Definition: Controller.h:144
void Controller::printDynamicsEnergies ( int  step)
protected

Definition at line 3015 of file Controller.C.

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

Referenced by integrate().

3015  {
3016 
3017  Node *node = Node::Object();
3018  Molecule *molecule = node->molecule;
3019  SimParameters *simParameters = node->simParameters;
3020  Lattice &lattice = state->lattice;
3021 
3022  compareChecksums(step);
3023 
3024  printEnergies(step,0);
3025 }
static Node * Object()
Definition: Node.h:86
void compareChecksums(int, int=0)
Definition: Controller.C:2767
Definition: Node.h:78
SimParameters * simParameters
Definition: Node.h:178
NamdState *const state
Definition: Controller.h:236
void printEnergies(int step, int minimize)
Definition: Controller.C:3027
Molecule * molecule
Definition: Node.h:176
void Controller::printEnergies ( int  step,
int  minimize 
)
protected

Definition at line 3027 of file Controller.C.

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

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

3028 {
3029  Node *node = Node::Object();
3030  Molecule *molecule = node->molecule;
3031  SimParameters *simParameters = node->simParameters;
3032  Lattice &lattice = state->lattice;
3033 
3034  // Drude model ANISO energy is added into BOND energy
3035  // and THOLE energy is added into ELECT energy
3036 
3037  BigReal bondEnergy;
3038  BigReal angleEnergy;
3039  BigReal dihedralEnergy;
3040  BigReal improperEnergy;
3041  BigReal crosstermEnergy;
3042  BigReal boundaryEnergy;
3043  BigReal miscEnergy;
3044  BigReal potentialEnergy;
3045  BigReal flatEnergy;
3046  BigReal smoothEnergy;
3047  BigReal work;
3048 
3049  Vector momentum;
3050  Vector angularMomentum;
3051  BigReal volume = lattice.volume();
3052 
3053  bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
3054  angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
3055  dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
3056  improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
3057  crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
3058  boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
3059  miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
3060 
3061  if ( minimize || ! ( step % nbondFreq ) )
3062  {
3065 
3066  // JLai
3069 
3073 
3074 //fepb
3078 
3085 //fepe
3086  }
3087 
3088  if ( minimize || ! ( step % slowFreq ) )
3089  {
3091 //fepb
3093 
3098 //fepe
3099  }
3100 
3101  if (simParameters->LJcorrection && volume) {
3102 #ifdef MEM_OPT_VERSION
3103  NAMD_bug("LJcorrection not supported in memory optimized build.");
3104 #else
3105  // Apply tail correction to energy.
3106  BigReal alchLambda = simParameters->getCurrentLambda(step);
3107  ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
3108 
3109  if (simParameters->alchOn) {
3110  if (simParameters->alchFepOn) {
3111  BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
3112  ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
3113  } else if (simParameters->alchThermIntOn) {
3114  ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
3115  ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
3116  }
3117  }
3118 #endif
3119  }
3120 
3121 //fepb BKR - Compute alchemical work if using dynamic lambda. This is here so
3122 // that the cumulative work can be given during a callback.
3123  if (simParameters->alchLambdaFreq > 0) {
3124  if (step <=
3125  simParameters->firstTimestep + simParameters->alchEquilSteps) {
3126  cumAlchWork = 0.0;
3127  }
3128  alchWork = computeAlchWork(step);
3129  cumAlchWork += alchWork;
3130  }
3131 //fepe
3132 
3133  momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
3134  momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
3135  momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
3136  angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
3137  angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
3138  angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
3139 
3140  // Ported by JLai
3141  potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
3142  + improperEnergy + electEnergy + electEnergySlow + ljEnergy
3143  + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy
3145  // End of port
3146  totalEnergy = potentialEnergy + kineticEnergy;
3147  flatEnergy = (totalEnergy
3149  if ( !(step%slowFreq) ) {
3150  // only adjust based on most accurate energies
3152  if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
3153  if ( step != simParams->firstTimestep ) {
3154  smooth2_avg *= 0.9375;
3155  smooth2_avg += 0.0625 * s;
3156  }
3157  }
3158  smoothEnergy = (flatEnergy + smooth2_avg
3160 
3161  // Reset values for accumulated heat and work.
3162  if (step <= simParameters->firstTimestep &&
3163  (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
3164  heat = 0.0;
3166  }
3167  if ( simParameters->outputMomenta && ! minimize &&
3168  ! ( step % simParameters->outputMomenta ) )
3169  {
3170  iout << "MOMENTUM: " << step
3171  << " P: " << momentum
3172  << " L: " << angularMomentum
3173  << "\n" << endi;
3174  }
3175 
3176  if ( simParameters->outputPressure ) {
3179  tavg_count += 1;
3180  if ( minimize || ! ( step % simParameters->outputPressure ) ) {
3181  iout << "PRESSURE: " << step << " "
3182  << PRESSUREFACTOR * pressure << "\n"
3183  << "GPRESSURE: " << step << " "
3184  << PRESSUREFACTOR * groupPressure << "\n";
3185  if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
3186  << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
3187  << "GPRESSAVG: " << step << " "
3189  iout << endi;
3190  pressure_tavg = 0;
3191  groupPressure_tavg = 0;
3192  tavg_count = 0;
3193  }
3194  }
3195 
3196  // pressure profile reductions
3197  if (pressureProfileSlabs) {
3198  const int freq = simParams->pressureProfileFreq;
3199  const int arraysize = 3*pressureProfileSlabs;
3200 
3201  BigReal *total = new BigReal[arraysize];
3202  memset(total, 0, arraysize*sizeof(BigReal));
3203  const int first = simParams->firstTimestep;
3204 
3205  if (ppbonded) ppbonded->getData(first, step, lattice, total);
3206  if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
3207  if (ppint) ppint->getData(first, step, lattice, total);
3208  for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
3210 
3211  if (!(step % freq)) {
3212  // convert NAMD internal virial to pressure in units of bar
3214 
3215  iout << "PRESSUREPROFILE: " << step << " ";
3216  if (step == first) {
3217  // output pressure profile for this step
3218  for (int i=0; i<arraysize; i++) {
3219  iout << total[i]