Controller Class Reference

#include <Controller.h>

Inheritance diagram for Controller:
ControllerState

List of all members.

Classes

struct  checkpoint

Public Member Functions

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

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

Friends

class ScriptTcl
class Node
class CheckpointMsg

Detailed Description

Definition at line 39 of file Controller.h.


Constructor & Destructor Documentation

Controller::Controller ( NamdState s  ) 

Definition at line 128 of file Controller.C.

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

00128                                    :
00129         computeChecksum(0), marginViolations(0), pairlistWarnings(0),
00130         simParams(Node::Object()->simParameters),
00131         state(s),
00132         collection(CollectionMaster::Object()),
00133         startCTime(0),
00134         firstCTime(CmiTimer()),
00135         startWTime(0),
00136         firstWTime(CmiWallTimer()),
00137         startBenchTime(0),
00138         stepInFullRun(0),
00139         ldbSteps(0),
00140         fflush_count(3)
00141 {
00142     broadcast = new ControllerBroadcasts;
00143     min_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MINIMIZER,1);
00144     // for accelMD
00145     if (simParams->accelMDOn) {
00146        // REDUCTIONS_BASIC wil contain all potential energies and arrive first
00147        amd_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
00148        // contents of amd_reduction be submitted to REDUCTIONS_AMD
00149        submit_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_AMD);
00150        // REDUCTIONS_AMD will contain Sequencer contributions and arrive second
00151        reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_AMD);
00152     } else {
00153        amd_reduction = NULL;
00154        submit_reduction = NULL;
00155        reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
00156     }
00157     // pressure profile reductions
00158     pressureProfileSlabs = 0;
00159     pressureProfileCount = 0;
00160     ppbonded = ppnonbonded = ppint = NULL;
00161     if (simParams->pressureProfileOn || simParams->pressureProfileEwaldOn) {
00162       int ntypes = simParams->pressureProfileAtomTypes;
00163       int nslabs = pressureProfileSlabs = simParams->pressureProfileSlabs;
00164       // number of partitions for pairwise interactions
00165       int npairs = (ntypes * (ntypes+1))/2;
00166       pressureProfileAverage = new BigReal[3*nslabs];
00167       memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
00168       int freq = simParams->pressureProfileFreq;
00169       if (simParams->pressureProfileOn) {
00170         ppbonded = new PressureProfileReduction(REDUCTIONS_PPROF_BONDED, 
00171             nslabs, npairs, "BONDED", freq);
00172         ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
00173             nslabs, npairs, "NONBONDED", freq);
00174         // internal partitions by atom type, but not in a pairwise fashion
00175         ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL, 
00176             nslabs, ntypes, "INTERNAL", freq);
00177       } else {
00178         // just doing Ewald, so only need nonbonded
00179         ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
00180             nslabs, npairs, "NONBONDED", freq);
00181       }
00182     }
00183     random = new Random(simParams->randomSeed);
00184     random->split(0,PatchMap::Object()->numPatches()+1);
00185 
00186     heat = totalEnergy0 = 0.0;
00187     rescaleVelocities_sumTemps = 0;  
00188     rescaleVelocities_numTemps = 0;
00189     stochRescale_count = 0;
00190     if (simParams->stochRescaleOn) {
00191       stochRescaleTimefactor = exp(-simParams->stochRescaleFreq *
00192           simParams->dt * 0.001 / simParams->stochRescalePeriod);
00193     }
00194     berendsenPressure_avg = 0; berendsenPressure_count = 0;
00195     // strainRate tensor is symmetric to avoid rotation
00196     langevinPiston_strainRate =
00197         Tensor::symmetric(simParams->strainRate,simParams->strainRate2);
00198     if ( ! simParams->useFlexibleCell ) {
00199       BigReal avg = trace(langevinPiston_strainRate) / 3.;
00200       langevinPiston_strainRate = Tensor::identity(avg);
00201     } else if ( simParams->useConstantRatio ) {
00202 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
00203                  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
00204       AVGXY(langevinPiston_strainRate);
00205 #undef AVGXY
00206     }
00207     langevinPiston_origStrainRate = langevinPiston_strainRate;
00208     if (simParams->multigratorOn) {
00209       multigratorXi = 0.0;
00210       int n = simParams->multigratorNoseHooverChainLength;
00211       BigReal tau = simParams->multigratorTemperatureRelaxationTime;
00212       Node *node = Node::Object();
00213       Molecule *molecule = node->molecule;
00214       BigReal Nf = molecule->num_deg_freedom();
00215       BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
00216       multigratorNu.resize(n);
00217       multigratorNuT.resize(n);
00218       multigratorZeta.resize(n);
00219       multigratorOmega.resize(n);
00220       for (int i=0;i < n;i++) {
00221         multigratorNu[i] = 0.0;
00222         multigratorZeta[i] = 0.0;
00223         if (i == 0) {
00224           multigratorOmega[i] = Nf*kT0*tau*tau;
00225         } else {
00226           multigratorOmega[i] = kT0*tau*tau;
00227         }
00228       }
00229       multigratorReduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
00230     } else {
00231       multigratorReduction = NULL;
00232     }
00233     origLattice = state->lattice;
00234     smooth2_avg = XXXBIGREAL;
00235     temp_avg = 0;
00236     pressure_avg = 0;
00237     groupPressure_avg = 0;
00238     avg_count = 0;
00239     pressure_tavg = 0;
00240     groupPressure_tavg = 0;
00241     tavg_count = 0;
00242     checkpoint_stored = 0;
00243     drudeBondTemp = 0;
00244     drudeBondTempAvg = 0;
00245     cumAlchWork = 0;
00246 }

