Controller Class Reference

#include <Controller.h>

Inheritance diagram for Controller:

ControllerState List of all members.

Public Member Functions

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

Protected Member Functions

virtual void algorithm (void)
void integrate (int)
void minimize ()
void receivePressure (int step, int minimize=0)
void calcPressure (int step, int minimize, const Tensor &virial_normal_in, const Tensor &virial_nbond_in, const Tensor &virial_slow_in, const Tensor &intVirial_normal, const Tensor &intVirial_nbond, const Tensor &intVirial_slow, const Vector &extForce_normal, const Vector &extForce_nbond, const Vector &extForce_slow)
void compareChecksums (int, int=0)
void printTiming (int)
void printMinimizeEnergies (int)
void printDynamicsEnergies (int)
void printEnergies (int step, int minimize)
void printFepMessage (int)
void printTiMessage (int)
void enqueueCollections (int)
void correctMomentum (int step)
void rescaleVelocities (int)
void reassignVelocities (int)
void tcoupleVelocities (int)
void berendsenPressure (int)
void langevinPiston1 (int)
void langevinPiston2 (int)
void multigratorPressure (int step, int callNumber)
void multigratorTemperature (int step, int callNumber)
BigReal multigatorCalcEnthalpy (BigReal potentialEnergy, int step, int minimize)
void rebalanceLoad (int)
void cycleBarrier (int, int)
void traceBarrier (int, int)
void terminate (void)
void outputExtendedSystem (int step)
void writeExtendedSystemLabels (ofstream_namd &file)
void writeExtendedSystemData (int step, ofstream_namd &file)
void outputFepEnergy (int step)
void writeFepEnergyData (int step, ofstream_namd &file)
void outputTiEnergy (int step)
BigReal computeAlchWork (const int step)
void writeTiEnergyData (int step, ofstream_namd &file)
void recvCheckpointReq (const char *key, int task, checkpoint &cp)
void recvCheckpointAck (checkpoint &cp)
void calc_accelMDG_mean_std (BigReal testV, int step_n, BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV)
void calc_accelMDG_E_k (int iE, int V_n, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal *k0, BigReal *k, BigReal *E, int *iEused, char *warn)
void calc_accelMDG_force_factor (BigReal k, BigReal E, BigReal testV, Tensor vir_orig, BigReal *dV, BigReal *factor, Tensor *vir)
void write_accelMDG_rest_file (int step_n, char type, int V_n, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal M2, BigReal E, BigReal k, bool write_topic, bool lasttime)
void rescaleaccelMD (int step, int minimize=0)
void adaptTempInit (int step)
void adaptTempUpdate (int step, int minimize=0)
void adaptTempWriteRestart (int step)

Protected Attributes

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

Friends

class ScriptTcl
class Node
class CheckpointMsg

Classes

struct  checkpoint

Detailed Description

Definition at line 39 of file Controller.h.


Constructor & Destructor Documentation

Controller::Controller ( NamdState s  ) 

Definition at line 128 of file Controller.C.

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

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

Controller::~Controller ( void   )  [virtual]

Definition at line 241 of file Controller.C.

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

00242 {
00243     delete broadcast;
00244     delete reduction;
00245     delete min_reduction;
00246     delete amd_reduction;
00247     delete submit_reduction;
00248     delete ppbonded;
00249     delete ppnonbonded;
00250     delete ppint;
00251     delete [] pressureProfileAverage;
00252     delete random;
00253     if (multigratorReduction) delete multigratorReduction;
00254 }


Member Function Documentation

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

Definition at line 2231 of file Controller.C.

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

Referenced by adaptTempUpdate().

02231                                        {
02232     if (!simParams->adaptTempOn) return;
02233     iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
02234     adaptTempDtMin = 0;
02235     adaptTempDtMax = 0;
02236     adaptTempAutoDt = false;
02237     if (simParams->adaptTempBins == 0) {
02238       iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
02239       std::ifstream adaptTempRead(simParams->adaptTempInFile);
02240       if (adaptTempRead) {
02241         int readInt;
02242         BigReal readReal;
02243         bool readBool;
02244         // step
02245         adaptTempRead >> readInt;
02246         // Start with min and max temperatures
02247         adaptTempRead >> adaptTempT;     // KELVIN
02248         adaptTempRead >> adaptTempBetaMin;  // KELVIN
02249         adaptTempRead >> adaptTempBetaMax;  // KELVIN
02250         adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
02251         adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
02252         // In case file is manually edited
02253         if (adaptTempBetaMin > adaptTempBetaMax){
02254             readReal = adaptTempBetaMax;
02255             adaptTempBetaMax = adaptTempBetaMin;
02256             adaptTempBetaMin = adaptTempBetaMax;
02257         }
02258         adaptTempRead >> adaptTempBins;     
02259         adaptTempRead >> adaptTempCg;
02260         adaptTempRead >> adaptTempDt;
02261         adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02262         adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02263         adaptTempPotEnergySamples = new int[adaptTempBins];
02264         adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02265         adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02266         adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02267         adaptTempBetaN           = new BigReal[adaptTempBins];
02268         adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02269         for(int j = 0; j < adaptTempBins; ++j) {
02270           adaptTempRead >> adaptTempPotEnergyAve[j];
02271           adaptTempRead >> adaptTempPotEnergyVar[j];
02272           adaptTempRead >> adaptTempPotEnergySamples[j];
02273           adaptTempRead >> adaptTempPotEnergyAveNum[j];
02274           adaptTempRead >> adaptTempPotEnergyVarNum[j];
02275           adaptTempRead >> adaptTempPotEnergyAveDen[j];
02276           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02277         } 
02278         adaptTempRead.close();
02279       }
02280       else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
02281     } 
02282     else {
02283       adaptTempBins = simParams->adaptTempBins;
02284       adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
02285       adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
02286       adaptTempPotEnergySamples = new int[adaptTempBins];
02287       adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
02288       adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
02289       adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
02290       adaptTempBetaN           = new BigReal[adaptTempBins];
02291       adaptTempBetaMax = 1./simParams->adaptTempTmin;
02292       adaptTempBetaMin = 1./simParams->adaptTempTmax;
02293       adaptTempCg = simParams->adaptTempCgamma;   
02294       adaptTempDt  = simParams->adaptTempDt;
02295       adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
02296       adaptTempT = simParams->initialTemp; 
02297       if (simParams->langevinOn)
02298         adaptTempT = simParams->langevinTemp;
02299       else if (simParams->rescaleFreq > 0)
02300         adaptTempT = simParams->rescaleTemp;
02301       for(int j = 0; j < adaptTempBins; ++j){
02302           adaptTempPotEnergyAveNum[j] = 0.;
02303           adaptTempPotEnergyAveDen[j] = 0.;
02304           adaptTempPotEnergySamples[j] = 0;
02305           adaptTempPotEnergyVarNum[j] = 0.;
02306           adaptTempPotEnergyVar[j] = 0.;
02307           adaptTempPotEnergyAve[j] = 0.;
02308           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
02309       }
02310     }
02311     if (simParams->adaptTempAutoDt > 0.0) {
02312        adaptTempAutoDt = true;
02313        adaptTempDtMin =  simParams->adaptTempAutoDt - 0.01;
02314        adaptTempDtMax =  simParams->adaptTempAutoDt + 0.01;
02315     }
02316     adaptTempDTave = 0;
02317     adaptTempDTavenum = 0;
02318     iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
02319     iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
02320     iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
02321     iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
02322     iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
02323     if (simParams->adaptTempRestartFreq > 0) {
02324       NAMD_backup_file(simParams->adaptTempRestartFile);
02325       adaptTempRestartFile.open(simParams->adaptTempRestartFile);
02326        if(!adaptTempRestartFile)
02327         NAMD_die("Error opening restart file");
02328     }
02329 }

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

Definition at line 2357 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 2331 of file Controller.C.

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

02331                                                {
02332     if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
02333         adaptTempRestartFile.seekp(std::ios::beg);        
02334         iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
02335         adaptTempRestartFile << step << " ";
02336         // Start with min and max temperatures
02337         adaptTempRestartFile << adaptTempT<< " ";     // KELVIN
02338         adaptTempRestartFile << 1./adaptTempBetaMin << " ";  // KELVIN
02339         adaptTempRestartFile << 1./adaptTempBetaMax << " ";  // KELVIN
02340         adaptTempRestartFile << adaptTempBins << " ";     
02341         adaptTempRestartFile << adaptTempCg << " ";
02342         adaptTempRestartFile << adaptTempDt ;
02343         adaptTempRestartFile << "\n" ;
02344         for(int j = 0; j < adaptTempBins; ++j) {
02345           adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
02346           adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
02347           adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
02348           adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
02349           adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
02350           adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
02351           adaptTempRestartFile << "\n";          
02352         }
02353         adaptTempRestartFile.flush(); 
02354     }
02355 }    

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

Definition at line 274 of file Controller.C.

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

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

