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:47
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
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:524
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
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ cudaTextureObject_t const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t 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__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ const int numPatches
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:114
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 2351 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().

2351  {
2352  if (!simParams->adaptTempOn) return;
2353  iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
2354  adaptTempDtMin = 0;
2355  adaptTempDtMax = 0;
2356  adaptTempAutoDt = false;
2357  if (simParams->adaptTempBins == 0) {
2358  iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
2359  std::ifstream adaptTempRead(simParams->adaptTempInFile);
2360  if (adaptTempRead) {
2361  int readInt;
2362  BigReal readReal;
2363  bool readBool;
2364  // step
2365  adaptTempRead >> readInt;
2366  // Start with min and max temperatures
2367  adaptTempRead >> adaptTempT; // KELVIN
2368  adaptTempRead >> adaptTempBetaMin; // KELVIN
2369  adaptTempRead >> adaptTempBetaMax; // KELVIN
2370  adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
2371  adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
2372  // In case file is manually edited
2373  if (adaptTempBetaMin > adaptTempBetaMax){
2374  readReal = adaptTempBetaMax;
2375  adaptTempBetaMax = adaptTempBetaMin;
2376  adaptTempBetaMin = adaptTempBetaMax;
2377  }
2378  adaptTempRead >> adaptTempBins;
2379  adaptTempRead >> adaptTempCg;
2380  adaptTempRead >> adaptTempDt;
2388  adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
2389  for(int j = 0; j < adaptTempBins; ++j) {
2390  adaptTempRead >> adaptTempPotEnergyAve[j];
2391  adaptTempRead >> adaptTempPotEnergyVar[j];
2392  adaptTempRead >> adaptTempPotEnergySamples[j];
2393  adaptTempRead >> adaptTempPotEnergyAveNum[j];
2394  adaptTempRead >> adaptTempPotEnergyVarNum[j];
2395  adaptTempRead >> adaptTempPotEnergyAveDen[j];
2396  adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
2397  }
2398  adaptTempRead.close();
2399  }
2400  else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
2401  }
2402  else {
2403  adaptTempBins = simParams->adaptTempBins;
2411  adaptTempBetaMax = 1./simParams->adaptTempTmin;
2412  adaptTempBetaMin = 1./simParams->adaptTempTmax;
2413  adaptTempCg = simParams->adaptTempCgamma;
2414  adaptTempDt = simParams->adaptTempDt;
2415  adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
2416  adaptTempT = simParams->initialTemp;
2417  if (simParams->langevinOn)
2418  adaptTempT = simParams->langevinTemp;
2419  else if (simParams->rescaleFreq > 0)
2420  adaptTempT = simParams->rescaleTemp;
2421  for(int j = 0; j < adaptTempBins; ++j){
2422  adaptTempPotEnergyAveNum[j] = 0.;
2423  adaptTempPotEnergyAveDen[j] = 0.;
2425  adaptTempPotEnergyVarNum[j] = 0.;
2426  adaptTempPotEnergyVar[j] = 0.;
2427  adaptTempPotEnergyAve[j] = 0.;
2428  adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
2429  }
2430  }
2431  if (simParams->adaptTempAutoDt > 0.0) {
2432  adaptTempAutoDt = true;
2435  }
2436  adaptTempDTave = 0;
2437  adaptTempDTavenum = 0;
2438  iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
2439  iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
2440  iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
2441  iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
2442  iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
2443  if (simParams->adaptTempRestartFreq > 0) {
2447  NAMD_die("Error opening restart file");
2448  }
2449 }
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:81
BigReal adaptTempCgamma
BigReal adaptTempTmin
BigReal adaptTempT
Definition: Controller.h:314
BigReal adaptTempDt
Definition: Controller.h:323
char adaptTempInFile[NAMD_FILENAME_BUFFER_SIZE]
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
BigReal * adaptTempPotEnergyAveNum
Definition: Controller.h:307
#define iout
Definition: InfoStream.h:51
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
BigReal adaptTempDtMax
Definition: Controller.h:326
char adaptTempRestartFile[NAMD_FILENAME_BUFFER_SIZE]
BigReal langevinTemp
Bool adaptTempAutoDt
Definition: Controller.h:324
void NAMD_die(const char *err_msg)
Definition: common.C:85
BigReal adaptTempDt
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:167
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
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:114
void Controller::adaptTempUpdate ( int  step,
int  minimize = 0 
)
protected

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

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

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