Controller::~Controller ( void   )  [virtual]

Definition at line 248 of file Controller.C.

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

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


Member Function Documentation

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

Definition at line 2343 of file Controller.C.

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

Referenced by adaptTempUpdate().

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

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

Definition at line 2469 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 2443 of file Controller.C.

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

Referenced by adaptTempUpdate().

02443                                                {
02444     if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
02445         adaptTempRestartFile.seekp(std::ios::beg);        
02446         iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
02447         adaptTempRestartFile << step << " ";
02448         // Start with min and max temperatures
02449         adaptTempRestartFile << adaptTempT<< " ";     // KELVIN
02450         adaptTempRestartFile << 1./adaptTempBetaMin << " ";  // KELVIN
02451         adaptTempRestartFile << 1./adaptTempBetaMax << " ";  // KELVIN
02452         adaptTempRestartFile << adaptTempBins << " ";     
02453         adaptTempRestartFile << adaptTempCg << " ";
02454         adaptTempRestartFile << adaptTempDt ;
02455         adaptTempRestartFile << "\n" ;
02456         for(int j = 0; j < adaptTempBins; ++j) {
02457           adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
02458           adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
02459           adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
02460           adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
02461           adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
02462           adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
02463           adaptTempRestartFile << "\n";          
02464         }
02465         adaptTempRestartFile.flush(); 
02466     }
02467 }    

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

Definition at line 281 of file Controller.C.

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

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

void Controller::awaken ( void   )  [inline]

Definition at line 45 of file Controller.h.

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

00045 { CthAwaken(thread); };

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

Definition at line 950 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 1728 of file Controller.C.

Referenced by rescaleaccelMD().

01729                                                               {
01730    switch(iE){
01731       case 2:
01732          *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
01733          // if k0 <=0 OR k0 > 1, switch to iE=1 mode
01734          if(*k0 > 1.)
01735             *warn |= 1;
01736          else if(*k0 <= 0.)
01737             *warn |= 2;
01738          //else stay in iE=2 mode
01739          else{
01740             *E = Vmin + (Vmax-Vmin)/(*k0);
01741             *iEused = 2;
01742             break;
01743          }
01744       case 1:
01745          *E = Vmax;
01746          *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
01747          if(*k0 > 1.)
01748             *k0 = 1.;
01749          *iEused = 1;
01750          break;
01751    }
01752 
01753    *k = *k0 / (Vmax - Vmin);
01754 }

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