void Controller::awaken ( void   )  [inline]

Definition at line 45 of file Controller.h.

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

00045 { CthAwaken(thread); };

void Controller::berendsenPressure ( int   )  [protected]

Definition at line 938 of file Controller.C.

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

Referenced by integrate().

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

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

Referenced by rescaleaccelMD().

01662                                                               {
01663    switch(iE){
01664       case 2:
01665          *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
01666          // if k0 <=0 OR k0 > 1, switch to iE=1 mode
01667          if(*k0 > 1.)
01668             *warn |= 1;
01669          else if(*k0 <= 0.)
01670             *warn |= 2;
01671          //else stay in iE=2 mode
01672          else{
01673             *E = Vmin + (Vmax-Vmin)/(*k0);
01674             *iEused = 2;
01675             break;
01676          }
01677       case 1:
01678          *E = Vmax;
01679          *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
01680          if(*k0 > 1.)
01681             *k0 = 1.;
01682          *iEused = 1;
01683          break;
01684    }
01685 
01686    *k = *k0 / (Vmax - Vmin);
01687 }

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

Referenced by rescaleaccelMD().

01691                                            {
01692    BigReal VE = testV - E;
01693    //if V < E, apply boost
01694    if(VE < 0){
01695       *factor = k * VE;
01696       *vir += vir_orig * (*factor);
01697       *dV += (*factor) * VE/2.;
01698       *factor += 1.;
01699    }
01700    else{
01701       *factor = 1.;
01702    }
01703 }

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

Referenced by rescaleaccelMD().

01642                                                                            {
01643    BigReal delta; 
01644 
01645    if(testV > *Vmax){
01646       *Vmax = testV;
01647    }
01648    else if(testV < *Vmin){
01649       *Vmin = testV;
01650    }
01651 
01652    //mean and std calculated by Online Algorithm
01653    delta = testV - *Vavg;
01654    *Vavg += delta / (BigReal)n;
01655    *M2 += delta * (testV - (*Vavg));
01656 
01657    *sigmaV = sqrt(*M2/n);
01658 }

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

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

Referenced by multigratorPressure(), and receivePressure().

01526                                                                                             {
01527 
01528   Tensor virial_normal = virial_normal_in;
01529   Tensor virial_nbond = virial_nbond_in;
01530   Tensor virial_slow = virial_slow_in;
01531 
01532   // Tensor tmp = virial_normal;
01533   // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
01534   //   tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
01535 
01536   Node *node = Node::Object();
01537   Molecule *molecule = node->molecule;
01538   SimParameters *simParameters = node->simParameters;
01539   Lattice &lattice = state->lattice;
01540 
01541   BigReal volume;
01542 
01543   Vector extPosition = lattice.origin();
01544   virial_normal -= outer(extForce_normal,extPosition);
01545   virial_nbond -= outer(extForce_nbond,extPosition);
01546   virial_slow -= outer(extForce_slow,extPosition);
01547 
01548   if ( (volume=lattice.volume()) != 0. )
01549   {
01550 
01551     if (simParameters->LJcorrection && volume) {
01552 #ifdef MEM_OPT_VERSION
01553       NAMD_bug("LJcorrection not supported in memory optimized build.");
01554 #else
01555       // Apply tail correction to pressure
01556       BigReal alchLambda = simParameters->getCurrentLambda(step);
01557       virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
01558 #endif
01559     }
01560 
01561     // kinetic energy component included in virials
01562     pressure_normal = virial_normal / volume;
01563     groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01564 
01565     if (simParameters->accelMDOn) {
01566       pressure_amd = virial_amd / volume;
01567       pressure_normal += pressure_amd;
01568       groupPressure_normal +=  pressure_amd;
01569     }
01570 
01571     if ( minimize || ! ( step % nbondFreq ) )
01572     {
01573       pressure_nbond = virial_nbond / volume;
01574       groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01575     }
01576 
01577     if ( minimize || ! ( step % slowFreq ) )
01578     {
01579       pressure_slow = virial_slow / volume;
01580       groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01581     }
01582 
01583     pressure = pressure_normal + pressure_nbond + pressure_slow; 
01584     groupPressure = groupPressure_normal + groupPressure_nbond +
01585           groupPressure_slow;
01586   }
01587   else
01588   {
01589     pressure = Tensor();
01590     groupPressure = Tensor();
01591   }
01592 
01593   if ( simParameters->useGroupPressure )
01594   {
01595     controlPressure_normal = groupPressure_normal;
01596     controlPressure_nbond = groupPressure_nbond;
01597     controlPressure_slow = groupPressure_slow;
01598     controlPressure = groupPressure;
01599   }
01600   else
01601   {
01602     controlPressure_normal = pressure_normal;
01603     controlPressure_nbond = pressure_nbond;
01604     controlPressure_slow = pressure_slow;
01605     controlPressure = pressure;
01606   }
01607 
01608   if ( simParameters->useFlexibleCell ) {
01609     // use symmetric pressure to control rotation
01610     // controlPressure_normal = symmetric(controlPressure_normal);
01611     // controlPressure_nbond = symmetric(controlPressure_nbond);
01612     // controlPressure_slow = symmetric(controlPressure_slow);
01613     // controlPressure = symmetric(controlPressure);
01614     // only use on-diagonal components for now
01615     controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01616     controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01617     controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01618     controlPressure = Tensor::diagonal(diagonal(controlPressure));
01619     if ( simParameters->useConstantRatio ) {
01620 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01621    T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01622       AVGXY(controlPressure_normal);
01623       AVGXY(controlPressure_nbond);
01624       AVGXY(controlPressure_slow);
01625       AVGXY(controlPressure);
01626 #undef AVGXY
01627     }
01628   } else {
01629     controlPressure_normal =
01630   Tensor::identity(trace(controlPressure_normal)/3.);
01631     controlPressure_nbond =
01632   Tensor::identity(trace(controlPressure_nbond)/3.);
01633     controlPressure_slow =
01634   Tensor::identity(trace(controlPressure_slow)/3.);
01635     controlPressure =
01636   Tensor::identity(trace(controlPressure)/3.);
01637   }
01638 }

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

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

