Controller Class Reference

#include <Controller.h>

Inheritance diagram for Controller:
ControllerState

List of all members.

Classes

struct  checkpoint

Public Member Functions

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

Protected Member Functions

virtual void algorithm (void)
void integrate (int)
void minimize ()
void receivePressure (int step, int minimize=0)
void calcPressure (int step, int minimize, const Tensor &virial_normal_in, const Tensor &virial_nbond_in, const Tensor &virial_slow_in, const Tensor &intVirial_normal, const Tensor &intVirial_nbond, const Tensor &intVirial_slow, const Vector &extForce_normal, const Vector &extForce_nbond, const Vector &extForce_slow)
void compareChecksums (int, int=0)
void printTiming (int)
void printMinimizeEnergies (int)
void printDynamicsEnergies (int)
void printEnergies (int step, int minimize)
void printFepMessage (int)
void printTiMessage (int)
void enqueueCollections (int)
void correctMomentum (int step)
void rescaleVelocities (int)
void reassignVelocities (int)
void tcoupleVelocities (int)
void stochRescaleVelocities (int)
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 2339 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().

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

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

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

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

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

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

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

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

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

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

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

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

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

References broadcast.

Referenced by integrate().

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

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

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

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

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

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

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

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

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

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

Definition at line 3011 of file Controller.C.

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

Referenced by integrate().

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

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

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

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

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

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

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

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

Referenced by integrate().

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

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

04104                                                  {  // initiating replica
04105   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
04106     state->lattice = cp.lattice;
04107     *(ControllerState*)this = cp.state;
04108   }
04109   CkpvAccess(_qd)->process();
04110 }

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

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

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

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

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

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

Referenced by Node::resumeAfterTraceBarrier().

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

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

References awaken().

Referenced by algorithm().

04165                                {
04166   BackEnd::awaken();
04167   CthFree(thread);
04168   CthSuspend();
04169 }

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

Definition at line 4137 of file Controller.C.

Referenced by integrate().

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

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

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

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

Definition at line 3617 of file Controller.C.

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

Referenced by outputExtendedSystem().

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

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

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

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

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

Friends And Related Function Documentation

friend class CheckpointMsg [friend]

Definition at line 54 of file Controller.h.

friend class Node [friend]

Definition at line 53 of file Controller.h.

friend class ScriptTcl [friend]

Definition at line 52 of file Controller.h.


Member Data Documentation

Definition at line 300 of file Controller.h.

Referenced by rescaleaccelMD().

Definition at line 323 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 317 of file Controller.h.

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

Definition at line 316 of file Controller.h.

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

Definition at line 312 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

int Controller::adaptTempBin [protected]

Definition at line 318 of file Controller.h.

Referenced by adaptTempUpdate().

int Controller::adaptTempBins [protected]

Definition at line 319 of file Controller.h.

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

Definition at line 321 of file Controller.h.

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

Definition at line 320 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 322 of file Controller.h.

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

Definition at line 314 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 315 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 324 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

Definition at line 309 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 306 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 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 326 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempWriteRestart().

Definition at line 150 of file Controller.h.

Referenced by outputTiEnergy(), and printEnergies().

Definition at line 237 of file Controller.h.

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

int Controller::avg_count [protected]

Definition at line 83 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 125 of file Controller.h.

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

Definition at line 126 of file Controller.h.

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

Definition at line 112 of file Controller.h.

Referenced by outputFepEnergy(), and printEnergies().

Definition at line 268 of file Controller.h.

Referenced by algorithm().

Definition at line 269 of file Controller.h.

Referenced by algorithm().

Definition at line 267 of file Controller.h.

Referenced by algorithm(), and Controller().

int Controller::checkpoint_task [protected]

Definition at line 276 of file Controller.h.

Referenced by recvCheckpointAck().

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

Definition at line 275 of file Controller.h.

Referenced by algorithm(), and recvCheckpointReq().

Definition at line 248 of file Controller.h.

Referenced by enqueueCollections().

int Controller::computeChecksum [protected]

Definition at line 88 of file Controller.h.

Referenced by compareChecksums().

Definition at line 76 of file Controller.h.

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

Definition at line 75 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 77 of file Controller.h.

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

Definition at line 139 of file Controller.h.

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