Definition at line 1757 of file Controller.C.

Referenced by rescaleaccelMD().

01758                                            {
01759    BigReal VE = testV - E;
01760    //if V < E, apply boost
01761    if(VE < 0){
01762       *factor = k * VE;
01763       *vir += vir_orig * (*factor);
01764       *dV += (*factor) * VE/2.;
01765       *factor += 1.;
01766    }
01767    else{
01768       *factor = 1.;
01769    }
01770 }

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

Definition at line 1708 of file Controller.C.

Referenced by rescaleaccelMD().

01709                                                                            {
01710    BigReal delta; 
01711 
01712    if(testV > *Vmax){
01713       *Vmax = testV;
01714    }
01715    else if(testV < *Vmin){
01716       *Vmin = testV;
01717    }
01718 
01719    //mean and std calculated by Online Algorithm
01720    delta = testV - *Vavg;
01721    *Vavg += delta / (BigReal)n;
01722    *M2 += delta * (testV - (*Vavg));
01723 
01724    *sigmaV = sqrt(*M2/n);
01725 }

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

Definition at line 1590 of file Controller.C.

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

Referenced by multigratorPressure(), and receivePressure().

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

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

Definition at line 2767 of file Controller.C.

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

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

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

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

Definition at line 3870 of file Controller.C.

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

Referenced by printEnergies().

03870                                                   {
03871   // alchemical scaling factors for groups 1/2 at the previous lambda
03872   const BigReal oldLambda = simParams->getCurrentLambda(step-1);
03873   const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
03874   const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
03875   const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
03876   const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
03877   const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
03878   const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
03879   // alchemical scaling factors for groups 1/2 at the new lambda
03880   const BigReal alchLambda = simParams->getCurrentLambda(step);
03881   const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
03882   const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
03883   const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
03884   const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
03885   const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
03886   const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda); 
03887 
03888   return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
03889           (elec_lambda_12 - elec_lambda_1)*
03890           (electEnergy_ti_1 + electEnergySlow_ti_1 +  electEnergyPME_ti_1) +
03891           (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
03892           (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
03893           (elec_lambda_22 - elec_lambda_2)*
03894           (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) + 
03895           (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
03896   );
03897 }

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

Definition at line 1255 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 4131 of file Controller.C.

References broadcast.

Referenced by integrate().

04131                                                      {
04132 #if USE_BARRIER
04133         if (doBarrier) {
04134           broadcast->cycleBarrier.publish(step,1);
04135           CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
04136                   CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04137         }
04138 #endif
04139 }

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

Definition at line 429 of file Controller.C.

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

Referenced by algorithm().

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

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

Definition at line 987 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 1142 of file Controller.C.

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

Referenced by integrate().

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

void Controller::minimize (  )  [protected]

Definition at line 594 of file Controller.C.

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

Referenced by algorithm().

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

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

Definition at line 919 of file Controller.C.

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

Referenced by printEnergies().

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

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

Definition at line 778 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 853 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 3959 of file Controller.C.

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

Referenced by algorithm(), and integrate().

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

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

Definition at line 3667 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 3741 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 3015 of file Controller.C.

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

Referenced by integrate().

03015                                                {
03016 
03017     Node *node = Node::Object();
03018     Molecule *molecule = node->molecule;
03019     SimParameters *simParameters = node->simParameters;
03020     Lattice &lattice = state->lattice;
03021 
03022     compareChecksums(step);
03023 
03024     printEnergies(step,0);
03025 }

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

Definition at line 3027 of file Controller.C.

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

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

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

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

Definition at line 1274 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 1294 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 2959 of file Controller.C.

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

Referenced by printEnergies().

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

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

Definition at line 1308 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 4117 of file Controller.C.

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

Referenced by integrate().

04118 {
04119   if ( ! ldbSteps ) { 
04120     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
04121   }
04122   if ( ! --ldbSteps ) {
04123     startBenchTime -= CmiWallTimer();
04124         Node::Object()->outputPatchComputeMaps("before_ldb", step);
04125     LdbCoordinator::Object()->rebalance(this);  
04126         startBenchTime += CmiWallTimer();
04127     fflush_count = 3;
04128   }
04129 }

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

Definition at line 1412 of file Controller.C.

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

Referenced by integrate(), and printMinimizeEnergies().

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

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

Definition at line 4108 of file Controller.C.

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

Referenced by Node::recvCheckpointAck().

04108                                                  {  // initiating replica
04109   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
04110     state->lattice = cp.lattice;
04111     *(ControllerState*)this = cp.state;
04112   }
04113   CkpvAccess(_qd)->process();
04114 }

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