02655                                                          {
02656     Node *node = Node::Object();
02657     Molecule *molecule = node->molecule;
02658 
02659     // Some basic correctness checking
02660     BigReal checksum, checksum_b;
02661     char errmsg[256];
02662 
02663     checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
02664     if ( ((int)checksum) != molecule->numAtoms ) {
02665       sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
02666               (int)checksum, molecule->numAtoms);
02667       NAMD_bug(errmsg);
02668     }
02669 
02670     checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
02671     if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
02672       if ( ! computeChecksum ) {
02673         computesPartitioned = 0;
02674       } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
02675         computesPartitioned = 1;
02676       } else {
02677         NAMD_bug("Bad global compute count!\n");
02678       }
02679       computeChecksum = ((int)checksum);
02680     }
02681 
02682     checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
02683     if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
02684       sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
02685               (int)checksum, molecule->numCalcBonds);
02686       if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
02687         iout << iWARN << errmsg << endi;
02688       else NAMD_bug(errmsg);
02689     }
02690 
02691     checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
02692     if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
02693       sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
02694               (int)checksum, molecule->numCalcAngles);
02695       if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
02696         iout << iWARN << errmsg << endi;
02697       else NAMD_bug(errmsg);
02698     }
02699 
02700     checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
02701     if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
02702       sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
02703               (int)checksum, molecule->numCalcDihedrals);
02704       if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
02705         iout << iWARN << errmsg << endi;
02706       else NAMD_bug(errmsg);
02707     }
02708 
02709     checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
02710     if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
02711       sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
02712               (int)checksum, molecule->numCalcImpropers);
02713       if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
02714         iout << iWARN << errmsg << endi;
02715       else NAMD_bug(errmsg);
02716     }
02717 
02718     checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
02719     if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
02720       sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
02721               (int)checksum, molecule->numCalcTholes);
02722       if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
02723         iout << iWARN << errmsg << endi;
02724       else NAMD_bug(errmsg);
02725     }
02726 
02727     checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
02728     if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
02729       sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
02730               (int)checksum, molecule->numCalcAnisos);
02731       if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
02732         iout << iWARN << errmsg << endi;
02733       else NAMD_bug(errmsg);
02734     }
02735 
02736     checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
02737     if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
02738       sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
02739               (int)checksum, molecule->numCalcCrossterms);
02740       if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
02741         iout << iWARN << errmsg << endi;
02742       else NAMD_bug(errmsg);
02743     }
02744 
02745     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
02746     if ( ((int)checksum) > molecule->numCalcExclusions &&
02747          ( ! simParams->mollyOn || step % slowFreq ) ) {
02748       if ( forgiving )
02749         iout << iWARN << "High global exclusion count ("
02750                       << ((int)checksum) << " vs "
02751                       << molecule->numCalcExclusions << "), possible error!\n"
02752           << iWARN << "This warning is not unusual during minimization.\n"
02753           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02754       else {
02755         char errmsg[256];
02756         sprintf(errmsg, "High global exclusion count!  (%d vs %d)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02757                 (int)checksum, molecule->numCalcExclusions);
02758         NAMD_bug(errmsg);
02759       }
02760     }
02761 
02762     if ( ((int)checksum) && ((int)checksum) < molecule->numCalcExclusions &&
02763         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02764       if ( forgiving || ! simParams->fullElectFrequency ) {
02765         iout << iWARN << "Low global exclusion count!  ("
02766           << ((int)checksum) << " vs " << molecule->numCalcExclusions << ")";
02767         if ( forgiving ) {
02768           iout << "\n"
02769             << iWARN << "This warning is not unusual during minimization.\n"
02770             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02771         } else {
02772           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02773         }
02774         iout << endi;
02775       } else {
02776         char errmsg[256];
02777         sprintf(errmsg, "Low global exclusion count!  (%d vs %d)  System unstable or pairlistdist or cutoff too small.\n",
02778                 (int)checksum, molecule->numCalcExclusions);
02779         NAMD_bug(errmsg);
02780       }
02781     }
02782 
02783 #ifdef NAMD_CUDA
02784     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
02785     if ( ((int)checksum) > molecule->numCalcFullExclusions &&
02786          ( ! simParams->mollyOn || step % slowFreq ) ) {
02787       if ( forgiving )
02788         iout << iWARN << "High global CUDA exclusion count ("
02789                       << ((int)checksum) << " vs "
02790                       << molecule->numCalcFullExclusions << "), possible error!\n"
02791           << iWARN << "This warning is not unusual during minimization.\n"
02792           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
02793       else {
02794         char errmsg[256];
02795         sprintf(errmsg, "High global CUDA exclusion count!  (%d vs %d)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
02796                 (int)checksum, molecule->numCalcFullExclusions);
02797         NAMD_bug(errmsg);
02798       }
02799     }
02800 
02801     if ( ((int)checksum) && ((int)checksum) < molecule->numCalcFullExclusions &&
02802         ! simParams->goGroPair && ! simParams->qmForcesOn) {
02803       if ( forgiving || ! simParams->fullElectFrequency ) {
02804         iout << iWARN << "Low global CUDA exclusion count!  ("
02805           << ((int)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
02806         if ( forgiving ) {
02807           iout << "\n"
02808             << iWARN << "This warning is not unusual during minimization.\n"
02809             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
02810         } else {
02811           iout << "  System unstable or pairlistdist or cutoff too small.\n";
02812         }
02813         iout << endi;
02814       } else {
02815         char errmsg[256];
02816         sprintf(errmsg, "Low global CUDA exclusion count!  (%d vs %d)  System unstable or pairlistdist or cutoff too small.\n",
02817                 (int)checksum, molecule->numCalcFullExclusions);
02818         NAMD_bug(errmsg);
02819       }
02820     }
02821 #endif
02822 
02823     checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
02824     if ( ((int)checksum) && ! marginViolations ) {
02825       iout << iERROR << "Margin is too small for " << ((int)checksum) <<
02826         " atoms during timestep " << step << ".\n" << iERROR <<
02827       "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
02828     }
02829     marginViolations += (int)checksum;
02830 
02831     checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
02832     if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
02833       iout << iINFO <<
02834         "Pairlistdist is too small for " << ((int)checksum) <<
02835         " computes during timestep " << step << ".\n" << endi;
02836     }
02837     if ( simParams->outputPairlists )  pairlistWarnings += (int)checksum;
02838 
02839     checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
02840     if ( checksum ) {
02841       if ( forgiving )
02842         iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
02843       else NAMD_bug("Stray PME grid charges detected!\n");
02844     }
02845 }

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

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

03731                                                   {
03732   // alchemical scaling factors for groups 1/2 at the previous lambda
03733   const BigReal oldLambda = simParams->getCurrentLambda(step-1);
03734   const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
03735   const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
03736   const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
03737   const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
03738   const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
03739   const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
03740   // alchemical scaling factors for groups 1/2 at the new lambda
03741   const BigReal alchLambda = simParams->getCurrentLambda(step);
03742   const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
03743   const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
03744   const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
03745   const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
03746   const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
03747   const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda); 
03748 
03749   return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
03750           (elec_lambda_12 - elec_lambda_1)*
03751           (electEnergy_ti_1 + electEnergySlow_ti_1 +  electEnergyPME_ti_1) +
03752           (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
03753           (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
03754           (elec_lambda_22 - elec_lambda_2)*
03755           (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) + 
03756           (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
03757   );
03758 }

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

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

01243                                          {
01244 
01245     Vector momentum;
01246     momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
01247     momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
01248     momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
01249     const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
01250 
01251     if ( momentum.length2() == 0. )
01252       iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
01253     if ( mass == 0. )
01254       NAMD_die("Total mass is zero in Controller::correctMomentum");
01255 
01256     momentum *= (-1./mass);
01257 
01258     broadcast->momentumCorrection.publish(step+slowFreq,momentum);
01259 }

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

Definition at line 4031 of file Controller.C.

References broadcast.

Referenced by integrate().

04031                                                      {
04032 #if USE_BARRIER
04033         if (doBarrier) {
04034           broadcast->cycleBarrier.publish(step,1);
04035           CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
04036                   CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
04037         }
04038 #endif
04039 }

void Controller::enqueueCollections ( int   )  [protected]

Definition at line 3516 of file Controller.C.

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

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

03517 {
03518   if ( Output::coordinateNeeded(timestep) ){
03519     collection->enqueuePositions(timestep,state->lattice);
03520   }
03521   if ( Output::velocityNeeded(timestep) )
03522     collection->enqueueVelocities(timestep);
03523   if ( Output::forceNeeded(timestep) )
03524     collection->enqueueForces(timestep);
03525 }

void Controller::integrate ( int   )  [protected]

Definition at line 418 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, tcoupleVelocities(), traceBarrier(), SimParameters::traceStartStep, and SimParameters::zeroMomentum.

Referenced by algorithm().

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

void Controller::langevinPiston1 ( int   )  [protected]

Definition at line 975 of file Controller.C.

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

Referenced by integrate().

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

void Controller::langevinPiston2 ( int   )  [protected]

Definition at line 1130 of file Controller.C.

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

Referenced by integrate().

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

void Controller::minimize (  )  [protected]

Definition at line 582 of file Controller.C.

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

Referenced by algorithm().

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

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

Definition at line 907 of file Controller.C.

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

Referenced by printEnergies().

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

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

Definition at line 766 of file Controller.C.

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

Referenced by integrate().

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

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

Definition at line 841 of file Controller.C.

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

Referenced by integrate().

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

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

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

03860 {
03861 
03862   if ( step >= 0 ) {
03863 
03864     // Write out eXtended System Trajectory (XST) file
03865     if ( simParams->xstFrequency &&
03866          ((step % simParams->xstFrequency) == 0) )
03867     {
03868       if ( ! xstFile.is_open() )
03869       {
03870         iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
03871         NAMD_backup_file(simParams->xstFilename);
03872         xstFile.open(simParams->xstFilename);
03873         while (!xstFile) {
03874           if ( errno == EINTR ) {
03875             CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
03876             xstFile.clear();
03877             xstFile.open(simParams->xstFilename);
03878             continue;
03879           }
03880           char err_msg[257];
03881           sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
03882           NAMD_err(err_msg);
03883         }
03884         xstFile << "# NAMD extended system trajectory file" << std::endl;
03885         writeExtendedSystemLabels(xstFile);
03886       }
03887       writeExtendedSystemData(step,xstFile);
03888       xstFile.flush();
03889     }
03890 
03891     // Write out eXtended System Configuration (XSC) files
03892     //  Output a restart file
03893     if ( simParams->restartFrequency &&
03894          ((step % simParams->restartFrequency) == 0) &&
03895          (step != simParams->firstTimestep) )
03896     {
03897       iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
03898                 << step << "\n" << endi;
03899       char fname[140];
03900       const char *bsuffix = ".old";
03901       strcpy(fname, simParams->restartFilename);
03902       if ( simParams->restartSave ) {
03903         char timestepstr[20];
03904         sprintf(timestepstr,".%d",step);
03905         strcat(fname, timestepstr);
03906         bsuffix = ".BAK";
03907       }
03908       strcat(fname, ".xsc");
03909       NAMD_backup_file(fname,bsuffix);
03910       ofstream_namd xscFile(fname);
03911       while (!xscFile) {
03912         if ( errno == EINTR ) {
03913           CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
03914           xscFile.clear();
03915           xscFile.open(fname);
03916           continue;
03917         }
03918         char err_msg[257];
03919         sprintf(err_msg, "Error opening XSC restart file %s",fname);
03920         NAMD_err(err_msg);
03921       } 
03922       xscFile << "# NAMD extended system configuration restart file" << std::endl;
03923       writeExtendedSystemLabels(xscFile);
03924       writeExtendedSystemData(step,xscFile);
03925       if (!xscFile) {
03926         char err_msg[257];
03927         sprintf(err_msg, "Error writing XSC restart file %s",fname);
03928         NAMD_err(err_msg);
03929       } 
03930     }
03931 
03932   }
03933 
03934   //  Output final coordinates
03935   if (step == FILE_OUTPUT || step == END_OF_RUN)
03936   {
03937     int realstep = ( step == FILE_OUTPUT ?
03938         simParams->firstTimestep : simParams->N );
03939     iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
03940                 << realstep << "\n" << endi;
03941     static char fname[140];
03942     strcpy(fname, simParams->outputFilename);
03943     strcat(fname, ".xsc");
03944     NAMD_backup_file(fname);
03945     ofstream_namd xscFile(fname);
03946     while (!xscFile) {
03947       if ( errno == EINTR ) {
03948         CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
03949         xscFile.clear();
03950         xscFile.open(fname);
03951         continue;
03952       }
03953       char err_msg[257];
03954       sprintf(err_msg, "Error opening XSC output file %s",fname);
03955       NAMD_err(err_msg);
03956     } 
03957     xscFile << "# NAMD extended system configuration output file" << std::endl;
03958     writeExtendedSystemLabels(xscFile);
03959     writeExtendedSystemData(realstep,xscFile);
03960     if (!xscFile) {
03961       char err_msg[257];
03962       sprintf(err_msg, "Error writing XSC output file %s",fname);
03963       NAMD_err(err_msg);
03964     } 
03965   }
03966 
03967   //  Close trajectory file
03968   if (step == END_OF_RUN) {
03969     if ( xstFile.is_open() ) {
03970       xstFile.close();
03971       iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
03972     }
03973   }
03974 
03975 }

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

Definition at line 3528 of file Controller.C.

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

Referenced by integrate().

03528                                          {
03529  if (simParams->alchOn && simParams->alchFepOn) {
03530   const int stepInRun = step - simParams->firstTimestep;
03531   const int alchEquilSteps = simParams->alchEquilSteps;
03532   const BigReal alchLambda = simParams->alchLambda;
03533   const BigReal alchLambda2 = simParams->alchLambda2;
03534   const bool alchEnsembleAvg = simParams->alchEnsembleAvg; 
03535   const bool FepWhamOn = simParams->alchFepWhamOn;
03536 
03537   if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
03538     FepNo = 0;
03539     exp_dE_ByRT = 0.0;
03540     net_dE = 0.0;
03541   }
03542   dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
03543        electEnergy - electEnergySlow - ljEnergy;
03544   BigReal RT = BOLTZMANN * simParams->alchTemp;
03545 
03546   if (alchEnsembleAvg){
03547     FepNo++;
03548     exp_dE_ByRT += exp(-dE/RT);
03549     net_dE += dE;
03550   }
03551 
03552   if (simParams->alchOutFreq) { 
03553     if (stepInRun == 0) {
03554       if (!fepFile.is_open()) {
03555         NAMD_backup_file(simParams->alchOutFile);
03556         fepFile.open(simParams->alchOutFile);
03557         iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
03558         if(alchEnsembleAvg){
03559           fepSum = 0.0;
03560           fepFile << "#            STEP                 Elec                            "
03561                   << "vdW                    dE           dE_avg         Temp             dG\n"
03562                   << "#                           l             l+dl      "
03563                   << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03564         }
03565         else{
03566           if(!FepWhamOn){ 
03567             fepFile << "#            STEP                 Elec                            "
03568                     << "vdW                    dE         Temp\n"
03569                     << "#                           l             l+dl      "
03570                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
03571           } 
03572         }
03573       }
03574       if(!step){
03575         if(!FepWhamOn){ 
03576           fepFile << "#NEW FEP WINDOW: "
03577                   << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
03578                   << alchLambda2 << std::endl;
03579         }
03580       }
03581     }
03582     if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03583       fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03584               << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03585               << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03586     }
03587     if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
03588       writeFepEnergyData(step, fepFile);
03589       fepFile.flush();
03590     }
03591     if (alchEnsembleAvg && (step == simParams->N)) {
03592       fepSum = fepSum + dG;
03593       fepFile << "#Free energy change for lambda window [ " << alchLambda
03594               << " " << alchLambda2 << " ] is " << dG 
03595               << " ; net change until now is " << fepSum << std::endl;
03596       fepFile.flush();
03597     }
03598   }
03599  }
03600 }

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