2451  {
2453  adaptTempRestartFile.seekp(std::ios::beg);
2454  iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
2455  adaptTempRestartFile << step << " ";
2456  // Start with min and max temperatures
2457  adaptTempRestartFile << adaptTempT<< " "; // KELVIN
2458  adaptTempRestartFile << 1./adaptTempBetaMin << " "; // KELVIN
2459  adaptTempRestartFile << 1./adaptTempBetaMax << " "; // KELVIN
2463  adaptTempRestartFile << "\n" ;
2464  for(int j = 0; j < adaptTempBins; ++j) {
2471  adaptTempRestartFile << "\n";
2472  }
2474  }
2475 }
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
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
ofstream_namd & flush()
Definition: fstream_namd.C:17
BigReal * adaptTempPotEnergyAveNum
Definition: Controller.h:307
#define iout
Definition: InfoStream.h:51
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
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:3664
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
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
if(ComputeNonbondedUtil::goMethod==2)
#define iout
Definition: InfoStream.h:51
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:129
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:85
void terminate(void)
Definition: Controller.C:4178
BigReal initialTemp
SimParameters *const simParams
Definition: Controller.h:235
void sendCheckpointReq(int remote, const char *key, int task, Lattice &lat, ControllerState &cs)
Definition: Node.C:1280
void outputExtendedSystem(int step)
Definition: Controller.C:3968
#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
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
int berendsenPressureFreq
#define iout
Definition: InfoStream.h:51
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:83
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 1736 of file Controller.C.

Referenced by rescaleaccelMD().

1737  {
1738  switch(iE){
1739  case 2:
1740  *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
1741  // if k0 <=0 OR k0 > 1, switch to iE=1 mode
1742  if(*k0 > 1.)
1743  *warn |= 1;
1744  else if(*k0 <= 0.)
1745  *warn |= 2;
1746  //else stay in iE=2 mode
1747  else{
1748  *E = Vmin + (Vmax-Vmin)/(*k0);
1749  *iEused = 2;
1750  break;
1751  }
1752  case 1:
1753  *E = Vmax;
1754  *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
1755  if(*k0 > 1.)
1756  *k0 = 1.;
1757  *iEused = 1;
1758  break;
1759  }
1760 
1761  *k = *k0 / (Vmax - Vmin);
1762 }
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 1765 of file Controller.C.

Referenced by rescaleaccelMD().

1766  {
1767  BigReal VE = testV - E;
1768  //if V < E, apply boost
1769  if(VE < 0){
1770  *factor = k * VE;
1771  *vir += vir_orig * (*factor);
1772  *dV += (*factor) * VE/2.;
1773  *factor += 1.;
1774  }
1775  else{
1776  *factor = 1.;
1777  }
1778 }
double BigReal
Definition: common.h:114
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 1716 of file Controller.C.

Referenced by rescaleaccelMD().

1717  {
1718  BigReal delta;
1719 
1720  if(testV > *Vmax){
1721  *Vmax = testV;
1722  }
1723  else if(testV < *Vmin){
1724  *Vmin = testV;
1725  }
1726 
1727  //mean and std calculated by Online Algorithm
1728  delta = testV - *Vavg;
1729  *Vavg += delta / (BigReal)n;
1730  *M2 += delta * (testV - (*Vavg));
1731 
1732  *sigmaV = sqrt(*M2/n);
1733 }
double BigReal
Definition: common.h:114
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 1598 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().