Definition at line 4078 of file Controller.C.

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

Referenced by Node::recvCheckpointReq().

04078                                                                             {  // responding replica
04079   switch ( task ) {
04080     case SCRIPT_CHECKPOINT_STORE:
04081       if ( ! checkpoints.count(key) ) {
04082         checkpoints[key] = new checkpoint;
04083       }
04084       *checkpoints[key] = cp;
04085       break;
04086     case SCRIPT_CHECKPOINT_LOAD:
04087       if ( ! checkpoints.count(key) ) {
04088         NAMD_die("Unable to load checkpoint, requested key was never stored.");
04089       }
04090       cp = *checkpoints[key];
04091       break;
04092     case SCRIPT_CHECKPOINT_SWAP:
04093       if ( ! checkpoints.count(key) ) {
04094         NAMD_die("Unable to swap checkpoint, requested key was never stored.");
04095       }
04096       std::swap(cp,*checkpoints[key]);
04097       break;
04098     case SCRIPT_CHECKPOINT_FREE:
04099       if ( ! checkpoints.count(key) ) {
04100         NAMD_die("Unable to free checkpoint, requested key was never stored.");
04101       }
04102       delete checkpoints[key];
04103       checkpoints.erase(key);
04104       break;
04105   }
04106 }

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

Definition at line 1831 of file Controller.C.

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

Referenced by integrate(), and printMinimizeEnergies().

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

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

Definition at line 1233 of file Controller.C.

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

Referenced by integrate().

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

void Controller::resumeAfterTraceBarrier ( int  step  ) 

Definition at line 4148 of file Controller.C.

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

Referenced by Node::resumeAfterTraceBarrier().

04148                                                 {
04149         broadcast->traceBarrier.publish(step,1);
04150         CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime); 
04151         awaken();
04152 }

void Controller::run ( void   ) 

Definition at line 268 of file Controller.C.

References awaken(), CTRL_STK_SZ, and DebugM.

Referenced by NamdState::runController().

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

double Controller::stochRescaleCoefficient (  )  [protected]

Calculate new coefficient for stochastic velocity rescaling and update heat.

Definition at line 1364 of file Controller.C.

References BOLTZMANN, Random::gaussian(), heat, numDegFreedom, random, simParams, SimParameters::stochRescaleTemp, stochRescaleTimefactor, Random::sum_of_squared_gaussians(), and temperature.

Referenced by stochRescaleVelocities().