Definition at line 3602 of file Controller.C.

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

Referenced by integrate().

03602                                         {
03603  if (simParams->alchOn && simParams->alchThermIntOn) {
03604   const int stepInRun = step - simParams->firstTimestep;
03605   const int alchEquilSteps = simParams->alchEquilSteps;
03606   const int alchLambdaFreq = simParams->alchLambdaFreq;
03607 
03608   if (stepInRun == 0 || stepInRun == alchEquilSteps) {
03609     TiNo = 0;
03610     net_dEdl_bond_1 = 0;
03611     net_dEdl_bond_2 = 0;
03612     net_dEdl_elec_1 = 0;
03613     net_dEdl_elec_2 = 0;
03614     net_dEdl_lj_1 = 0;
03615     net_dEdl_lj_2 = 0;
03616   }
03617   TiNo++;
03618   net_dEdl_bond_1 += bondedEnergy_ti_1;
03619   net_dEdl_bond_2 += bondedEnergy_ti_2;
03620   net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
03621                       + electEnergyPME_ti_1);
03622   net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
03623                       + electEnergyPME_ti_2);
03624   net_dEdl_lj_1 += ljEnergy_ti_1;
03625   net_dEdl_lj_2 += ljEnergy_ti_2;
03626 
03627   // Don't accumulate block averages (you'll get a divide by zero!) or even 
03628   // open the TI output if alchOutFreq is zero.
03629   if (simParams->alchOutFreq) {
03630     if (stepInRun == 0 || stepInRun == alchEquilSteps 
03631         || (! ((step - 1) % simParams->alchOutFreq))) {
03632       // output of instantaneous dU/dl now replaced with running average
03633       // over last alchOutFreq steps (except for step 0)
03634       recent_TiNo = 0;
03635       recent_dEdl_bond_1 = 0;
03636       recent_dEdl_bond_2 = 0;
03637       recent_dEdl_elec_1 = 0;
03638       recent_dEdl_elec_2 = 0;
03639       recent_dEdl_lj_1 = 0;
03640       recent_dEdl_lj_2 = 0;
03641       recent_alchWork = 0;
03642     }
03643     recent_TiNo++;
03644     recent_dEdl_bond_1 += bondedEnergy_ti_1;
03645     recent_dEdl_bond_2 += bondedEnergy_ti_2;
03646     recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 
03647                            + electEnergyPME_ti_1); 
03648     recent_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 
03649                            + electEnergyPME_ti_2); 
03650     recent_dEdl_lj_1 += ljEnergy_ti_1;
03651     recent_dEdl_lj_2 += ljEnergy_ti_2;
03652     recent_alchWork += alchWork;
03653 
03654     if (stepInRun == 0) {
03655       if (!tiFile.is_open()) {
03656         NAMD_backup_file(simParams->alchOutFile);
03657         tiFile.open(simParams->alchOutFile);
03658         /* BKR - This has been rather drastically updated to better match 
03659            stdout. This was necessary for several reasons:
03660            1) PME global scaling is obsolete (now removed)
03661            2) scaling of bonded terms was added
03662            3) alchemical work is now accumulated when switching is active
03663          */
03664         iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
03665         tiFile << "#TITITLE:    TS";
03666         tiFile << FORMAT("BOND1");
03667         tiFile << FORMAT("AVGBOND1");
03668         tiFile << FORMAT("ELECT1");
03669         tiFile << FORMAT("AVGELECT1");
03670         tiFile << "     ";
03671         tiFile << FORMAT("VDW1");
03672         tiFile << FORMAT("AVGVDW1");
03673         tiFile << FORMAT("BOND2");
03674         tiFile << FORMAT("AVGBOND2");
03675         tiFile << FORMAT("ELECT2");
03676         tiFile << "     ";
03677         tiFile << FORMAT("AVGELECT2");
03678         tiFile << FORMAT("VDW2");
03679         tiFile << FORMAT("AVGVDW2");
03680         if (alchLambdaFreq > 0) {
03681           tiFile << FORMAT("ALCHWORK");
03682           tiFile << FORMAT("CUMALCHWORK");
03683         }
03684         tiFile << std::endl;
03685       }
03686 
03687       if (alchLambdaFreq > 0) {
03688         tiFile << "#ALCHEMICAL SWITCHING ACTIVE " 
03689                << simParams->alchLambda << " --> " << simParams->alchLambda2
03690                << "\n#LAMBDA SCHEDULE: " 
03691                << "dL: " << simParams->getLambdaDelta() 
03692                << " Freq: " << alchLambdaFreq;
03693       }
03694       else {
03695         const BigReal alchLambda = simParams->alchLambda;    
03696         const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
03697         const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
03698         const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
03699         const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
03700         const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
03701         const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
03702         tiFile << "#NEW TI WINDOW: "
03703                << "LAMBDA " << alchLambda 
03704                << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
03705                << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
03706                << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
03707                << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
03708       }
03709       tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
03710              << std::endl;
03711     }
03712 
03713     if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
03714       tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
03715              << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
03716              << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
03717     }
03718     if ((step%simParams->alchOutFreq) == 0) {
03719       writeTiEnergyData(step, tiFile);
03720       tiFile.flush();
03721     }
03722   }
03723  }
03724 }