1601  {
1602 
1603  Tensor virial_normal = virial_normal_in;
1604  Tensor virial_nbond = virial_nbond_in;
1605  Tensor virial_slow = virial_slow_in;
1606 
1607  // Tensor tmp = virial_normal;
1608  // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
1609  // tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
1610 
1611  Node *node = Node::Object();
1612  Molecule *molecule = node->molecule;
1613  SimParameters *simParameters = node->simParameters;
1614  Lattice &lattice = state->lattice;
1615 
1616  BigReal volume;
1617 
1618  Vector extPosition = lattice.origin();
1619  virial_normal -= outer(extForce_normal,extPosition);
1620  virial_nbond -= outer(extForce_nbond,extPosition);
1621  virial_slow -= outer(extForce_slow,extPosition);
1622 
1623  if ( (volume=lattice.volume()) != 0. )
1624  {
1625 
1626  if (simParameters->LJcorrection && volume) {
1627 #ifdef MEM_OPT_VERSION
1628  NAMD_bug("LJcorrection not supported in memory optimized build.");
1629 #else
1630  // Apply tail correction to pressure
1631  BigReal alchLambda = simParameters->getCurrentLambda(step);
1632  virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
1633 #endif
1634  }
1635 
1636  // kinetic energy component included in virials
1637  pressure_normal = virial_normal / volume;
1638  groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
1639 
1640  if (simParameters->accelMDOn) {
1641  pressure_amd = virial_amd / volume;
1644  }
1645 
1646  if ( minimize || ! ( step % nbondFreq ) )
1647  {
1648  pressure_nbond = virial_nbond / volume;
1649  groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
1650  }
1651 
1652  if ( minimize || ! ( step % slowFreq ) )
1653  {
1654  pressure_slow = virial_slow / volume;
1655  groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
1656  }
1657 
1661  }
1662  else
1663  {
1664  pressure = Tensor();
1665  groupPressure = Tensor();
1666  }
1667 
1668  if ( simParameters->useGroupPressure )
1669  {
1674  }
1675  else
1676  {
1681  }
1682 
1683  if ( simParameters->useFlexibleCell ) {
1684  // use symmetric pressure to control rotation
1685  // controlPressure_normal = symmetric(controlPressure_normal);
1686  // controlPressure_nbond = symmetric(controlPressure_nbond);
1687  // controlPressure_slow = symmetric(controlPressure_slow);
1688  // controlPressure = symmetric(controlPressure);
1689  // only use on-diagonal components for now
1694  if ( simParameters->useConstantRatio ) {
1695 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
1696  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
1701 #undef AVGXY
1702  }
1703  } else {
1710  controlPressure =
1711  Tensor::identity(trace(controlPressure)/3.);
1712  }
1713 }
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:129
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:114
void Controller::compareChecksums ( int  step,
int  forgiving = 0 
)
protected

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

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

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