01364                                            {
01365   const double stochRescaleTemp = simParams->stochRescaleTemp;
01366   double coefficient = 1;
01367   if ( temperature > 0 ) {
01368     double R1 = random->gaussian();
01369     // double gammaShape = 0.5*(numDegFreedom - 1);
01370     // R2sum is the sum of (numDegFreedom - 1) squared normal variables,
01371     // which is chi-squared distributed.
01372     // This is in turn a special case of the Gamma distribution,
01373     // which converges to a normal distribution in the limit of a
01374     // large shape parameter.
01375     // double R2sum = 2*(gammaShape+sqrt(gammaShape)*random->gaussian())+R1*R1;
01376     double R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
01377     double tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
01378 
01379     coefficient = sqrt(stochRescaleTimefactor +
01380         (1 - stochRescaleTimefactor)*tempfactor*R2sum +
01381         2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*
01382           stochRescaleTimefactor)); 
01383   }
01384   heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient-1);
01385   return coefficient;
01386 }

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

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

Generate and broadcast the scale factor for stochastic velocity rescaling.

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

Parameters:
step the current timestep

Definition at line 1348 of file Controller.C.

References broadcast, SimpleBroadcastObject< T >::publish(), simParams, stochRescale_count, ControllerBroadcasts::stochRescaleCoefficient, stochRescaleCoefficient(), SimParameters::stochRescaleFreq, and SimParameters::stochRescaleOn.

Referenced by integrate().

01349 {
01350   if ( simParams->stochRescaleOn ) {
01351     ++stochRescale_count;
01352     if ( stochRescale_count == simParams->stochRescaleFreq ) {
01353       double coefficient = stochRescaleCoefficient();
01354       broadcast->stochRescaleCoefficient.publish(step,coefficient);
01355       stochRescale_count = 0;
01356     }
01357   }
01358 }

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

Definition at line 1326 of file Controller.C.

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

Referenced by integrate().

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

void Controller::terminate ( void   )  [protected]

Definition at line 4169 of file Controller.C.

References awaken().

Referenced by algorithm().

04169                                {
04170   BackEnd::awaken();
04171   CthFree(thread);
04172   CthSuspend();
04173 }

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

Definition at line 4141 of file Controller.C.

Referenced by integrate().

04141                                                        {
04142         CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);       
04143         CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04144         nd.traceBarrier(turnOnTrace, step);
04145         CthSuspend();
04146 }

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

Definition at line 1773 of file Controller.C.

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

Referenced by rescaleaccelMD().

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

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

Definition at line 3634 of file Controller.C.

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

Referenced by outputExtendedSystem().

03634                                                                       {
03635   Lattice &lattice = state->lattice;
03636   file.precision(12);
03637   file << step;
03638     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
03639     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
03640     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
03641     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
03642   if ( simParams->langevinPistonOn ) {
03643     Vector strainRate = diagonal(langevinPiston_strainRate);
03644     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
03645     file << " " << strainRate.x;
03646     file << " " << strainRate.y;
03647     file << " " << strainRate.z;
03648     file << " " << strainRate2.x;
03649     file << " " << strainRate2.y;
03650     file << " " << strainRate2.z;
03651   }
03652   file << std::endl;
03653 }

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

Definition at line 3621 of file Controller.C.

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

Referenced by outputExtendedSystem().

03621                                                               {
03622   Lattice &lattice = state->lattice;
03623   file << "#$LABELS step";
03624   if ( lattice.a_p() ) file << " a_x a_y a_z";
03625   if ( lattice.b_p() ) file << " b_x b_y b_z";
03626   if ( lattice.c_p() ) file << " c_x c_y c_z";
03627   file << " o_x o_y o_z";
03628   if ( simParams->langevinPistonOn ) {
03629     file << " s_x s_y s_z s_u s_v s_w";
03630   }
03631   file << std::endl;
03632 }

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

Definition at line 3899 of file Controller.C.

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

Referenced by outputFepEnergy().