BigReal Controller::dE [protected]

Definition at line 118 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::dG [protected]

Definition at line 120 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 154 of file Controller.h.

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

Definition at line 155 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 103 of file Controller.h.

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

Definition at line 113 of file Controller.h.

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

Definition at line 127 of file Controller.h.

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

Definition at line 130 of file Controller.h.

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

Definition at line 140 of file Controller.h.

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

Definition at line 141 of file Controller.h.

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

Definition at line 104 of file Controller.h.

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

Definition at line 114 of file Controller.h.

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

Definition at line 128 of file Controller.h.

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

Definition at line 131 of file Controller.h.

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

Definition at line 117 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 257 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

int Controller::FepNo [protected]

Definition at line 121 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 123 of file Controller.h.

Referenced by outputFepEnergy().

int Controller::fflush_count [protected]

Definition at line 221 of file Controller.h.

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

Definition at line 108 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 109 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 110 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 107 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 106 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Definition at line 167 of file Controller.h.

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

Definition at line 82 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 73 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 72 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 74 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 85 of file Controller.h.

Referenced by Controller(), and printEnergies().

heat exchanged with the thermostat since firstTimestep

Definition at line 162 of file Controller.h.

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

Definition at line 159 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Definition at line 158 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Definition at line 203 of file Controller.h.

Referenced by Controller().

int Controller::ldbSteps [protected]

Definition at line 219 of file Controller.h.

Referenced by rebalanceLoad().

Definition at line 105 of file Controller.h.

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

Definition at line 115 of file Controller.h.

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

Definition at line 116 of file Controller.h.

Definition at line 129 of file Controller.h.

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

Definition at line 132 of file Controller.h.

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

Definition at line 89 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Definition at line 93 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 94 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 95 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

int Controller::min_huge_count [protected]

Definition at line 97 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 59 of file Controller.h.

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

Definition at line 96 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Definition at line 210 of file Controller.h.

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

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

Definition at line 212 of file Controller.h.

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

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

Definition at line 213 of file Controller.h.

Referenced by Controller(), and multigratorTemperature().

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

Definition at line 214 of file Controller.h.

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

Definition at line 216 of file Controller.h.

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

Definition at line 208 of file Controller.h.

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

Definition at line 209 of file Controller.h.

Referenced by multigratorPressure().

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

Definition at line 215 of file Controller.h.

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

int Controller::nbondFreq [protected]

Definition at line 119 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

Definition at line 133 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 134 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 135 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 136 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 137 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 138 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int64_t Controller::numDegFreedom [protected]

Definition at line 280 of file Controller.h.

Referenced by Controller().

Definition at line 90 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Definition at line 205 of file Controller.h.

Referenced by langevinPiston1().

Definition at line 241 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 242 of file Controller.h.

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

Definition at line 166 of file Controller.h.

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

Definition at line 70 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 81 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 68 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 67 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 69 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Definition at line 84 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 246 of file Controller.h.

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

Definition at line 245 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 244 of file Controller.h.

Referenced by Controller(), and printEnergies().

Random* Controller::random [protected]

Definition at line 149 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 143 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 144 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 145 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 146 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 147 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 148 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::recent_TiNo [protected]

Definition at line 151 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 174 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

Definition at line 173 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

int Controller::slowFreq [protected]

Definition at line 165 of file Controller.h.

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

Definition at line 101 of file Controller.h.

Referenced by printEnergies().

Count time steps until next stochastic velocity rescaling.

Definition at line 192 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 195 of file Controller.h.

Referenced by Controller(), and stochRescaleCoefficient().

Definition at line 204 of file Controller.h.

Referenced by langevinPiston1().

Definition at line 238 of file Controller.h.

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

int Controller::tavg_count [protected]

Definition at line 86 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 80 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 261 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::TiNo [protected]

Definition at line 142 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

Definition at line 102 of file Controller.h.

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

totalEnergy at firstTimestep

Definition at line 163 of file Controller.h.

Referenced by Controller(), and printEnergies().

Definition at line 71 of file Controller.h.

Referenced by calcPressure(), and rescaleaccelMD().

Definition at line 251 of file Controller.h.

Referenced by outputExtendedSystem().


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

Generated on 6 Dec 2019 for NAMD by  doxygen 1.6.1