void Controller::printDynamicsEnergies ( int   )  [protected]

Definition at line 2903 of file Controller.C.

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

Referenced by integrate().

02903                                                {
02904 
02905     Node *node = Node::Object();
02906     Molecule *molecule = node->molecule;
02907     SimParameters *simParameters = node->simParameters;
02908     Lattice &lattice = state->lattice;
02909 
02910     compareChecksums(step);
02911 
02912     printEnergies(step,0);
02913 }

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

Definition at line 2915 of file Controller.C.

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

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

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

void Controller::printFepMessage ( int   )  [protected]

Definition at line 1262 of file Controller.C.

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

Referenced by integrate().

01263 {
01264   if (simParams->alchOn && simParams->alchFepOn 
01265       && !simParams->alchLambdaFreq) {
01266     const BigReal alchLambda = simParams->alchLambda;
01267     const BigReal alchLambda2 = simParams->alchLambda2;
01268     const BigReal alchTemp = simParams->alchTemp;
01269     const int alchEquilSteps = simParams->alchEquilSteps;
01270     iout << "FEP: RESETTING FOR NEW FEP WINDOW "
01271          << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2
01272          << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
01273          << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
01274          << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp 
01275          << " K FOR FEP CALCULATION\n" << endi;
01276   }
01277 } 

void Controller::printMinimizeEnergies ( int   )  [protected]

Definition at line 2884 of file Controller.C.

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

02884                                                {
02885 
02886     rescaleaccelMD(step,1);
02887     receivePressure(step,1);
02888 
02889     Node *node = Node::Object();
02890     Molecule *molecule = node->molecule;
02891     compareChecksums(step,1);
02892 
02893     printEnergies(step,1);
02894 
02895     min_energy = totalEnergy;
02896     min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
02897     min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
02898     min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
02899     min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
02900 
02901 }

void Controller::printTiMessage ( int   )  [protected]

Definition at line 1278 of file Controller.C.

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

Referenced by integrate().

01279 {
01280   if (simParams->alchOn && simParams->alchThermIntOn 
01281       && !simParams->alchLambdaFreq) {
01282     const BigReal alchLambda = simParams->alchLambda;
01283     const int alchEquilSteps = simParams->alchEquilSteps;
01284     iout << "TI: RESETTING FOR NEW WINDOW "
01285          << "LAMBDA SET TO " << alchLambda 
01286          << "\nTI: WINDOW TO HAVE " << alchEquilSteps 
01287          << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
01288   }
01289 } 

void Controller::printTiming ( int   )  [protected]

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

02847                                      {
02848 
02849     if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
02850     {
02851       const double endWTime = CmiWallTimer() - firstWTime;
02852       const double endCTime = CmiTimer() - firstCTime;
02853 
02854       // fflush about once per minute
02855       if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
02856 
02857       const double elapsedW = 
02858         (endWTime - startWTime) / simParams->outputTiming;
02859       const double elapsedC = 
02860         (endCTime - startCTime) / simParams->outputTiming;
02861 
02862       const double remainingW = elapsedW * (simParams->N - step);
02863       const double remainingW_hours = remainingW / 3600;
02864 
02865       startWTime = endWTime;
02866       startCTime = endCTime;
02867 
02868 #ifdef NAMD_CUDA
02869       if ( simParams->outputEnergies < 60 &&
02870            step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
02871         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
02872       }
02873 #endif
02874       if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
02875         CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
02876                   ", %g hours remaining, %f MB of memory in use.\n",
02877                   step, endCTime, elapsedC, endWTime, elapsedW,
02878                   remainingW_hours, memusage_MB());
02879         if ( fflush_count ) { --fflush_count; fflush(stdout); }
02880       }
02881     }
02882 }

void Controller::reassignVelocities ( int   )  [protected]

Definition at line 1292 of file Controller.C.

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

Referenced by integrate().

01293 {
01294   const int reassignFreq = simParams->reassignFreq;
01295   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
01296     BigReal newTemp = simParams->reassignTemp;
01297     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
01298     if ( simParams->reassignIncr > 0.0 ) {
01299       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
01300         newTemp = simParams->reassignHold;
01301     } else {
01302       if ( newTemp < simParams->reassignHold )
01303         newTemp = simParams->reassignHold;
01304     }
01305     iout << "REASSIGNING VELOCITIES AT STEP " << step
01306          << " TO " << newTemp << " KELVIN.\n" << endi;
01307   }
01308 }

void Controller::rebalanceLoad ( int   )  [protected]

Definition at line 4017 of file Controller.C.

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

Referenced by integrate().

04018 {
04019   if ( ! ldbSteps ) { 
04020     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
04021   }
04022   if ( ! --ldbSteps ) {
04023     startBenchTime -= CmiWallTimer();
04024         Node::Object()->outputPatchComputeMaps("before_ldb", step);
04025     LdbCoordinator::Object()->rebalance(this);  
04026         startBenchTime += CmiWallTimer();
04027     fflush_count = 3;
04028   }
04029 }

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

Definition at line 1345 of file Controller.C.

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

Referenced by integrate(), and printMinimizeEnergies().

01346 {
01347     Node *node = Node::Object();
01348     Molecule *molecule = node->molecule;
01349     SimParameters *simParameters = node->simParameters;
01350     Lattice &lattice = state->lattice;
01351 
01352     reduction->require();
01353 
01354     Tensor virial_normal;
01355     Tensor virial_nbond;
01356     Tensor virial_slow;
01357 #ifdef ALTVIRIAL
01358     Tensor altVirial_normal;
01359     Tensor altVirial_nbond;
01360     Tensor altVirial_slow;
01361 #endif
01362     Tensor intVirial_normal;
01363     Tensor intVirial_nbond;
01364     Tensor intVirial_slow;
01365     Vector extForce_normal;
01366     Vector extForce_nbond;
01367     Vector extForce_slow;
01368 
01369 #if 1
01370     numDegFreedom = molecule->num_deg_freedom();
01371     int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
01372     int numFixedGroups = molecule->num_fixed_groups();
01373     int numFixedAtoms = molecule->num_fixed_atoms();
01374 #endif
01375 #if 0
01376     int numAtoms = molecule->numAtoms;
01377     numDegFreedom = 3 * numAtoms;
01378     int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
01379     int numFixedAtoms =
01380         ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
01381     int numLonepairs = molecule->numLonepairs;
01382     int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
01383     if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
01384     if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01385     if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01386     if ( ! ( numFixedAtoms || molecule->numConstraints
01387         || simParameters->comMove || simParameters->langevinOn ) ) {
01388       numDegFreedom -= 3;
01389       numGroupDegFreedom -= 3;
01390     }
01391     if (simParameters->pairInteractionOn) {
01392       // this doesn't attempt to deal with fixed atoms or constraints
01393       numDegFreedom = 3 * molecule->numFepInitial;
01394     }
01395     int numRigidBonds = molecule->numRigidBonds;
01396     int numFixedRigidBonds =
01397         ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01398 
01399     // numLonepairs is subtracted here because all lonepairs have a rigid bond
01400     // to oxygen, but all of the LP degrees of freedom are dealt with above
01401     numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01402 #endif
01403 
01404     kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01405     kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01406 
01407     BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01408         reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01409     BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01410         reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01411 
01412     BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01413                                         / ( numDegFreedom * BOLTZMANN );
01414     BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01415                                         / ( numDegFreedom * BOLTZMANN );
01416     BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01417                                         / ( numGroupDegFreedom * BOLTZMANN );
01418     BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01419                                         / ( numGroupDegFreedom * BOLTZMANN );
01420 
01421     /*  test code for comparing different temperatures  
01422     iout << "TEMPTEST: " << step << " " << 
01423         atomTempHalfstep << " " <<
01424         atomTempCentered << " " <<
01425         groupTempHalfstep << " " <<
01426         groupTempCentered << "\n" << endi;
01427   iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
01428     numGroupDegFreedom << "\n" << endi;
01429      */
01430 
01431     GET_TENSOR(momentumSqrSum,reduction,REDUCTION_MOMENTUM_SQUARED);
01432 
01433     GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01434     GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01435     GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01436 
01437 #ifdef ALTVIRIAL
01438     GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01439     GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01440     GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01441 #endif
01442 
01443     GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01444     GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01445     GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01446 
01447     GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01448     GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01449     GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01450     // APH NOTE: These four lines are now done in calcPressure()
01451     // Vector extPosition = lattice.origin();
01452     // virial_normal -= outer(extForce_normal,extPosition);
01453     // virial_nbond -= outer(extForce_nbond,extPosition);
01454     // virial_slow -= outer(extForce_slow,extPosition);
01455 
01456     kineticEnergy = kineticEnergyCentered;
01457     temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );
01458 
01459     if (simParameters->drudeOn) {
01460       BigReal drudeComKE
01461         = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01462       BigReal drudeBondKE
01463         = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01464       int g_bond = 3 * molecule->numDrudeAtoms;
01465       int g_com = numDegFreedom - g_bond;
01466       temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
01467       drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
01468     }
01469 
01470     // Calculate number of degrees of freedom (controlNumDegFreedom)
01471     if ( simParameters->useGroupPressure )
01472     {
01473       controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01474       if ( ! ( numFixedAtoms || molecule->numConstraints
01475   || simParameters->comMove || simParameters->langevinOn ) ) {
01476         controlNumDegFreedom -= 1;
01477       }
01478     }
01479     else
01480     {
01481       controlNumDegFreedom = numDegFreedom / 3;
01482     }
01483     if (simParameters->fixCellDims) {
01484       if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01485       if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01486       if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01487     }
01488 
01489     // Calculate pressure tensors using virials
01490     calcPressure(step, minimize,
01491       virial_normal, virial_nbond, virial_slow,
01492       intVirial_normal, intVirial_nbond, intVirial_slow,
01493       extForce_normal, extForce_nbond, extForce_slow);
01494 
01495 #ifdef DEBUG_PRESSURE
01496     iout << iINFO << "Control pressure = " << controlPressure <<
01497       " = " << controlPressure_normal << " + " <<
01498       controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01499 #endif
01500     if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
01501         iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
01502              << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
01503              << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
01504              << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
01505              << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
01506              << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
01507              << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
01508              << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
01509              << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
01510              << endi;
01511    }
01512 }

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