03899                                                                  {
03900   BigReal eeng = electEnergy + electEnergySlow;
03901   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
03902   BigReal RT = BOLTZMANN * simParams->alchTemp;
03903 
03904   const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03905   const int stepInRun = step - simParams->firstTimestep;
03906   const BigReal alchLambda = simParams->alchLambda;
03907 
03908   if(stepInRun){
03909     if ( simParams->alchLambdaIDWS >= 0. &&
03910         (step / simParams->alchIDWSFreq) % 2 == 0 ) {
03911       // IDWS is active and we are on a "backward" timestep
03912       fepFile << FEPTITLE_BACK(step);
03913     } else {
03914       // "forward" timestep
03915       fepFile << FEPTITLE(step);
03916     }
03917     fepFile << FORMAT(eeng);
03918     fepFile << FORMAT(eeng_f);
03919     fepFile << FORMAT(ljEnergy);
03920     fepFile << FORMAT(ljEnergy_f);
03921     fepFile << FORMAT(dE);
03922     if(alchEnsembleAvg){
03923       BigReal dE_avg = net_dE / FepNo;
03924       fepFile << FORMAT(dE_avg);
03925     }
03926     fepFile << FORMAT(temperature);
03927     if(alchEnsembleAvg){
03928       dG = -(RT * log(exp_dE_ByRT / FepNo));
03929       fepFile << FORMAT(dG);
03930     }
03931     fepFile << std::endl;
03932   }
03933 }

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

Friends And Related Function Documentation

friend class CheckpointMsg [friend]

Definition at line 55 of file Controller.h.

friend class Node [friend]

Definition at line 54 of file Controller.h.

friend class ScriptTcl [friend]

Definition at line 53 of file Controller.h.


Member Data Documentation

Definition at line 50 of file Controller.h.

Referenced by rescaleaccelMD().

Definition at line 301 of file Controller.h.

Referenced by rescaleaccelMD().

Definition at line 324 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 318 of file Controller.h.

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

Definition at line 317 of file Controller.h.

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

Definition at line 313 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

int Controller::adaptTempBin [protected]

Definition at line 319 of file Controller.h.

Referenced by adaptTempUpdate().

int Controller::adaptTempBins [protected]

Definition at line 320 of file Controller.h.

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

Definition at line 322 of file Controller.h.

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

Definition at line 321 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 323 of file Controller.h.

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

Definition at line 315 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 316 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 326 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 325 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 310 of file Controller.h.

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

Definition at line 308 of file Controller.h.

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

Definition at line 307 of file Controller.h.

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

Definition at line 312 of file Controller.h.

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

Definition at line 311 of file Controller.h.

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

Definition at line 309 of file Controller.h.

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

Definition at line 327 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempWriteRestart().

Definition at line 151 of file Controller.h.

Referenced by outputTiEnergy(), and printEnergies().

Definition at line 238 of file Controller.h.

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

int Controller::avg_count [protected]

Definition at line 84 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 126 of file Controller.h.

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

Definition at line 127 of file Controller.h.

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

Definition at line 113 of file Controller.h.

Referenced by outputFepEnergy(), and printEnergies().

Definition at line 269 of file Controller.h.

Referenced by algorithm().

Definition at line 270 of file Controller.h.

Referenced by algorithm().

Definition at line 268 of file Controller.h.

Referenced by algorithm(), and Controller().

int Controller::checkpoint_task [protected]

Definition at line 277 of file Controller.h.

Referenced by recvCheckpointAck().

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

Definition at line 276 of file Controller.h.

Referenced by algorithm(), and recvCheckpointReq().

Definition at line 249 of file Controller.h.

Referenced by enqueueCollections().

int Controller::computeChecksum [protected]

Definition at line 89 of file Controller.h.

Referenced by compareChecksums().

Definition at line 77 of file Controller.h.

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

Definition at line 76 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 78 of file Controller.h.

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

Definition at line 140 of file Controller.h.

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

BigReal Controller::dE [protected]

Definition at line 119 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::dG [protected]

Definition at line 121 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 155 of file Controller.h.

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

Definition at line 156 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 104 of file Controller.h.

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

Definition at line 114 of file Controller.h.

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

Definition at line 128 of file Controller.h.

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

Definition at line 131 of file Controller.h.

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

Definition at line 141 of file Controller.h.

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