3879  {
3880  // alchemical scaling factors for groups 1/2 at the previous lambda
3881  const BigReal oldLambda = simParams->getCurrentLambda(step-1);
3882  const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
3883  const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
3884  const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
3885  const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
3886  const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
3887  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
3888  // alchemical scaling factors for groups 1/2 at the new lambda
3889  const BigReal alchLambda = simParams->getCurrentLambda(step);
3890  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
3891  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
3892  const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
3893  const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
3894  const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
3895  const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda);
3896 
3897  return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
3898  (elec_lambda_12 - elec_lambda_1)*
3900  (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
3901  (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
3902  (elec_lambda_22 - elec_lambda_2)*
3904  (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
3905  );
3906 }
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:114
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
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
BigReal z
Definition: Vector.h:66
#define iout
Definition: InfoStream.h:51
ControllerBroadcasts * broadcast
Definition: Controller.h:251
BigReal x
Definition: Vector.h:66
void NAMD_die(const char *err_msg)
Definition: common.C:85
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:83
RequireReduction * reduction
Definition: Controller.h:237
double BigReal
Definition: common.h:114
void Controller::cycleBarrier ( int  doBarrier,
int  step 
)
protected

Definition at line 4140 of file Controller.C.

References broadcast.

Referenced by integrate().

4140  {
4141 #if USE_BARRIER
4142  if (doBarrier) {
4143  broadcast->cycleBarrier.publish(step,1);
4144  CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
4145  CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4146  }
4147 #endif
4148 }
ControllerBroadcasts * broadcast
Definition: Controller.h:251
void Controller::enqueueCollections ( int  timestep)
protected

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

3665 {
3666  if ( Output::coordinateNeeded(timestep) ){
3667  collection->enqueuePositions(timestep,state->lattice);
3668  }
3669  if ( Output::velocityNeeded(timestep) )
3670  collection->enqueueVelocities(timestep);
3671  if ( Output::forceNeeded(timestep) )
3672  collection->enqueueForces(timestep);
3673 }
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:3664
void cycleBarrier(int, int)
Definition: Controller.C:4140
void rescaleVelocities(int)
Definition: Controller.C:1233
void rescaleaccelMD(int step, int minimize=0)
Definition: Controller.C:1839
int nbondFreq
Definition: Controller.h:79
int eventEndOfTimeStep
Definition: Node.C:286
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:3750
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:1420
void reassignVelocities(int)
Definition: Controller.C:1308
void outputFepEnergy(int step)
Definition: Controller.C:3676
void printFepMessage(int)
Definition: Controller.C:1274
int slowFreq
Definition: Controller.h:80
void traceBarrier(int, int)
Definition: Controller.C:4150
void printDynamicsEnergies(int)
Definition: Controller.C:3023
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:4126
void outputExtendedSystem(int step)
Definition: Controller.C:3968
void adaptTempUpdate(int step, int minimize=0)
Definition: Controller.C:2477
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:81
int nbondFreq
Definition: Controller.h:79
BigReal langevinPistonTemp
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
#define BOLTZMANN
Definition: common.h:47
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:51
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:114
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:81
int nbondFreq
Definition: Controller.h:79
BigReal langevinPistonTemp
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
#define BOLTZMANN
Definition: common.h:47
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:51
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:114
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:3664
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
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
BigReal minLineGoal
if(ComputeNonbondedUtil::goMethod==2)
BigReal u
Definition: Controller.C:590
#define iout
Definition: InfoStream.h:51
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
gridSize x
double BigReal
Definition: common.h:114
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:47
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:114
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:1598
#define BOLTZMANN
Definition: common.h:47
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:114
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:47
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:114
void Controller::outputExtendedSystem ( int  step)
protected

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

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

3676  {
3677  if (simParams->alchOn && simParams->alchFepOn) {
3678  const int stepInRun = step - simParams->firstTimestep;
3679  const int alchEquilSteps = simParams->alchEquilSteps;
3680  const BigReal alchLambda = simParams->alchLambda;
3681  const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
3682 
3683  if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
3684  FepNo = 0;
3685  exp_dE_ByRT = 0.0;
3686  net_dE = 0.0;
3687  }
3691 
3692  if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
3693  (step / simParams->alchIDWSFreq) % 2 == 1 )) {
3694  // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
3695  FepNo++;
3696  exp_dE_ByRT += exp(-dE/RT);
3697  net_dE += dE;
3698  }
3699 
3700  if (simParams->alchOutFreq) {
3701  if (stepInRun == 0) {
3702  if (!fepFile.is_open()) {
3705  iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
3706  if(alchEnsembleAvg){
3707  fepSum = 0.0;
3708  fepFile << "# STEP Elec "
3709  << "vdW dE dE_avg Temp dG\n"
3710  << "# l l+dl "
3711  << " l l+dl E(l+dl)-E(l)" << std::endl;
3712  }
3713  else{
3714  fepFile << "# STEP Elec "
3715  << "vdW dE Temp\n"
3716  << "# l l+dl "
3717  << " l l+dl E(l+dl)-E(l)" << std::endl;
3718  }
3719  }
3720  if(!step){
3721  fepFile << "#NEW FEP WINDOW: "
3722  << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
3723  << simParams->alchLambda2;
3724  if ( simParams->alchLambdaIDWS >= 0. ) {
3725  fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
3726  }
3727  fepFile << std::endl;
3728  }
3729  }
3730  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3731  fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3732  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3733  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3734  }
3735  if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
3736  writeFepEnergyData(step, fepFile);
3737  fepFile.flush();
3738  }
3739  if (alchEnsembleAvg && (step == simParams->N)) {
3740  fepSum = fepSum + dG;
3741  fepFile << "#Free energy change for lambda window [ " << alchLambda
3742  << " " << simParams->alchLambda2 << " ] is " << dG
3743  << " ; net change until now is " << fepSum << std::endl;
3744  fepFile.flush();
3745  }
3746  }
3747  }
3748 }
BigReal net_dE
Definition: Controller.h:120
BigReal dG
Definition: Controller.h:121
BigReal fepSum
Definition: Controller.h:124
#define BOLTZMANN
Definition: common.h:47
BigReal ljEnergy
Definition: Controller.h:106
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
BigReal alchLambda2
BigReal electEnergySlow_f
Definition: Controller.h:115
ofstream_namd & flush()
Definition: fstream_namd.C:17
#define iout
Definition: InfoStream.h:51
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
char alchOutFile[NAMD_FILENAME_BUFFER_SIZE]
BigReal exp_dE_ByRT
Definition: Controller.h:118
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:167
BigReal alchTemp
BigReal electEnergy
Definition: Controller.h:104
BigReal ljEnergy_f
Definition: Controller.h:116
SimParameters *const simParams
Definition: Controller.h:235
ofstream_namd fepFile
Definition: Controller.h:258
void writeFepEnergyData(int step, ofstream_namd &file)
Definition: Controller.C:3908
void open(const char *_fname, std::ios_base::openmode _mode=std::ios_base::out)
Definition: fstream_namd.C:11
double BigReal
Definition: common.h:114
void Controller::outputTiEnergy ( int  step)
protected

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