Definition at line 4008 of file Controller.C.

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

Referenced by Node::recvCheckpointAck().

04008                                                  {  // initiating replica
04009   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
04010     state->lattice = cp.lattice;
04011     *(ControllerState*)this = cp.state;
04012   }
04013   CkpvAccess(_qd)->process();
04014 }

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

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

03978                                                                             {  // responding replica
03979   switch ( task ) {
03980     case SCRIPT_CHECKPOINT_STORE:
03981       if ( ! checkpoints.count(key) ) {
03982         checkpoints[key] = new checkpoint;
03983       }
03984       *checkpoints[key] = cp;
03985       break;
03986     case SCRIPT_CHECKPOINT_LOAD:
03987       if ( ! checkpoints.count(key) ) {
03988         NAMD_die("Unable to load checkpoint, requested key was never stored.");
03989       }
03990       cp = *checkpoints[key];
03991       break;
03992     case SCRIPT_CHECKPOINT_SWAP:
03993       if ( ! checkpoints.count(key) ) {
03994         NAMD_die("Unable to swap checkpoint, requested key was never stored.");
03995       }
03996       std::swap(cp,*checkpoints[key]);
03997       break;
03998     case SCRIPT_CHECKPOINT_FREE:
03999       if ( ! checkpoints.count(key) ) {
04000         NAMD_die("Unable to free checkpoint, requested key was never stored.");
04001       }
04002       delete checkpoints[key];
04003       checkpoints.erase(key);
04004       break;
04005   }
04006 }

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

Definition at line 1764 of file Controller.C.

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

Referenced by integrate(), and printMinimizeEnergies().

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

void Controller::rescaleVelocities ( int   )  [protected]

Definition at line 1221 of file Controller.C.

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

Referenced by integrate().

01222 {
01223   const int rescaleFreq = simParams->rescaleFreq;
01224   if ( rescaleFreq > 0 ) {
01225     rescaleVelocities_sumTemps += temperature;  ++rescaleVelocities_numTemps;
01226     if ( rescaleVelocities_numTemps == rescaleFreq ) {
01227       BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
01228       BigReal rescaleTemp = simParams->rescaleTemp;
01229       if ( simParams->adaptTempOn && simParams->adaptTempRescale && (step > simParams->adaptTempFirstStep ) && 
01230                 (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
01231         rescaleTemp = adaptTempT;
01232       }
01233       BigReal factor = sqrt(rescaleTemp/avgTemp);
01234       broadcast->velocityRescaleFactor.publish(step,factor);
01235       //iout << "RESCALING VELOCITIES AT STEP " << step
01236       //     << " FROM AVERAGE TEMPERATURE OF " << avgTemp
01237       //     << " TO " << rescaleTemp << " KELVIN.\n" << endi;
01238       rescaleVelocities_sumTemps = 0;  rescaleVelocities_numTemps = 0;
01239     }
01240   }
01241 }

void Controller::resumeAfterTraceBarrier ( int   ) 

Definition at line 4048 of file Controller.C.

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

Referenced by Node::resumeAfterTraceBarrier().

04048                                                 {
04049         broadcast->traceBarrier.publish(step,1);
04050         CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime); 
04051         awaken();
04052 }

void Controller::run ( void   ) 

Definition at line 261 of file Controller.C.

References awaken(), CTRL_STK_SZ, and DebugM.

Referenced by NamdState::runController().

00262 {
00263     // create a Thread and invoke it
00264     DebugM(4, "Starting thread in controller on this=" << this << "\n");
00265     thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
00266     CthSetStrategyDefault(thread);
00267 #if CMK_BLUEGENE_CHARM
00268     BgAttach(thread);
00269 #endif
00270     awaken();
00271 }

void Controller::tcoupleVelocities ( int   )  [protected]

Definition at line 1310 of file Controller.C.

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

Referenced by integrate().

01311 {
01312   if ( simParams->tCoupleOn )
01313   {
01314     const BigReal tCoupleTemp = simParams->tCoupleTemp;
01315     BigReal coefficient = 1.;
01316     if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
01317     broadcast->tcoupleCoefficient.publish(step,coefficient);
01318   }
01319 }

void Controller::terminate ( void   )  [protected]

Definition at line 4069 of file Controller.C.

References BackEnd::awaken().

Referenced by algorithm().

04069                                {
04070   BackEnd::awaken();
04071   CthFree(thread);
04072   CthSuspend();
04073 }

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

Definition at line 4041 of file Controller.C.

Referenced by integrate().

04041                                                        {
04042         CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);       
04043         CProxy_Node nd(CkpvAccess(BOCclass_group).node);
04044         nd.traceBarrier(turnOnTrace, step);
04045         CthSuspend();
04046 }

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

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

Referenced by rescaleaccelMD().