Definition at line 142 of file Controller.h.

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

Definition at line 105 of file Controller.h.

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

Definition at line 115 of file Controller.h.

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

Definition at line 129 of file Controller.h.

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

Definition at line 132 of file Controller.h.

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

Definition at line 118 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 258 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

int Controller::FepNo [protected]

Definition at line 122 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 124 of file Controller.h.

Referenced by outputFepEnergy().

int Controller::fflush_count [protected]

Definition at line 222 of file Controller.h.

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

Definition at line 109 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 110 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 111 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 108 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 107 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 168 of file Controller.h.

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

Definition at line 83 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 74 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 73 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 75 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 86 of file Controller.h.

Referenced by Controller(), and printEnergies().

heat exchanged with the thermostat since firstTimestep

Definition at line 163 of file Controller.h.

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

Definition at line 160 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Definition at line 159 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Definition at line 204 of file Controller.h.

Referenced by Controller().

int Controller::ldbSteps [protected]

Definition at line 220 of file Controller.h.

Referenced by rebalanceLoad().

Definition at line 106 of file Controller.h.

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

Definition at line 116 of file Controller.h.

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

Definition at line 117 of file Controller.h.

Definition at line 130 of file Controller.h.

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

Definition at line 133 of file Controller.h.

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

Definition at line 90 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Definition at line 94 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 95 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 96 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

int Controller::min_huge_count [protected]

Definition at line 98 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 60 of file Controller.h.

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

Definition at line 97 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 211 of file Controller.h.

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

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

Definition at line 213 of file Controller.h.

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

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

Definition at line 214 of file Controller.h.

Referenced by Controller(), and multigratorTemperature().

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

Definition at line 215 of file Controller.h.

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

Definition at line 217 of file Controller.h.

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

Definition at line 209 of file Controller.h.

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

Definition at line 210 of file Controller.h.

Referenced by multigratorPressure().

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

Definition at line 216 of file Controller.h.

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

int Controller::nbondFreq [protected]

Definition at line 120 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 134 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 135 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 136 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 137 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 138 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 139 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int64_t Controller::numDegFreedom [protected]

Definition at line 281 of file Controller.h.

Referenced by Controller().

Definition at line 91 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Definition at line 206 of file Controller.h.

Referenced by langevinPiston1().

Definition at line 242 of file Controller.h.

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

Definition at line 244 of file Controller.h.

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

Definition at line 243 of file Controller.h.

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

Definition at line 167 of file Controller.h.

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

Definition at line 71 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 82 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 69 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 68 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 70 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 85 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 247 of file Controller.h.

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

Definition at line 246 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 245 of file Controller.h.

Referenced by Controller(), and printEnergies().

Random* Controller::random [protected]

Definition at line 150 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 144 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 145 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 146 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 147 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 148 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 149 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::recent_TiNo [protected]

Definition at line 152 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 175 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

Definition at line 174 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

int Controller::slowFreq [protected]

Definition at line 166 of file Controller.h.

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

Definition at line 102 of file Controller.h.

Referenced by printEnergies().

Count time steps until next stochastic velocity rescaling.

Definition at line 193 of file Controller.h.

Referenced by Controller(), and stochRescaleVelocities().

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

Definition at line 196 of file Controller.h.

Referenced by Controller(), and stochRescaleCoefficient().

Definition at line 205 of file Controller.h.

Referenced by langevinPiston1().

Definition at line 239 of file Controller.h.

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

int Controller::tavg_count [protected]

Definition at line 87 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 81 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 262 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::TiNo [protected]

Definition at line 143 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 103 of file Controller.h.

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

totalEnergy at firstTimestep

Definition at line 164 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 72 of file Controller.h.

Referenced by calcPressure(), and rescaleaccelMD().

Definition at line 252 of file Controller.h.

Referenced by outputExtendedSystem().


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

Generated on 3 Aug 2020 for NAMD by  doxygen 1.6.1