3750  {
3752  const int stepInRun = step - simParams->firstTimestep;
3753  const int alchEquilSteps = simParams->alchEquilSteps;
3754  const int alchLambdaFreq = simParams->alchLambdaFreq;
3755 
3756  if (stepInRun == 0 || stepInRun == alchEquilSteps) {
3757  TiNo = 0;
3758  net_dEdl_bond_1 = 0;
3759  net_dEdl_bond_2 = 0;
3760  net_dEdl_elec_1 = 0;
3761  net_dEdl_elec_2 = 0;
3762  net_dEdl_lj_1 = 0;
3763  net_dEdl_lj_2 = 0;
3764  }
3765  TiNo++;
3774 
3775  // Don't accumulate block averages (you'll get a divide by zero!) or even
3776  // open the TI output if alchOutFreq is zero.
3777  if (simParams->alchOutFreq) {
3778  if (stepInRun == 0 || stepInRun == alchEquilSteps
3779  || (! ((step - 1) % simParams->alchOutFreq))) {
3780  // output of instantaneous dU/dl now replaced with running average
3781  // over last alchOutFreq steps (except for step 0)
3782  recent_TiNo = 0;
3783  recent_dEdl_bond_1 = 0;
3784  recent_dEdl_bond_2 = 0;
3785  recent_dEdl_elec_1 = 0;
3786  recent_dEdl_elec_2 = 0;
3787  recent_dEdl_lj_1 = 0;
3788  recent_dEdl_lj_2 = 0;
3789  recent_alchWork = 0;
3790  }
3791  recent_TiNo++;
3795  + electEnergyPME_ti_1);
3797  + electEnergyPME_ti_2);
3801 
3802  if (stepInRun == 0) {
3803  if (!tiFile.is_open()) {
3806  /* BKR - This has been rather drastically updated to better match
3807  stdout. This was necessary for several reasons:
3808  1) PME global scaling is obsolete (now removed)
3809  2) scaling of bonded terms was added
3810  3) alchemical work is now accumulated when switching is active
3811  */
3812  iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
3813  tiFile << "#TITITLE: TS";
3814  tiFile << FORMAT("BOND1");
3815  tiFile << FORMAT("AVGBOND1");
3816  tiFile << FORMAT("ELECT1");
3817  tiFile << FORMAT("AVGELECT1");
3818  tiFile << " ";
3819  tiFile << FORMAT("VDW1");
3820  tiFile << FORMAT("AVGVDW1");
3821  tiFile << FORMAT("BOND2");
3822  tiFile << FORMAT("AVGBOND2");
3823  tiFile << FORMAT("ELECT2");
3824  tiFile << " ";
3825  tiFile << FORMAT("AVGELECT2");
3826  tiFile << FORMAT("VDW2");
3827  tiFile << FORMAT("AVGVDW2");
3828  if (alchLambdaFreq > 0) {
3829  tiFile << FORMAT("ALCHWORK");
3830  tiFile << FORMAT("CUMALCHWORK");
3831  }
3832  tiFile << std::endl;
3833  }
3834 
3835  if (alchLambdaFreq > 0) {
3836  tiFile << "#ALCHEMICAL SWITCHING ACTIVE "
3837  << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
3838  << "\n#LAMBDA SCHEDULE: "
3839  << "dL: " << simParams->getLambdaDelta()
3840  << " Freq: " << alchLambdaFreq;
3841  }
3842  else {
3843  const BigReal alchLambda = simParams->alchLambda;
3844  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
3845  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
3846  const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
3847  const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
3848  const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
3849  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
3850  tiFile << "#NEW TI WINDOW: "
3851  << "LAMBDA " << alchLambda
3852  << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
3853  << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
3854  << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
3855  << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
3856  }
3857  tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
3858  << std::endl;
3859  }
3860 
3861  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3862  tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3863  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3864  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3865  }
3866  if ((step%simParams->alchOutFreq) == 0) {
3867  writeTiEnergyData(step, tiFile);
3868  tiFile.flush();
3869  }
3870  }
3871  }
3872 }
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
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
ofstream_namd & flush()
Definition: fstream_namd.C:17
BigReal recent_alchWork
Definition: Controller.h:150
#define iout
Definition: InfoStream.h:51
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
char alchOutFile[NAMD_FILENAME_BUFFER_SIZE]
static char * FORMAT(BigReal X)
Definition: ComputeQM.C:367
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:167
BigReal electEnergy_ti_2
Definition: Controller.h:131
BigReal net_dEdl_elec_1
Definition: Controller.h:136
BigReal alchTemp
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
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:3944
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:114
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 3023 of file Controller.C.

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

Referenced by integrate().

3023  {
3024 
3025  Node *node = Node::Object();
3026  Molecule *molecule = node->molecule;
3027  SimParameters *simParameters = node->simParameters;
3028  Lattice &lattice = state->lattice;
3029 
3030  compareChecksums(step);
3031 
3032  printEnergies(step,0);
3033 }
static Node * Object()
Definition: Node.h:86
void compareChecksums(int, int=0)
Definition: Controller.C:2775
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:3035
Molecule * molecule
Definition: Node.h:176
void Controller::printEnergies ( int  step,
int  minimize 
)
protected

Definition at line 3035 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::outputEnergiesPrecision, 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().

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