01706                                                                                                                                                              {
01707    FILE *rest;
01708    char timestepstr[20];
01709    char *filename;
01710    int baselen;
01711    char *restart_name;
01712    const char *bsuffix;
01713 
01714    if(lasttime || simParams->restartFrequency == 0){
01715            filename = simParams->outputFilename;
01716            bsuffix = ".BAK";
01717    }
01718    else{
01719            filename = simParams->restartFilename;
01720            bsuffix = ".old";
01721    }
01722 
01723    baselen = strlen(filename);
01724    restart_name = new char[baselen+26];
01725 
01726    strcpy(restart_name, filename);
01727    if ( simParams->restartSave && !lasttime) {
01728       sprintf(timestepstr,".%d",step_n);
01729       strcat(restart_name, timestepstr);
01730       bsuffix = ".BAK";
01731    }
01732    strcat(restart_name, ".gamd");
01733 
01734    if(write_topic){
01735       NAMD_backup_file(restart_name,bsuffix);
01736 
01737       rest = fopen(restart_name, "w");
01738       if(rest){
01739          fprintf(rest, "# NAMD accelMDG restart file\n"
01740                "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
01741                "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01742                type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01743          fclose(rest);
01744          iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
01745       }
01746       else
01747          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01748    }
01749    else{
01750       rest = fopen(restart_name, "a");
01751       if(rest){
01752          fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
01753                           type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
01754          fclose(rest);
01755       }
01756       else
01757          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
01758    }
01759 
01760    delete[] restart_name;
01761 }

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

Definition at line 3495 of file Controller.C.

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

Referenced by outputExtendedSystem().

03495                                                                       {
03496   Lattice &lattice = state->lattice;
03497   file.precision(12);
03498   file << step;
03499     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
03500     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
03501     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
03502     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
03503   if ( simParams->langevinPistonOn ) {
03504     Vector strainRate = diagonal(langevinPiston_strainRate);
03505     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
03506     file << " " << strainRate.x;
03507     file << " " << strainRate.y;
03508     file << " " << strainRate.z;
03509     file << " " << strainRate2.x;
03510     file << " " << strainRate2.y;
03511     file << " " << strainRate2.z;
03512   }
03513   file << std::endl;
03514 }

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

Definition at line 3482 of file Controller.C.

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

Referenced by outputExtendedSystem().

03482                                                               {
03483   Lattice &lattice = state->lattice;
03484   file << "#$LABELS step";
03485   if ( lattice.a_p() ) file << " a_x a_y a_z";
03486   if ( lattice.b_p() ) file << " b_x b_y b_z";
03487   if ( lattice.c_p() ) file << " c_x c_y c_z";
03488   file << " o_x o_y o_z";
03489   if ( simParams->langevinPistonOn ) {
03490     file << " s_x s_y s_z s_u s_v s_w";
03491   }
03492   file << std::endl;
03493 }

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

Definition at line 3760 of file Controller.C.

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

Referenced by outputFepEnergy().

03760                                                                  {
03761   BigReal eeng = electEnergy + electEnergySlow;
03762   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
03763   BigReal dE_Left = eeng_f + ljEnergy_f_left - eeng - ljEnergy;
03764   BigReal RT = BOLTZMANN * simParams->alchTemp;
03765 
03766   const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
03767   const int stepInRun = step - simParams->firstTimestep;
03768   const bool FepWhamOn = simParams->alchFepWhamOn;
03769   const bool WCARepuOn = simParams->alchFepWCARepuOn;
03770   const BigReal WCArcut1 = simParams->alchFepWCArcut1;
03771   const BigReal WCArcut2 = simParams->alchFepWCArcut2;
03772   const BigReal WCArcut3 = simParams->alchFepWCArcut3;
03773   const BigReal alchLambda = simParams->alchLambda;
03774         
03775   const BigReal alchRepLambda = simParams->alchRepLambda;
03776   const BigReal alchDispLambda = simParams->alchDispLambda;
03777   const BigReal alchElecLambda = simParams->alchElecLambda;
03778 
03779   if(stepInRun){
03780     if(!FepWhamOn){
03781       fepFile << FEPTITLE(step);
03782       fepFile << FORMAT(eeng);
03783       fepFile << FORMAT(eeng_f);
03784       fepFile << FORMAT(ljEnergy);
03785       fepFile << FORMAT(ljEnergy_f);
03786     }
03787     else{ // FepWhamOn = ON
03788       if(WCARepuOn){
03789         if(WCArcut1<WCArcut2) { // [s1,s2]
03790           fepFile << "FEP_WCA_REP  ";
03791           fepFile << FORMAT(WCArcut1);
03792           fepFile << FORMAT(WCArcut2);
03793           fepFile << FORMAT(1.0);
03794           fepFile << FORMAT(dE_Left);
03795         }
03796         if(WCArcut2<WCArcut3) { // [s2,s3]
03797           if(WCArcut1<WCArcut2) fepFile << " BREAK ";
03798           fepFile << "FEP_WCA_REP  ";
03799           fepFile << FORMAT(WCArcut2);
03800           fepFile << FORMAT(WCArcut3);
03801           fepFile << FORMAT(0.0);
03802           fepFile << FORMAT(dE);
03803         }
03804         fepFile << std::endl;
03805       }
03806       else if(simParams->alchFepWCADispOn) {
03807         fepFile << "FEP_WCA_DISP ";
03808         fepFile << FORMAT(alchDispLambda);
03809       }
03810       else if(simParams->alchFepElecOn) {
03811         fepFile << "FEP_ELEC     ";
03812         fepFile << FORMAT(alchElecLambda);
03813       }
03814     }
03815     if( ! WCARepuOn ) {
03816       fepFile << FORMAT(dE);
03817     }
03818     if(alchEnsembleAvg){
03819       BigReal dE_avg = net_dE/FepNo;
03820       fepFile << FORMAT(dE_avg);
03821     }
03822     if(!FepWhamOn){
03823       fepFile << FORMAT(temperature);
03824     }
03825     if(alchEnsembleAvg){
03826       dG = -(RT * log(exp_dE_ByRT/FepNo));
03827       fepFile << FORMAT(dG);
03828     } 
03829     if( ! WCARepuOn ) {
03830       fepFile << std::endl;
03831     }
03832   }
03833 }

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

Definition at line 3835 of file Controller.C.

References SimParameters::alchLambdaFreq, cumAlchWork, FORMAT(), net_dEdl_bond_1, net_dEdl_bond_2, net_dEdl_elec_1, net_dEdl_elec_2, net_dEdl_lj_1, net_dEdl_lj_2, recent_alchWork, recent_dEdl_bond_1, recent_dEdl_bond_2, recent_dEdl_elec_1, recent_dEdl_elec_2, recent_dEdl_lj_1, recent_dEdl_lj_2, recent_TiNo, simParams, tiFile, TiNo, and TITITLE().

Referenced by outputTiEnergy().


Friends And Related Function Documentation

friend class CheckpointMsg [friend]

Definition at line 54 of file Controller.h.

friend class Node [friend]

Definition at line 53 of file Controller.h.

friend class ScriptTcl [friend]

Definition at line 52 of file Controller.h.


Member Data Documentation

BigReal Controller::accelMDdVAverage [protected]

Definition at line 276 of file Controller.h.

Referenced by rescaleaccelMD().

Bool Controller::adaptTempAutoDt [protected]

Definition at line 299 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempBetaMax [protected]

Definition at line 293 of file Controller.h.

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

BigReal Controller::adaptTempBetaMin [protected]

Definition at line 292 of file Controller.h.

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

BigReal* Controller::adaptTempBetaN [protected]

Definition at line 288 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

int Controller::adaptTempBin [protected]

Definition at line 294 of file Controller.h.

Referenced by adaptTempUpdate().

int Controller::adaptTempBins [protected]

Definition at line 295 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempCg [protected]

Definition at line 297 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempDBeta [protected]

Definition at line 296 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal Controller::adaptTempDt [protected]

Definition at line 298 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempDTave [protected]

Definition at line 290 of file Controller.h.

BigReal Controller::adaptTempDTavenum [protected]

Definition at line 291 of file Controller.h.

BigReal Controller::adaptTempDtMax [protected]

Definition at line 301 of file Controller.h.

Referenced by adaptTempInit().

BigReal Controller::adaptTempDtMin [protected]

Definition at line 300 of file Controller.h.

Referenced by adaptTempInit().

BigReal* Controller::adaptTempPotEnergyAve [protected]

Definition at line 285 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyAveDen [protected]

Definition at line 283 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyAveNum [protected]

Definition at line 282 of file Controller.h.

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

int* Controller::adaptTempPotEnergySamples [protected]

Definition at line 287 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyVar [protected]

Definition at line 286 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyVarNum [protected]

Definition at line 284 of file Controller.h.

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

ofstream_namd Controller::adaptTempRestartFile [protected]

Definition at line 302 of file Controller.h.

Referenced by adaptTempWriteRestart().

BigReal Controller::adaptTempT [protected]

Definition at line 289 of file Controller.h.

Referenced by adaptTempInit(), and rescaleVelocities().

BigReal Controller::alchWork [protected]

Definition at line 150 of file Controller.h.

Referenced by outputTiEnergy(), and printEnergies().

RequireReduction* Controller::amd_reduction [protected]

Definition at line 213 of file Controller.h.

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

int Controller::avg_count [protected]

Definition at line 83 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::bondedEnergy_ti_1 [protected]

Definition at line 125 of file Controller.h.

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

BigReal Controller::bondedEnergy_ti_2 [protected]

Definition at line 126 of file Controller.h.

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

BigReal Controller::bondedEnergyDiff_f [protected]

Definition at line 112 of file Controller.h.

Referenced by outputFepEnergy(), and printEnergies().

ControllerBroadcasts* Controller::broadcast [protected]

Definition at line 226 of file Controller.h.

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

Lattice Controller::checkpoint_lattice [protected]

Definition at line 244 of file Controller.h.

Referenced by algorithm().

ControllerState Controller::checkpoint_state [protected]

Definition at line 245 of file Controller.h.

Referenced by algorithm().

int Controller::checkpoint_stored [protected]

Definition at line 243 of file Controller.h.

Referenced by algorithm(), and Controller().

int Controller::checkpoint_task [protected]

Definition at line 252 of file Controller.h.

Referenced by recvCheckpointAck().

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

Definition at line 251 of file Controller.h.

Referenced by algorithm(), and recvCheckpointReq().

CollectionMaster* const Controller::collection [protected]

Definition at line 224 of file Controller.h.

Referenced by enqueueCollections().

int Controller::computeChecksum [protected]

Definition at line 88 of file Controller.h.

Referenced by compareChecksums().

int Controller::controlNumDegFreedom [protected]

Definition at line 165 of file Controller.h.

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

Tensor Controller::controlPressure [protected]

Definition at line 166 of file Controller.h.

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

Tensor Controller::controlPressure_nbond [protected]

Definition at line 76 of file Controller.h.

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

Tensor Controller::controlPressure_normal [protected]

Definition at line 75 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::controlPressure_slow [protected]

Definition at line 77 of file Controller.h.

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

BigReal Controller::cumAlchWork [protected]

Definition at line 139 of file Controller.h.

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

BigReal Controller::dE [protected]

Definition at line 118 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::dG [protected]

Definition at line 120 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::drudeBondTemp [protected]

Definition at line 154 of file Controller.h.

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

BigReal Controller::drudeBondTempAvg [protected]

Definition at line 155 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::electEnergy [protected]

Definition at line 103 of file Controller.h.

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

BigReal Controller::electEnergy_f [protected]

Definition at line 113 of file Controller.h.

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

BigReal Controller::electEnergy_ti_1 [protected]

Definition at line 127 of file Controller.h.

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

BigReal Controller::electEnergy_ti_2 [protected]

Definition at line 130 of file Controller.h.

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

BigReal Controller::electEnergyPME_ti_1 [protected]

Definition at line 140 of file Controller.h.

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

BigReal Controller::electEnergyPME_ti_2 [protected]

Definition at line 141 of file Controller.h.

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

BigReal Controller::electEnergySlow [protected]

Definition at line 104 of file Controller.h.

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

BigReal Controller::electEnergySlow_f [protected]

Definition at line 114 of file Controller.h.

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

BigReal Controller::electEnergySlow_ti_1 [protected]

Definition at line 128 of file Controller.h.

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

BigReal Controller::electEnergySlow_ti_2 [protected]

Definition at line 131 of file Controller.h.

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

BigReal Controller::exp_dE_ByRT [protected]

Definition at line 117 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

ofstream_namd Controller::fepFile [protected]

Definition at line 233 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

int Controller::FepNo [protected]

Definition at line 121 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::fepSum [protected]

Definition at line 123 of file Controller.h.

Referenced by outputFepEnergy().

int Controller::fflush_count [protected]

Definition at line 197 of file Controller.h.

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

BigReal Controller::goNativeEnergy [protected]

Definition at line 108 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::goNonnativeEnergy [protected]

Definition at line 109 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::goTotalEnergy [protected]

Definition at line 110 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::groGaussEnergy [protected]

Definition at line 107 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::groLJEnergy [protected]

Definition at line 106 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Tensor Controller::groupPressure [protected]

Definition at line 164 of file Controller.h.

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

BigReal Controller::groupPressure_avg [protected]

Definition at line 82 of file Controller.h.

Referenced by Controller(), and printEnergies().

Tensor Controller::groupPressure_nbond [protected]

Definition at line 73 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_normal [protected]

Definition at line 72 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_slow [protected]

Definition at line 74 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_tavg [protected]

Definition at line 85 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::kineticEnergy [protected]

Definition at line 157 of file Controller.h.

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

BigReal Controller::kineticEnergyCentered [protected]

Definition at line 159 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

BigReal Controller::kineticEnergyHalfstep [protected]

Definition at line 158 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Tensor Controller::langevinPiston_origStrainRate [protected]

Definition at line 179 of file Controller.h.

Referenced by Controller().

int Controller::ldbSteps [protected]

Definition at line 195 of file Controller.h.

Referenced by rebalanceLoad().

BigReal Controller::ljEnergy [protected]

Definition at line 105 of file Controller.h.

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

BigReal Controller::ljEnergy_f [protected]

Definition at line 115 of file Controller.h.

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

BigReal Controller::ljEnergy_f_left [protected]

Definition at line 116 of file Controller.h.

Referenced by printEnergies(), and writeFepEnergyData().

BigReal Controller::ljEnergy_ti_1 [protected]

Definition at line 129 of file Controller.h.

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

BigReal Controller::ljEnergy_ti_2 [protected]

Definition at line 132 of file Controller.h.

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

int Controller::marginViolations [protected]

Definition at line 89 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

BigReal Controller::min_energy [protected]

Definition at line 93 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

BigReal Controller::min_f_dot_f [protected]

Definition at line 94 of file Controller.h.

Referenced by minimize().

BigReal Controller::min_f_dot_v [protected]

Definition at line 95 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

int Controller::min_huge_count [protected]

Definition at line 97 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

RequireReduction* Controller::min_reduction [protected]

Definition at line 59 of file Controller.h.

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

BigReal Controller::min_v_dot_v [protected]

Definition at line 96 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Tensor Controller::momentumSqrSum [protected]

Definition at line 186 of file Controller.h.

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

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

Definition at line 188 of file Controller.h.

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

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

Definition at line 189 of file Controller.h.

Referenced by Controller(), and multigratorTemperature().

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

Definition at line 190 of file Controller.h.

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

RequireReduction* Controller::multigratorReduction [protected]

Definition at line 192 of file Controller.h.

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

BigReal Controller::multigratorXi [protected]

Definition at line 184 of file Controller.h.

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

BigReal Controller::multigratorXiT [protected]

Definition at line 185 of file Controller.h.

Referenced by multigratorPressure().

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

Definition at line 191 of file Controller.h.

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

int Controller::nbondFreq [protected]

Definition at line 78 of file Controller.h.

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

BigReal Controller::net_dE [protected]

Definition at line 119 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::net_dEdl_bond_1 [protected]

Definition at line 133 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_bond_2 [protected]

Definition at line 134 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_elec_1 [protected]

Definition at line 135 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_elec_2 [protected]

Definition at line 136 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_lj_1 [protected]

Definition at line 137 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_lj_2 [protected]

Definition at line 138 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int64_t Controller::numDegFreedom [protected]

Definition at line 100 of file Controller.h.

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

Lattice Controller::origLattice [protected]

Definition at line 256 of file Controller.h.

Referenced by Controller().

int Controller::pairlistWarnings [protected]

Definition at line 90 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Tensor Controller::positionRescaleFactor [protected]

Definition at line 181 of file Controller.h.

Referenced by langevinPiston1().

PressureProfileReduction* Controller::ppbonded [protected]

Definition at line 217 of file Controller.h.

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

PressureProfileReduction* Controller::ppint [protected]

Definition at line 219 of file Controller.h.

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

PressureProfileReduction* Controller::ppnonbonded [protected]

Definition at line 218 of file Controller.h.

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

Tensor Controller::pressure [protected]

Definition at line 163 of file Controller.h.

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

Tensor Controller::pressure_amd [protected]

Definition at line 70 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

BigReal Controller::pressure_avg [protected]

Definition at line 81 of file Controller.h.

Referenced by Controller(), and printEnergies().

Tensor Controller::pressure_nbond [protected]

Definition at line 68 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_normal [protected]

Definition at line 67 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_slow [protected]

Definition at line 69 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_tavg [protected]

Definition at line 84 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal* Controller::pressureProfileAverage [protected]

Definition at line 222 of file Controller.h.

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

int Controller::pressureProfileCount [protected]

Definition at line 221 of file Controller.h.

Referenced by Controller(), and printEnergies().

int Controller::pressureProfileSlabs [protected]

Definition at line 220 of file Controller.h.

Referenced by Controller(), and printEnergies().

Random* Controller::random [protected]

Definition at line 209 of file Controller.h.

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

BigReal Controller::recent_alchWork [protected]

Definition at line 149 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_bond_1 [protected]

Definition at line 143 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_bond_2 [protected]

Definition at line 144 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_elec_1 [protected]

Definition at line 145 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_elec_2 [protected]

Definition at line 146 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_lj_1 [protected]

Definition at line 147 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_lj_2 [protected]

Definition at line 148 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::recent_TiNo [protected]

Definition at line 151 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

RequireReduction* Controller::reduction [protected]

Definition at line 212 of file Controller.h.

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

int Controller::rescaleVelocities_numTemps [protected]

Definition at line 171 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

BigReal Controller::rescaleVelocities_sumTemps [protected]

Definition at line 170 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

SimParameters* const Controller::simParams [protected]

Definition at line 210 of file Controller.h.

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

int Controller::slowFreq [protected]

Definition at line 79 of file Controller.h.

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

BigReal Controller::smooth2_avg2 [protected]

Definition at line 162 of file Controller.h.

NamdState* const Controller::state [protected]

Definition at line 211 of file Controller.h.

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

int Controller::stepInFullRun [protected]

Definition at line 101 of file Controller.h.

Referenced by printEnergies().

Tensor Controller::strainRate_old [protected]

Definition at line 180 of file Controller.h.

Referenced by langevinPiston1().

SubmitReduction* Controller::submit_reduction [protected]

Definition at line 214 of file Controller.h.

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

int Controller::tavg_count [protected]

Definition at line 86 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::temp_avg [protected]

Definition at line 80 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::temperature [protected]

Definition at line 160 of file Controller.h.

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

ofstream_namd Controller::tiFile [protected]

Definition at line 237 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::TiNo [protected]

Definition at line 142 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::totalEnergy [protected]

Definition at line 102 of file Controller.h.

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

Tensor Controller::virial_amd [protected]

Definition at line 71 of file Controller.h.

Referenced by calcPressure(), and rescaleaccelMD().

ofstream_namd Controller::xstFile [protected]

Definition at line 227 of file Controller.h.

Referenced by outputExtendedSystem().


The documentation for this class was generated from the following files:
Generated on Mon Sep 25 01:17:17 2017 for NAMD by  doxygen 1.4.7