#include <Controller.h>

Inheritance diagram for Controller:

Classes
struct	checkpoint

Public Member Functions
	Controller (NamdState *s)

virtual	~Controller (void)

void	run (void)

void	awaken (void)

void	resumeAfterTraceBarrier (int)

Public Attributes
BigReal	accelMDdV

Protected Member Functions
virtual void	algorithm (void)

void	integrate (int)

void	minimize ()

void	receivePressure (int step, int minimize=0)

void	calcPressure (int step, int minimize, const Tensor &virial_normal_in, const Tensor &virial_nbond_in, const Tensor &virial_slow_in, const Tensor &intVirial_normal, const Tensor &intVirial_nbond, const Tensor &intVirial_slow, const Vector &extForce_normal, const Vector &extForce_nbond, const Vector &extForce_slow)

void	compareChecksums (int, int=0)

void	printTiming (int)

void	printMinimizeEnergies (int)

void	printDynamicsEnergies (int)

void	printEnergies (int step, int minimize)

void	printFepMessage (int)

void	printTiMessage (int)

void	enqueueCollections (int)

void	correctMomentum (int step)

void	rescaleVelocities (int)

void	reassignVelocities (int)

void	tcoupleVelocities (int)

void	stochRescaleVelocities (int)

double	stochRescaleCoefficient ()

void	berendsenPressure (int)

void	langevinPiston1 (int)

void	langevinPiston2 (int)

void	multigratorPressure (int step, int callNumber)

void	multigratorTemperature (int step, int callNumber)

BigReal	multigatorCalcEnthalpy (BigReal potentialEnergy, int step, int minimize)

void	rebalanceLoad (int)

void	cycleBarrier (int, int)

void	traceBarrier (int, int)

void	terminate (void)

void	outputExtendedSystem (int step)

void	writeExtendedSystemLabels (ofstream_namd &file)

void	writeExtendedSystemData (int step, ofstream_namd &file)

void	outputFepEnergy (int step)

void	writeFepEnergyData (int step, ofstream_namd &file)

void	outputTiEnergy (int step)

BigReal	computeAlchWork (const int step)

void	writeTiEnergyData (int step, ofstream_namd &file)

void	recvCheckpointReq (const char *key, int task, checkpoint &cp)

void	recvCheckpointAck (checkpoint &cp)

void	calc_accelMDG_mean_std (BigReal testV, int step_n, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal M2, BigReal *sigmaV)

void	calc_accelMDG_E_k (int iE, int V_n, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal k0, BigReal k, BigReal E, int iEused, char *warn)

void	calc_accelMDG_force_factor (BigReal k, BigReal E, BigReal testV, Tensor vir_orig, BigReal dV, BigReal factor, Tensor *vir)

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

void	rescaleaccelMD (int step, int minimize=0)

void	adaptTempInit (int step)

void	adaptTempUpdate (int step, int minimize=0)

void	adaptTempWriteRestart (int step)

Protected Attributes
RequireReduction *	min_reduction

Tensor	pressure_normal

Tensor	pressure_nbond

Tensor	pressure_slow

Tensor	pressure_amd

Tensor	virial_amd

Tensor	groupPressure_normal

Tensor	groupPressure_nbond

Tensor	groupPressure_slow

Tensor	controlPressure_normal

Tensor	controlPressure_nbond

Tensor	controlPressure_slow

int	nbondFreq

int	slowFreq

BigReal	temp_avg

BigReal	pressure_avg

BigReal	groupPressure_avg

int	avg_count

Tensor	pressure_tavg

Tensor	groupPressure_tavg

int	tavg_count

int	computeChecksum

int	marginViolations

int	pairlistWarnings

BigReal	min_energy

BigReal	min_f_dot_f

BigReal	min_f_dot_v

BigReal	min_v_dot_v

int	min_huge_count

int64_t	numDegFreedom

int	stepInFullRun

BigReal	totalEnergy

BigReal	electEnergy

BigReal	electEnergySlow

BigReal	ljEnergy

BigReal	groLJEnergy

BigReal	groGaussEnergy

BigReal	goNativeEnergy

BigReal	goNonnativeEnergy

BigReal	goTotalEnergy

BigReal	bondedEnergyDiff_f

BigReal	electEnergy_f

BigReal	electEnergySlow_f

BigReal	ljEnergy_f

BigReal	ljEnergy_f_left

BigReal	exp_dE_ByRT

BigReal	dE

BigReal	net_dE

BigReal	dG

int	FepNo

BigReal	fepSum

BigReal	bondedEnergy_ti_1

BigReal	bondedEnergy_ti_2

BigReal	electEnergy_ti_1

BigReal	electEnergySlow_ti_1

BigReal	ljEnergy_ti_1

BigReal	electEnergy_ti_2

BigReal	electEnergySlow_ti_2

BigReal	ljEnergy_ti_2

BigReal	net_dEdl_bond_1

BigReal	net_dEdl_bond_2

BigReal	net_dEdl_elec_1

BigReal	net_dEdl_elec_2

BigReal	net_dEdl_lj_1

BigReal	net_dEdl_lj_2

BigReal	cumAlchWork

BigReal	electEnergyPME_ti_1

BigReal	electEnergyPME_ti_2

int	TiNo

BigReal	recent_dEdl_bond_1

BigReal	recent_dEdl_bond_2

BigReal	recent_dEdl_elec_1

BigReal	recent_dEdl_elec_2

BigReal	recent_dEdl_lj_1

BigReal	recent_dEdl_lj_2

BigReal	recent_alchWork

BigReal	alchWork

int	recent_TiNo

BigReal	drudeBondTemp

BigReal	drudeBondTempAvg

BigReal	kineticEnergy

BigReal	kineticEnergyHalfstep

BigReal	kineticEnergyCentered

BigReal	temperature

BigReal	heat

BigReal	totalEnergy0

BigReal	smooth2_avg2

Tensor	pressure

Tensor	groupPressure

int	controlNumDegFreedom

Tensor	controlPressure

BigReal	rescaleVelocities_sumTemps

int	rescaleVelocities_numTemps

int	stochRescale_count

BigReal	stochRescaleTimefactor

Tensor	langevinPiston_origStrainRate

Tensor	strainRate_old

Tensor	positionRescaleFactor

BigReal	multigratorXi

BigReal	multigratorXiT

Tensor	momentumSqrSum

std::vector< BigReal >	multigratorNu

std::vector< BigReal >	multigratorNuT

std::vector< BigReal >	multigratorOmega

std::vector< BigReal >	multigratorZeta

RequireReduction *	multigratorReduction

int	ldbSteps

int	fflush_count

Random *	random

SimParameters *const	simParams

NamdState *const	state

RequireReduction *	reduction

RequireReduction *	amd_reduction

SubmitReduction *	submit_reduction

PressureProfileReduction *	ppbonded

PressureProfileReduction *	ppnonbonded

PressureProfileReduction *	ppint

int	pressureProfileSlabs

int	pressureProfileCount

BigReal *	pressureProfileAverage

CollectionMaster *const	collection

ControllerBroadcasts *	broadcast

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

BigReal *	adaptTempPotEnergyAveNum

BigReal *	adaptTempPotEnergyAveDen

BigReal *	adaptTempPotEnergyVarNum

BigReal *	adaptTempPotEnergyAve

BigReal *	adaptTempPotEnergyVar

int *	adaptTempPotEnergySamples

BigReal *	adaptTempBetaN

BigReal	adaptTempT

BigReal	adaptTempDTave

BigReal	adaptTempDTavenum

BigReal	adaptTempBetaMin

BigReal	adaptTempBetaMax

int	adaptTempBin

int	adaptTempBins

BigReal	adaptTempDBeta

BigReal	adaptTempCg

BigReal	adaptTempDt

Bool	adaptTempAutoDt

BigReal	adaptTempDtMin

BigReal	adaptTempDtMax

ofstream_namd	adaptTempRestartFile

Protected Attributes inherited from ControllerState
Tensor	langevinPiston_strainRate

Tensor	berendsenPressure_avg

int	berendsenPressure_count

BigReal	smooth2_avg

Friends
class	ScriptTcl

class	Node

class	CheckpointMsg

Detailed Description

Definition at line 39 of file Controller.h.

Constructor & Destructor Documentation

Controller::Controller ( NamdState * s )

Definition at line 128 of file Controller.C.

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

                                    :
         computeChecksum(0), marginViolations(0), pairlistWarnings(0),
         simParams(Node::Object()->simParameters),
         state(s),
         collection(CollectionMaster::Object()),
         startCTime(0),
         firstCTime(CmiTimer()),
         startWTime(0),
         firstWTime(CmiWallTimer()),
         startBenchTime(0),
         stepInFullRun(0),
         ldbSteps(0),
         fflush_count(3)
 {
     broadcast = new ControllerBroadcasts;
     min_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MINIMIZER,1);
     // for accelMD
     if (simParams->accelMDOn) {
        // REDUCTIONS_BASIC wil contain all potential energies and arrive first
        amd_reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
        // contents of amd_reduction be submitted to REDUCTIONS_AMD
        submit_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_AMD);
        // REDUCTIONS_AMD will contain Sequencer contributions and arrive second
        reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_AMD);
     } else {
        amd_reduction = NULL;
        submit_reduction = NULL;
        reduction = ReductionMgr::Object()->willRequire(REDUCTIONS_BASIC);
     }
     // pressure profile reductions
     pressureProfileSlabs = 0;
     pressureProfileCount = 0;
     ppbonded = ppnonbonded = ppint = NULL;
     if (simParams->pressureProfileOn || simParams->pressureProfileEwaldOn) {
       int ntypes = simParams->pressureProfileAtomTypes;
       int nslabs = pressureProfileSlabs = simParams->pressureProfileSlabs;
       // number of partitions for pairwise interactions
       int npairs = (ntypes * (ntypes+1))/2;
       pressureProfileAverage = new BigReal[3*nslabs];
       memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
       int freq = simParams->pressureProfileFreq;
       if (simParams->pressureProfileOn) {
         ppbonded = new PressureProfileReduction(REDUCTIONS_PPROF_BONDED, 
             nslabs, npairs, "BONDED", freq);
         ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
             nslabs, npairs, "NONBONDED", freq);
         // internal partitions by atom type, but not in a pairwise fashion
         ppint = new PressureProfileReduction(REDUCTIONS_PPROF_INTERNAL, 
             nslabs, ntypes, "INTERNAL", freq);
       } else {
         // just doing Ewald, so only need nonbonded
         ppnonbonded = new PressureProfileReduction(REDUCTIONS_PPROF_NONBONDED, 
             nslabs, npairs, "NONBONDED", freq);
       }
     }
     random = new Random(simParams->randomSeed);
     random->split(0,PatchMap::Object()->numPatches()+1);
 
     heat = totalEnergy0 = 0.0;
     rescaleVelocities_sumTemps = 0;  
     rescaleVelocities_numTemps = 0;
     stochRescale_count = 0;
     if (simParams->stochRescaleOn) {
       stochRescaleTimefactor = exp(-simParams->stochRescaleFreq *
           simParams->dt * 0.001 / simParams->stochRescalePeriod);
     }
     berendsenPressure_avg = 0; berendsenPressure_count = 0;
     // strainRate tensor is symmetric to avoid rotation
     langevinPiston_strainRate =
         Tensor::symmetric(simParams->strainRate,simParams->strainRate2);
     if ( ! simParams->useFlexibleCell ) {
       BigReal avg = trace(langevinPiston_strainRate) / 3.;
       langevinPiston_strainRate = Tensor::identity(avg);
     } else if ( simParams->useConstantRatio ) {
 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
                  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
       AVGXY(langevinPiston_strainRate);
 #undef AVGXY
     }
     langevinPiston_origStrainRate = langevinPiston_strainRate;
     if (simParams->multigratorOn) {
       multigratorXi = 0.0;
       int n = simParams->multigratorNoseHooverChainLength;
       BigReal tau = simParams->multigratorTemperatureRelaxationTime;
       Node *node = Node::Object();
       Molecule *molecule = node->molecule;
       BigReal Nf = molecule->num_deg_freedom();
       BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
       multigratorNu.resize(n);
       multigratorNuT.resize(n);
       multigratorZeta.resize(n);
       multigratorOmega.resize(n);
       for (int i=0;i < n;i++) {
         multigratorNu[i] = 0.0;
         multigratorZeta[i] = 0.0;
         if (i == 0) {
           multigratorOmega[i] = Nf*kT0*tau*tau;
         } else {
           multigratorOmega[i] = kT0*tau*tau;
         }
       }
       multigratorReduction = ReductionMgr::Object()->willRequire(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
     } else {
       multigratorReduction = NULL;
     }
     origLattice = state->lattice;
     smooth2_avg = XXXBIGREAL;
     temp_avg = 0;
     pressure_avg = 0;
     groupPressure_avg = 0;
     avg_count = 0;
     pressure_tavg = 0;
     groupPressure_tavg = 0;
     tavg_count = 0;
     checkpoint_stored = 0;
     drudeBondTemp = 0;
     drudeBondTempAvg = 0;
     cumAlchWork = 0;
 }

Controller::~Controller ( void )

virtual

Definition at line 248 of file Controller.C.

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

 {
     delete broadcast;
     delete reduction;
     delete min_reduction;
     delete amd_reduction;
     delete submit_reduction;
     delete ppbonded;
     delete ppnonbonded;
     delete ppint;
     delete [] pressureProfileAverage;
     delete random;
     if (multigratorReduction) delete multigratorReduction;
 }

Member Function Documentation

void Controller::adaptTempInit ( int step )

protected

Definition at line 2351 of file Controller.C.

Referenced by adaptTempUpdate().

                                        {
     if (!simParams->adaptTempOn) return;
     iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
     adaptTempDtMin = 0;
     adaptTempDtMax = 0;
     adaptTempAutoDt = false;
     if (simParams->adaptTempBins == 0) {
       iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
       std::ifstream adaptTempRead(simParams->adaptTempInFile);
       if (adaptTempRead) {
         int readInt;
         BigReal readReal;
         bool readBool;
         // step
         adaptTempRead >> readInt;
         // Start with min and max temperatures
         adaptTempRead >> adaptTempT;     // KELVIN
         adaptTempRead >> adaptTempBetaMin;  // KELVIN
         adaptTempRead >> adaptTempBetaMax;  // KELVIN
         adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
         adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
         // In case file is manually edited
         if (adaptTempBetaMin > adaptTempBetaMax){
             readReal = adaptTempBetaMax;
             adaptTempBetaMax = adaptTempBetaMin;
             adaptTempBetaMin = adaptTempBetaMax;
         }
         adaptTempRead >> adaptTempBins;     
         adaptTempRead >> adaptTempCg;
         adaptTempRead >> adaptTempDt;
         adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
         adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
         adaptTempPotEnergySamples = new int[adaptTempBins];
         adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
         adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
         adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
         adaptTempBetaN           = new BigReal[adaptTempBins];
         adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
         for(int j = 0; j < adaptTempBins; ++j) {
           adaptTempRead >> adaptTempPotEnergyAve[j];
           adaptTempRead >> adaptTempPotEnergyVar[j];
           adaptTempRead >> adaptTempPotEnergySamples[j];
           adaptTempRead >> adaptTempPotEnergyAveNum[j];
           adaptTempRead >> adaptTempPotEnergyVarNum[j];
           adaptTempRead >> adaptTempPotEnergyAveDen[j];
           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
         } 
         adaptTempRead.close();
       }
       else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
     } 
     else {
       adaptTempBins = simParams->adaptTempBins;
       adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
       adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
       adaptTempPotEnergySamples = new int[adaptTempBins];
       adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
       adaptTempPotEnergyVar    = new BigReal[adaptTempBins];
       adaptTempPotEnergyAve    = new BigReal[adaptTempBins];
       adaptTempBetaN           = new BigReal[adaptTempBins];
       adaptTempBetaMax = 1./simParams->adaptTempTmin;
       adaptTempBetaMin = 1./simParams->adaptTempTmax;
       adaptTempCg = simParams->adaptTempCgamma;   
       adaptTempDt  = simParams->adaptTempDt;
       adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
       adaptTempT = simParams->initialTemp; 
       if (simParams->langevinOn)
         adaptTempT = simParams->langevinTemp;
       else if (simParams->rescaleFreq > 0)
         adaptTempT = simParams->rescaleTemp;
       for(int j = 0; j < adaptTempBins; ++j){
           adaptTempPotEnergyAveNum[j] = 0.;
           adaptTempPotEnergyAveDen[j] = 0.;
           adaptTempPotEnergySamples[j] = 0;
           adaptTempPotEnergyVarNum[j] = 0.;
           adaptTempPotEnergyVar[j] = 0.;
           adaptTempPotEnergyAve[j] = 0.;
           adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
       }
     }
     if (simParams->adaptTempAutoDt > 0.0) {
        adaptTempAutoDt = true;
        adaptTempDtMin =  simParams->adaptTempAutoDt - 0.01;
        adaptTempDtMax =  simParams->adaptTempAutoDt + 0.01;
     }
     adaptTempDTave = 0;
     adaptTempDTavenum = 0;
     iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
     iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
     iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
     iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
     iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
     if (simParams->adaptTempRestartFreq > 0) {
       NAMD_backup_file(simParams->adaptTempRestartFile);
       adaptTempRestartFile.open(simParams->adaptTempRestartFile);
        if(!adaptTempRestartFile)
         NAMD_die("Error opening restart file");
     }
 }

void Controller::adaptTempUpdate	(	int	step,
		int	minimize = `0`
	)

protected

Definition at line 2477 of file Controller.C.

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

Referenced by integrate().

 {
     //Beta = 1./T
     if ( !simParams->adaptTempOn ) return;
     int j = 0;
     if (step == simParams->firstTimestep) {
         adaptTempInit(step);
         return;
     }
     if ( minimize || (step < simParams->adaptTempFirstStep ) || 
         ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
     const int adaptTempOutFreq  = simParams->adaptTempOutFreq;
     const bool adaptTempDebug  = simParams->adaptTempDebug;
     //Calculate Current inverse temperature and bin 
     BigReal adaptTempBeta = 1./adaptTempT;
     adaptTempBin   = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
 
     if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
         iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin  
                                << " adaptTempBeta: " << adaptTempBeta 
                               << " adaptTempDBeta: " << adaptTempDBeta 
                                << " betaMin:" << adaptTempBetaMin 
                                << " betaMax: " << adaptTempBetaMax << "\n";
     adaptTempPotEnergySamples[adaptTempBin] += 1;
     BigReal gammaAve = 1.-adaptTempCg/adaptTempPotEnergySamples[adaptTempBin];
 
     BigReal potentialEnergy;
     BigReal potEnergyAverage;
     BigReal potEnergyVariance;
     potentialEnergy = totalEnergy-kineticEnergy;
 
     //calculate new bin average and variance using adaptive averaging
     adaptTempPotEnergyAveNum[adaptTempBin] = adaptTempPotEnergyAveNum[adaptTempBin]*gammaAve + potentialEnergy;
     adaptTempPotEnergyAveDen[adaptTempBin] = adaptTempPotEnergyAveDen[adaptTempBin]*gammaAve + 1;
     adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
     
     potEnergyAverage = adaptTempPotEnergyAveNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
     potEnergyVariance = adaptTempPotEnergyVarNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
     potEnergyVariance -= potEnergyAverage*potEnergyAverage;
 
     adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
     adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
     
     // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
     // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
     // is constant for each bin. This is to estimate <E(beta)> where beta \in
     // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
     if ( ! ( step % simParams->adaptTempFreq ) ) {
       // If adaptTempBin not at the edge, calculate improved average:
       if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
           // Get Averaging Limits:
           BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
           BigReal betaPlus;
           BigReal betaMinus;
           int     nMinus =0;
           int     nPlus = 0;
           if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
           if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
           betaMinus = adaptTempBeta - deltaBeta;
           betaPlus =  adaptTempBeta + deltaBeta;
           BigReal betaMinus2 = betaMinus*betaMinus;
           BigReal betaPlus2  = betaPlus*betaPlus;
           nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
           nPlus  = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
           // Variables for <E(beta)> estimate:
           BigReal potEnergyAve0 = 0.0;
           BigReal potEnergyAve1 = 0.0;
           // Integral terms
           BigReal A0 = 0;
           BigReal A1 = 0;
           BigReal A2 = 0;
           //A0 phi_s integral for beta_minus < beta < beta_{i+1}
           BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1]; 
           BigReal DeltaE2Ave;
           BigReal DeltaE2AveSum = 0;
           BigReal betaj;
           for (j = nMinus+1; j <= adaptTempBin; ++j) {
               potEnergyAve0 += adaptTempPotEnergyAve[j]; 
               DeltaE2Ave = adaptTempPotEnergyVar[j];
               DeltaE2AveSum += DeltaE2Ave;
               A0 += j*DeltaE2Ave;
           }
           A0 *= adaptTempDBeta;
           A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
           A0 *= adaptTempDBeta;
           betaj = adaptTempBetaN[nMinus+1]-betaMinus; 
           betaj *= betaj;
           A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
           A0 /= (betaNp1-betaMinus);
 
           //A1 phi_s integral for beta_{i+1} < beta < beta_plus
           DeltaE2AveSum = 0;
           for (j = adaptTempBin+1; j < nPlus; j++) {
               potEnergyAve1 += adaptTempPotEnergyAve[j];
               DeltaE2Ave = adaptTempPotEnergyVar[j];
               DeltaE2AveSum += DeltaE2Ave;
               A1 += j*DeltaE2Ave;
           }
           A1 *= adaptTempDBeta;   
           A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
           A1 *= adaptTempDBeta;
           if ( nPlus < adaptTempBins ) {
             betaj = betaPlus-adaptTempBetaN[nPlus];
             betaj *= betaj;
             A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
           }
           A1 /= (betaPlus-betaNp1);
 
           //A2 phi_t integral for beta_i
           A2 = 0.5*adaptTempDBeta*potEnergyVariance;
 
           // Now calculate a+ and a-
           DeltaE2Ave = A0-A1;
           BigReal aplus = 0;
           BigReal aminus = 0;
           if (DeltaE2Ave != 0) {
             aplus = (A0-A2)/(A0-A1);
             if (aplus < 0) {
                     aplus = 0;
             }
             if (aplus > 1)  {
                     aplus = 1;
             }
             aminus = 1-aplus;
             potEnergyAve0 *= adaptTempDBeta;
             potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
             potEnergyAve0 /= (betaNp1-betaMinus);
             potEnergyAve1 *= adaptTempDBeta;
             if ( nPlus < adaptTempBins ) {
                 potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
             }
             potEnergyAve1 /= (betaPlus-betaNp1);
             potEnergyAverage = aminus*potEnergyAve0;
             potEnergyAverage += aplus*potEnergyAve1;
           }
           if (simParams->adaptTempDebug) {
        iout << "ADAPTEMP DEBUG:"  << "\n"
             << "     adaptTempBin:    " << adaptTempBin << "\n"
             << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
             << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
             << "     adaptTemp:   " << adaptTempT<< "\n"
             << "     betaMin:   " << adaptTempBetaMin << "\n"
             << "     betaMax:   " << adaptTempBetaMax << "\n"
             << "     gammaAve:  " << gammaAve << "\n"
             << "     deltaBeta: " << deltaBeta << "\n"
             << "     betaMinus: " << betaMinus << "\n"
             << "     betaPlus:  " << betaPlus << "\n"
             << "     nMinus:    " << nMinus << "\n"
             << "     nPlus:     " << nPlus << "\n"
             << "     A0:        " << A0 << "\n"
             << "     A1:        " << A1 << "\n"
             << "     A2:        " << A2 << "\n"
             << "     a+:        " << aplus << "\n"
             << "     a-:        " << aminus << "\n"
             << endi;
           }
       }
       else {
           if (simParams->adaptTempDebug) {
        iout << "ADAPTEMP DEBUG:"  << "\n"
             << "     adaptTempBin:    " << adaptTempBin << "\n"
             << "     Samples:   " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
             << "     adaptTempBeta:   " << adaptTempBeta << "\n" 
             << "     adaptTemp:   " << adaptTempT<< "\n"
             << "     betaMin:   " << adaptTempBetaMin << "\n"
             << "     betaMax:   " << adaptTempBetaMax << "\n"
             << "     gammaAve:  " << gammaAve << "\n"
             << "     deltaBeta:  N/A\n"
             << "     betaMinus:  N/A\n"
             << "     betaPlus:   N/A\n"
             << "     nMinus:     N/A\n"
             << "     nPlus:      N/A\n"
             << "     A0:         N/A\n"
             << "     A1:         N/A\n"
             << "     A2:         N/A\n"
             << "     a+:         N/A\n"
             << "     a-:         N/A\n"
             << endi;
           }
       }
       
       //dT is new temperature
       BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
       dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
       dT += adaptTempT;
       
      // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
      // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
       if ( dT > 1./adaptTempBetaMin || dT  < 1./adaptTempBetaMax ) {
         dT = ((potentialEnergy-adaptTempPotEnergyAve[adaptTempBin])/BOLTZMANN+adaptTempT)*adaptTempDt;
         dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
         dT += adaptTempT;
         // Check again, if not then keep original adaptTempTor assign random.
         if ( dT > 1./adaptTempBetaMin ) {
           if (!simParams->adaptTempRandom) {             
              //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
              //     << " K higher than adaptTempTmax."
              //     << " Keeping temperature at " 
              //     << adaptTempT<< "\n"<< endi;             
              dT = adaptTempT;
           }
           else {             
              //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT 
              //     << " K higher than adaptTempTmax."
              //     << " Assigning random temperature in range\n" << endi;
              dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);             
              dT = 1./dT;
           }
         } 
         else if ( dT  < 1./adaptTempBetaMax ) {
           if (!simParams->adaptTempRandom) {            
             //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
             //     << " K lower than adaptTempTmin."
             //     << " Keeping temperature at " << adaptTempT<< "\n" << endi; 
             dT = adaptTempT;
           }
           else {
             //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT 
             //     << " K lower than adaptTempTmin."
             //     << " Assigning random temperature in range\n" << endi;
             dT = adaptTempBetaMin +  random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
             dT = 1./dT;
           }
         }
         else if (adaptTempAutoDt) {
           //update temperature step size counter
           //FOR "TRUE" ADAPTIVE TEMPERING 
           BigReal adaptTempTdiff = fabs(dT-adaptTempT);
           if (adaptTempTdiff > 0) {
             adaptTempDTave += adaptTempTdiff; 
             adaptTempDTavenum++;
 //            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
           }
           if(adaptTempDTavenum == 100){
                 BigReal Frac;
                 adaptTempDTave /= adaptTempDTavenum;
                 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
                 Frac = adaptTempDTave/Frac;
                 //if average temperature jump is > 3% of temperature range,
                 //modify jump size to match 3%
                 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
                 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
                     Frac = adaptTempDtMax/Frac;
                     iout << "ADAPTEMP: Updating adaptTempDt to ";
                     adaptTempDt *= Frac;
                     iout << adaptTempDt << "\n" << endi;
                 }
                 adaptTempDTave = 0;
                 adaptTempDTavenum = 0;
           }
         }
       }
       else if (adaptTempAutoDt) {
           //update temperature step size counter
           // FOR "TRUE" ADAPTIVE TEMPERING
           BigReal adaptTempTdiff = fabs(dT-adaptTempT);
           if (adaptTempTdiff > 0) {
             adaptTempDTave += adaptTempTdiff; 
             adaptTempDTavenum++;
 //            iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
           }
           if(adaptTempDTavenum == 100){
                 BigReal Frac;
                 adaptTempDTave /= adaptTempDTavenum;
                 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
                 Frac = adaptTempDTave/Frac;
                 //if average temperature jump is > 3% of temperature range,
                 //modify jump size to match 3%
                 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n"; 
                 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
                     Frac = adaptTempDtMax/Frac;
                     iout << "ADAPTEMP: Updating adaptTempDt to ";
                     adaptTempDt *= Frac;
                     iout << adaptTempDt << "\n" << endi;
                 }
                 adaptTempDTave = 0;
                 adaptTempDTavenum = 0;
 
           }
           
       }
       
       adaptTempT = dT; 
       broadcast->adaptTemperature.publish(step,adaptTempT);
     }
     adaptTempWriteRestart(step);
     if ( ! (step % adaptTempOutFreq) ) {
         iout << "ADAPTEMP: STEP " << step
              << " TEMP "   << adaptTempT
              << " ENERGY " << std::setprecision(10) << potentialEnergy   
              << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
              << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
         iout << "\n" << endi;
    }
    
 }

void Controller::adaptTempWriteRestart ( int step )

protected

Definition at line 2451 of file Controller.C.

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

Referenced by adaptTempUpdate().

                                                {
     if (simParams->adaptTempOn && !(step%simParams->adaptTempRestartFreq)) {
         adaptTempRestartFile.seekp(std::ios::beg);        
         iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
         adaptTempRestartFile << step << " ";
         // Start with min and max temperatures
         adaptTempRestartFile << adaptTempT<< " ";     // KELVIN
         adaptTempRestartFile << 1./adaptTempBetaMin << " ";  // KELVIN
         adaptTempRestartFile << 1./adaptTempBetaMax << " ";  // KELVIN
         adaptTempRestartFile << adaptTempBins << " ";     
         adaptTempRestartFile << adaptTempCg << " ";
         adaptTempRestartFile << adaptTempDt ;
         adaptTempRestartFile << "\n" ;
         for(int j = 0; j < adaptTempBins; ++j) {
           adaptTempRestartFile << adaptTempPotEnergyAve[j] << " ";
           adaptTempRestartFile << adaptTempPotEnergyVar[j] << " ";
           adaptTempRestartFile << adaptTempPotEnergySamples[j] << " ";
           adaptTempRestartFile << adaptTempPotEnergyAveNum[j] << " ";
           adaptTempRestartFile << adaptTempPotEnergyVarNum[j] << " ";
           adaptTempRestartFile << adaptTempPotEnergyAveDen[j] << " ";
           adaptTempRestartFile << "\n";          
         }
         adaptTempRestartFile.flush(); 
     }
 }    

void Controller::algorithm ( void )

protectedvirtual

Definition at line 281 of file Controller.C.

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

 {
   int scriptTask;
   int scriptSeq = 0;
   BackEnd::awaken();
   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
     switch ( scriptTask ) {
       case SCRIPT_OUTPUT:
         enqueueCollections(FILE_OUTPUT);
         outputExtendedSystem(FILE_OUTPUT);
         break;
       case SCRIPT_FORCEOUTPUT:
         enqueueCollections(FORCE_OUTPUT);
         break;
       case SCRIPT_MEASURE:
         enqueueCollections(EVAL_MEASURE);
         break;
       case SCRIPT_REINITVELS:
         iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
           << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
         break;
       case SCRIPT_RESCALEVELS:
         iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
           << " BY " << simParams->scriptArg1 << "\n" << endi;
         break;
       case SCRIPT_RESCALESOLUTECHARGES:
         // Parameter setting already reported in ScriptTcl
         // Nothing to do!
         break;
       case SCRIPT_CHECKPOINT:
         iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
           << "\n" << endi;
         checkpoint_stored = 1;
         checkpoint_lattice = state->lattice;
         checkpoint_state = *(ControllerState*)this;
         break;
       case SCRIPT_REVERT:
         iout << "REVERTING AT STEP " << simParams->firstTimestep
           << "\n" << endi;
         if ( ! checkpoint_stored )
           NAMD_die("Unable to revert, checkpoint was never called!");
         state->lattice = checkpoint_lattice;
         *(ControllerState*)this = checkpoint_state;
         break;
       case SCRIPT_CHECKPOINT_STORE:
         iout << "STORING CHECKPOINT AT STEP " << simParams->firstTimestep
           << " TO KEY " << simParams->scriptStringArg1;
         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
         iout << "\n" << endi;
         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
             checkpoints[simParams->scriptStringArg1] = new checkpoint;
           }
           checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
           cp.lattice = state->lattice;
           cp.state = *(ControllerState*)this;
         } else {
           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                             scriptTask, state->lattice, *(ControllerState*)this);
         }
         break;
       case SCRIPT_CHECKPOINT_LOAD:
         iout << "LOADING CHECKPOINT AT STEP " << simParams->firstTimestep
           << " FROM KEY " << simParams->scriptStringArg1;
         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
         iout << "\n" << endi;
         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
             NAMD_die("Unable to load checkpoint, requested key was never stored.");
           }
           checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
           state->lattice = cp.lattice;
           *(ControllerState*)this = cp.state;
         } else {
           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                             scriptTask, state->lattice, *(ControllerState*)this);
         }
         break;
       case SCRIPT_CHECKPOINT_SWAP:
         iout << "SWAPPING CHECKPOINT AT STEP " << simParams->firstTimestep
           << " FROM KEY " << simParams->scriptStringArg1;
         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
         iout << "\n" << endi;
         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
             NAMD_die("Unable to swap checkpoint, requested key was never stored.");
           }
           checkpoint &cp = *checkpoints[simParams->scriptStringArg1];
           std::swap(state->lattice,cp.lattice);
           std::swap(*(ControllerState*)this,cp.state);
         } else {
           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                             scriptTask, state->lattice, *(ControllerState*)this);
         }
         break;
       case SCRIPT_CHECKPOINT_FREE:
         iout << "FREEING CHECKPOINT AT STEP " << simParams->firstTimestep
           << " FROM KEY " << simParams->scriptStringArg1;
         if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
         iout << "\n" << endi;
         if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
           if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
             NAMD_die("Unable to free checkpoint, requested key was never stored.");
           }
           delete checkpoints[simParams->scriptStringArg1];
           checkpoints.erase(simParams->scriptStringArg1);
         } else {
           Node::Object()->sendCheckpointReq(simParams->scriptIntArg1,simParams->scriptStringArg1,
                                             scriptTask, state->lattice, *(ControllerState*)this);
         }
         break;
       case SCRIPT_ATOMSENDRECV:
       case SCRIPT_ATOMSEND:
       case SCRIPT_ATOMRECV:
         break;
       case SCRIPT_MINIMIZE:
         minimize();
         break;
       case SCRIPT_RUN:
       case SCRIPT_CONTINUE:
         integrate(scriptTask);
         break;
       default:
         NAMD_bug("Unknown task in Controller::algorithm");
     }
     BackEnd::awaken();
   }
   enqueueCollections(END_OF_RUN);
   outputExtendedSystem(END_OF_RUN);
   terminate();
 }

void Controller::awaken ( void )

inline

Definition at line 45 of file Controller.h.

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

45 { CthAwaken(thread); };

void Controller::berendsenPressure ( int step )

protected

Definition at line 950 of file Controller.C.

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

Referenced by integrate().

 {
   if ( simParams->berendsenPressureOn ) {
    berendsenPressure_count += 1;
    berendsenPressure_avg += controlPressure;
    const int freq = simParams->berendsenPressureFreq;
    if ( ! (berendsenPressure_count % freq) ) {
     Tensor factor = berendsenPressure_avg / berendsenPressure_count;
     berendsenPressure_avg = 0;
     berendsenPressure_count = 0;
     // We only use on-diagonal terms (for now)
     factor = Tensor::diagonal(diagonal(factor));
     factor -= Tensor::identity(simParams->berendsenPressureTarget);
     factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
        simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
     factor += Tensor::identity(1.0);
 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
          if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
          if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
     int limited = 0;
     LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
     LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
     LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
 #undef LIMIT_SCALING
     if ( limited ) {
       iout << iERROR << "Step " << step <<
         " cell rescaling factor limited.\n" << endi;
     }
     broadcast->positionRescaleFactor.publish(step,factor);
     state->lattice.rescale(factor);
    }
   } else {
     berendsenPressure_avg = 0;
     berendsenPressure_count = 0;
   }
 }

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

inlineprotected

Definition at line 1736 of file Controller.C.

Referenced by rescaleaccelMD().

                                                               {
    switch(iE){
       case 2:
          *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
          // if k0 <=0 OR k0 > 1, switch to iE=1 mode
          if(*k0 > 1.)
             *warn |= 1;
          else if(*k0 <= 0.)
             *warn |= 2;
          //else stay in iE=2 mode
          else{
             *E = Vmin + (Vmax-Vmin)/(*k0);
             *iEused = 2;
             break;
          }
       case 1:
          *E = Vmax;
          *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
          if(*k0 > 1.)
             *k0 = 1.;
          *iEused = 1;
          break;
    }
 
    *k = *k0 / (Vmax - Vmin);
 }

void Controller::calc_accelMDG_force_factor	(	BigReal	k,
		BigReal	E,
		BigReal	testV,
		Tensor	vir_orig,
		BigReal *	dV,
		BigReal *	factor,
		Tensor *	vir
	)

inlineprotected

Definition at line 1765 of file Controller.C.

Referenced by rescaleaccelMD().

                                            {
    BigReal VE = testV - E;
    //if V < E, apply boost
    if(VE < 0){
       *factor = k * VE;
       *vir += vir_orig * (*factor);
       *dV += (*factor) * VE/2.;
       *factor += 1.;
    }
    else{
       *factor = 1.;
    }
 }

void Controller::calc_accelMDG_mean_std	(	BigReal	testV,
		int	step_n,
		BigReal *	Vmax,
		BigReal *	Vmin,
		BigReal *	Vavg,
		BigReal *	M2,
		BigReal *	sigmaV
	)

inlineprotected

Definition at line 1716 of file Controller.C.

Referenced by rescaleaccelMD().

                                                                            {
    BigReal delta; 
 
    if(testV > *Vmax){
       *Vmax = testV;
    }
    else if(testV < *Vmin){
       *Vmin = testV;
    }
 
    //mean and std calculated by Online Algorithm
    delta = testV - *Vavg;
    *Vavg += delta / (BigReal)n;
    *M2 += delta * (testV - (*Vavg));
 
    *sigmaV = sqrt(*M2/n);
 }

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

protected

Definition at line 1598 of file Controller.C.

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

Referenced by multigratorPressure(), and receivePressure().

                                                                                             {
 
   Tensor virial_normal = virial_normal_in;
   Tensor virial_nbond = virial_nbond_in;
   Tensor virial_slow = virial_slow_in;
 
   // Tensor tmp = virial_normal;
   // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
   //   tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
 
   Node *node = Node::Object();
   Molecule *molecule = node->molecule;
   SimParameters *simParameters = node->simParameters;
   Lattice &lattice = state->lattice;
 
   BigReal volume;
 
   Vector extPosition = lattice.origin();
   virial_normal -= outer(extForce_normal,extPosition);
   virial_nbond -= outer(extForce_nbond,extPosition);
   virial_slow -= outer(extForce_slow,extPosition);
 
   if ( (volume=lattice.volume()) != 0. )
   {
 
     if (simParameters->LJcorrection && volume) {
 #ifdef MEM_OPT_VERSION
       NAMD_bug("LJcorrection not supported in memory optimized build.");
 #else
       // Apply tail correction to pressure
       BigReal alchLambda = simParameters->getCurrentLambda(step);
       virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
 #endif
     }
 
     // kinetic energy component included in virials
     pressure_normal = virial_normal / volume;
     groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
 
     if (simParameters->accelMDOn) {
       pressure_amd = virial_amd / volume;
       pressure_normal += pressure_amd;
       groupPressure_normal +=  pressure_amd;
     }
 
     if ( minimize || ! ( step % nbondFreq ) )
     {
       pressure_nbond = virial_nbond / volume;
       groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
     }
 
     if ( minimize || ! ( step % slowFreq ) )
     {
       pressure_slow = virial_slow / volume;
       groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
     }
 
     pressure = pressure_normal + pressure_nbond + pressure_slow; 
     groupPressure = groupPressure_normal + groupPressure_nbond +
           groupPressure_slow;
   }
   else
   {
     pressure = Tensor();
     groupPressure = Tensor();
   }
 
   if ( simParameters->useGroupPressure )
   {
     controlPressure_normal = groupPressure_normal;
     controlPressure_nbond = groupPressure_nbond;
     controlPressure_slow = groupPressure_slow;
     controlPressure = groupPressure;
   }
   else
   {
     controlPressure_normal = pressure_normal;
     controlPressure_nbond = pressure_nbond;
     controlPressure_slow = pressure_slow;
     controlPressure = pressure;
   }
 
   if ( simParameters->useFlexibleCell ) {
     // use symmetric pressure to control rotation
     // controlPressure_normal = symmetric(controlPressure_normal);
     // controlPressure_nbond = symmetric(controlPressure_nbond);
     // controlPressure_slow = symmetric(controlPressure_slow);
     // controlPressure = symmetric(controlPressure);
     // only use on-diagonal components for now
     controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
     controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
     controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
     controlPressure = Tensor::diagonal(diagonal(controlPressure));
     if ( simParameters->useConstantRatio ) {
 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
    T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
       AVGXY(controlPressure_normal);
       AVGXY(controlPressure_nbond);
       AVGXY(controlPressure_slow);
       AVGXY(controlPressure);
 #undef AVGXY
     }
   } else {
     controlPressure_normal =
   Tensor::identity(trace(controlPressure_normal)/3.);
     controlPressure_nbond =
   Tensor::identity(trace(controlPressure_nbond)/3.);
     controlPressure_slow =
   Tensor::identity(trace(controlPressure_slow)/3.);
     controlPressure =
   Tensor::identity(trace(controlPressure)/3.);
   }
 }

void Controller::compareChecksums	(	int	step,
		int	forgiving = `0`
	)

protected

Definition at line 2775 of file Controller.C.

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

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

                                                          {
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
 
     // Some basic correctness checking
     BigReal checksum, checksum_b;
     char errmsg[256];
 
     checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
     if ( ((int)checksum) != molecule->numAtoms ) {
       sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
               (int)checksum, molecule->numAtoms);
       NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
     if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
       if ( ! computeChecksum ) {
         computesPartitioned = 0;
       } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
         computesPartitioned = 1;
       } else {
         NAMD_bug("Bad global compute count!\n");
       }
       computeChecksum = ((int)checksum);
     }
 
     checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
       sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcBonds);
       if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
         iout << iWARN << errmsg << endi;
       else NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
       sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcAngles);
       if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
         iout << iWARN << errmsg << endi;
       else NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
       sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcDihedrals);
       if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
         iout << iWARN << errmsg << endi;
       else NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
       sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcImpropers);
       if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
         iout << iWARN << errmsg << endi;
       else NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
       sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcTholes);
       if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
         iout << iWARN << errmsg << endi;
       else NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
       sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcAnisos);
       if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
         iout << iWARN << errmsg << endi;
       else NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
       sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcCrossterms);
       if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
         iout << iWARN << errmsg << endi;
       else NAMD_bug(errmsg);
     }
 
     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
     if ( ((int64)checksum) > molecule->numCalcExclusions &&
          ( ! simParams->mollyOn || step % slowFreq ) ) {
       if ( forgiving )
         iout << iWARN << "High global exclusion count ("
                       << ((int64)checksum) << " vs "
                       << molecule->numCalcExclusions << "), possible error!\n"
           << iWARN << "This warning is not unusual during minimization.\n"
           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
       else {
         char errmsg[256];
         sprintf(errmsg, "High global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
                 (long long)checksum, (long long)molecule->numCalcExclusions);
         NAMD_bug(errmsg);
       }
     }
 
     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
         ! simParams->goGroPair && ! simParams->qmForcesOn) {
       if ( forgiving || ! simParams->fullElectFrequency ) {
         iout << iWARN << "Low global exclusion count!  ("
           << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
         if ( forgiving ) {
           iout << "\n"
             << iWARN << "This warning is not unusual during minimization.\n"
             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
         } else {
           iout << "  System unstable or pairlistdist or cutoff too small.\n";
         }
         iout << endi;
       } else {
         char errmsg[256];
         sprintf(errmsg, "Low global exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
                 (long long)checksum, (long long)molecule->numCalcExclusions);
         NAMD_bug(errmsg);
       }
     }
 
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
     checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
     if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
          ( ! simParams->mollyOn || step % slowFreq ) ) {
       if ( forgiving )
         iout << iWARN << "High global CUDA exclusion count ("
                       << ((int64)checksum) << " vs "
                       << molecule->numCalcFullExclusions << "), possible error!\n"
           << iWARN << "This warning is not unusual during minimization.\n"
           << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
       else {
         char errmsg[256];
         sprintf(errmsg, "High global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
         NAMD_bug(errmsg);
       }
     }
 
     if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
         ! simParams->goGroPair && ! simParams->qmForcesOn) {
       if ( forgiving || ! simParams->fullElectFrequency ) {
         iout << iWARN << "Low global CUDA exclusion count!  ("
           << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
         if ( forgiving ) {
           iout << "\n"
             << iWARN << "This warning is not unusual during minimization.\n"
             << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
         } else {
           iout << "  System unstable or pairlistdist or cutoff too small.\n";
         }
         iout << endi;
       } else {
         char errmsg[256];
         sprintf(errmsg, "Low global CUDA exclusion count!  (%lld vs %lld)  System unstable or pairlistdist or cutoff too small.\n",
                 (long long)checksum, (long long)molecule->numCalcFullExclusions);
         NAMD_bug(errmsg);
       }
     }
 #endif
 
     checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
     if ( ((int)checksum) && ! marginViolations ) {
       iout << iERROR << "Margin is too small for " << ((int)checksum) <<
         " atoms during timestep " << step << ".\n" << iERROR <<
       "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
     }
     marginViolations += (int)checksum;
 
     checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
     if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
       iout << iINFO <<
         "Pairlistdist is too small for " << ((int)checksum) <<
         " computes during timestep " << step << ".\n" << endi;
     }
     if ( simParams->outputPairlists )  pairlistWarnings += (int)checksum;
 
     checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
     if ( checksum ) {
       if ( forgiving )
         iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
       else NAMD_bug("Stray PME grid charges detected!\n");
     }
 }

BigReal Controller::computeAlchWork ( const int step )

protected

Definition at line 3879 of file Controller.C.

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

Referenced by printEnergies().

                                                   {
   // alchemical scaling factors for groups 1/2 at the previous lambda
   const BigReal oldLambda = simParams->getCurrentLambda(step-1);
   const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
   const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
   const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
   const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
   const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
   const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
   // alchemical scaling factors for groups 1/2 at the new lambda
   const BigReal alchLambda = simParams->getCurrentLambda(step);
   const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
   const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
   const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
   const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
   const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
   const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda); 
 
   return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
           (elec_lambda_12 - elec_lambda_1)*
           (electEnergy_ti_1 + electEnergySlow_ti_1 +  electEnergyPME_ti_1) +
           (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
           (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
           (elec_lambda_22 - elec_lambda_2)*
           (electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2) + 
           (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
   );
 }

void Controller::correctMomentum ( int step )

protected

Definition at line 1255 of file Controller.C.

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

Referenced by integrate().

                                          {
 
     Vector momentum;
     momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
     momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
     momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
     const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
 
     if ( momentum.length2() == 0. )
       iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
     if ( mass == 0. )
       NAMD_die("Total mass is zero in Controller::correctMomentum");
 
     momentum *= (-1./mass);
 
     broadcast->momentumCorrection.publish(step+slowFreq,momentum);
 }

void Controller::cycleBarrier	(	int	doBarrier,
		int	step
	)

protected

Definition at line 4140 of file Controller.C.

References broadcast.

Referenced by integrate().

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

void Controller::enqueueCollections ( int timestep )

protected

Definition at line 3664 of file Controller.C.

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

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

 {
   if ( Output::coordinateNeeded(timestep) ){
     collection->enqueuePositions(timestep,state->lattice);
   }
   if ( Output::velocityNeeded(timestep) )
     collection->enqueueVelocities(timestep);
   if ( Output::forceNeeded(timestep) )
     collection->enqueueForces(timestep);
 }

void Controller::integrate ( int scriptTask )

protected

Definition at line 429 of file Controller.C.

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

Referenced by algorithm().

                                          {
     char traceNote[24];
   
     int step = simParams->firstTimestep;
 
     const int numberOfSteps = simParams->N;
     const int stepsPerCycle = simParams->stepsPerCycle;
 
     const int zeroMomentum = simParams->zeroMomentum;
 
     nbondFreq = simParams->nonbondedFrequency;
     const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
     if (dofull)
       slowFreq = simParams->fullElectFrequency;
     else
       slowFreq = simParams->nonbondedFrequency;
     if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
 
   if ( scriptTask == SCRIPT_RUN ) {
 
     reassignVelocities(step);  // only for full-step velecities
     rescaleaccelMD(step);
     adaptTempUpdate(step); // Init adaptive tempering;
 
     receivePressure(step);
     if ( zeroMomentum && dofull && ! (step % slowFreq) )
                                                 correctMomentum(step);
     printFepMessage(step);
     printTiMessage(step);
     printDynamicsEnergies(step);
     outputFepEnergy(step);
     outputTiEnergy(step);
     if(traceIsOn()){
         traceUserEvent(eventEndOfTimeStep);
         sprintf(traceNote, "s:%d", step);
         traceUserSuppliedNote(traceNote);
     }
     outputExtendedSystem(step);
     rebalanceLoad(step);
 
   }
 
     // Handling SIGINT doesn't seem to be working on Lemieux, and it
     // sometimes causes the net-xxx versions of NAMD to segfault on exit, 
     // so disable it for now.
     // namd_sighandler_t oldhandler = signal(SIGINT, 
     //  (namd_sighandler_t)my_sigint_handler);
     for ( ++step ; step <= numberOfSteps; ++step )
     {
 
         adaptTempUpdate(step);
         rescaleVelocities(step);
         tcoupleVelocities(step);
         stochRescaleVelocities(step);
         berendsenPressure(step);
         langevinPiston1(step);
         rescaleaccelMD(step);
         enqueueCollections(step);  // after lattice scaling!
         receivePressure(step);
         if ( zeroMomentum && dofull && ! (step % slowFreq) )
                                                 correctMomentum(step);
         langevinPiston2(step);
         reassignVelocities(step);
 
       multigratorTemperature(step, 1);
       multigratorPressure(step, 1);
       multigratorPressure(step, 2);
       multigratorTemperature(step, 2);
 
         printDynamicsEnergies(step);
         outputFepEnergy(step);
         outputTiEnergy(step);
         if(traceIsOn()){
             traceUserEvent(eventEndOfTimeStep);
             sprintf(traceNote, "s:%d", step);
             traceUserSuppliedNote(traceNote);
         }
   // if (gotsigint) {
   //   iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
   //   NAMD_quit();
   // }
         outputExtendedSystem(step);
 #if CYCLE_BARRIER
         cycleBarrier(!((step+1) % stepsPerCycle),step);
 #elif  PME_BARRIER
         cycleBarrier(dofull && !(step%slowFreq),step);
 #elif  STEP_BARRIER
         cycleBarrier(1, step);
 #endif
                 
         if(Node::Object()->specialTracing || simParams->statsOn){               
                  int bstep = simParams->traceStartStep;
                  int estep = bstep + simParams->numTraceSteps;          
                  if(step == bstep){
                          traceBarrier(1, step);
                  }
                  if(step == estep){                     
                          traceBarrier(0, step);
                  }
          }
 
 #ifdef MEASURE_NAMD_WITH_PAPI
         if(simParams->papiMeasure) {
                 int bstep = simParams->papiMeasureStartStep;
                 int estep = bstep + simParams->numPapiMeasureSteps;
                 if(step == bstep) {
                         papiMeasureBarrier(1, step);
                 }
                 if(step == estep) {
                         papiMeasureBarrier(0, step);
                 }
         }
 #endif
          
         rebalanceLoad(step);
 
 #if  PME_BARRIER
         cycleBarrier(dofull && !((step+1)%slowFreq),step);   // step before PME
 #endif
     }
     // signal(SIGINT, oldhandler);
 }

void Controller::langevinPiston1 ( int step )

protected

Definition at line 987 of file Controller.C.

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

Referenced by integrate().

 {
   if ( simParams->langevinPistonOn )
   {
     Tensor &strainRate = langevinPiston_strainRate;
     int cellDims = simParams->useFlexibleCell ? 1 : 3;
     BigReal dt = simParams->dt;
     BigReal dt_long = slowFreq * dt;
     BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
     BigReal tau = simParams->langevinPistonPeriod;
     BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
 
 #ifdef DEBUG_PRESSURE
     iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
 #endif
 
     if ( ! ( (step-1) % slowFreq ) )
     {
       BigReal gamma = 1 / simParams->langevinPistonDecay;
       BigReal f1 = exp( -0.5 * dt_long * gamma );
       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
       strainRate *= f1;
       if ( simParams->useFlexibleCell ) {
         // We only use on-diagonal terms (for now)
         if ( simParams->useConstantRatio ) {
           BigReal r = f2 * random->gaussian();
           strainRate.xx += r;
           strainRate.yy += r;
           strainRate.zz += f2 * random->gaussian();
         } else {
           strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
         }
       } else {
         strainRate += f2 * Tensor::identity(random->gaussian());
       }
 
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
 #endif
     }
 
     // Apply surface tension.  If surfaceTensionTarget is zero, we get
     // the default (isotropic pressure) case.
     
     Tensor ptarget;
     ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
     if ( simParams->surfaceTensionTarget != 0. ) {
       ptarget.xx = ptarget.yy = simParams->langevinPistonTarget - 
         simParams->surfaceTensionTarget / state->lattice.c().z;
     }
 
     strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
       ( controlPressure - ptarget );
 
     if (simParams->fixCellDims) {
       if (simParams->fixCellDimX) strainRate.xx = 0;
       if (simParams->fixCellDimY) strainRate.yy = 0;
       if (simParams->fixCellDimZ) strainRate.zz = 0;
     }
 
 
 #ifdef DEBUG_PRESSURE
     iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
 #endif
 
     if ( simParams->langevinPistonBarrier ) {
     if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
     {
       // We only use on-diagonal terms (for now)
       Tensor factor;
       if ( !simParams->useConstantArea ) {
         factor.xx = exp( dt_long * strainRate.xx );
         factor.yy = exp( dt_long * strainRate.yy );
       } else {
         factor.xx = factor.yy = 1;
       }
       factor.zz = exp( dt_long * strainRate.zz );
       broadcast->positionRescaleFactor.publish(step,factor);
       state->lattice.rescale(factor);
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "rescaling by: " << factor << "\n";
 #endif
     }
     } else { // langevinPistonBarrier
     if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
     {
       if ( ! positionRescaleFactor.xx ) {  // first pass without MTS
       // We only use on-diagonal terms (for now)
       Tensor factor;
       if ( !simParams->useConstantArea ) {
         factor.xx = exp( dt_long * strainRate.xx );
         factor.yy = exp( dt_long * strainRate.yy );
       } else {
         factor.xx = factor.yy = 1;
       }
       factor.zz = exp( dt_long * strainRate.zz );
       broadcast->positionRescaleFactor.publish(step,factor);
       positionRescaleFactor = factor;
       strainRate_old = strainRate;
       }
       state->lattice.rescale(positionRescaleFactor);
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "rescaling by: " << positionRescaleFactor << "\n";
 #endif
     }
     if ( ! ( (step-slowFreq/2) % slowFreq ) )
     {
       // We only use on-diagonal terms (for now)
       Tensor factor;
       if ( !simParams->useConstantArea ) {
         factor.xx = exp( (dt+dt_long) * strainRate.xx - dt * strainRate_old.xx );
         factor.yy = exp( (dt+dt_long) * strainRate.yy - dt * strainRate_old.yy );
       } else {
         factor.xx = factor.yy = 1;
       }
       factor.zz = exp( (dt+dt_long) * strainRate.zz - dt * strainRate_old.zz );
       broadcast->positionRescaleFactor.publish(step+1,factor);
       positionRescaleFactor = factor;
       strainRate_old = strainRate;
     }
     }
 
     // corrections to integrator
     if ( ! ( step % nbondFreq ) )
     {
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "correcting strain rate for nbond, ";
 #endif
       strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                 ( (nbondFreq - 1) * controlPressure_nbond );
 #ifdef DEBUG_PRESSURE
       iout << "strain rate: " << strainRate << "\n";
 #endif
     }
     if ( ! ( step % slowFreq ) )
     {
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "correcting strain rate for slow, ";
 #endif
       strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                 ( (slowFreq - 1) * controlPressure_slow );
 #ifdef DEBUG_PRESSURE
       iout << "strain rate: " << strainRate << "\n";
 #endif
     }
     if (simParams->fixCellDims) {
       if (simParams->fixCellDimX) strainRate.xx = 0;
       if (simParams->fixCellDimY) strainRate.yy = 0;
       if (simParams->fixCellDimZ) strainRate.zz = 0;
     }
 
   }
 
 }

void Controller::langevinPiston2 ( int step )

protected

Definition at line 1142 of file Controller.C.

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

Referenced by integrate().

 {
   if ( simParams->langevinPistonOn )
   {
     Tensor &strainRate = langevinPiston_strainRate;
     int cellDims = simParams->useFlexibleCell ? 1 : 3;
     BigReal dt = simParams->dt;
     BigReal dt_long = slowFreq * dt;
     BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
     BigReal tau = simParams->langevinPistonPeriod;
     BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
 
     // corrections to integrator
     if ( ! ( step % nbondFreq ) )
     {
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "correcting strain rate for nbond, ";
 #endif
       strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                 ( (nbondFreq - 1) * controlPressure_nbond );
 #ifdef DEBUG_PRESSURE
       iout << "strain rate: " << strainRate << "\n";
 #endif
     }
     if ( ! ( step % slowFreq ) )
     {
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "correcting strain rate for slow, ";
 #endif
       strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
                 ( (slowFreq - 1) * controlPressure_slow );
 #ifdef DEBUG_PRESSURE
       iout << "strain rate: " << strainRate << "\n";
 #endif
     }
 
     // Apply surface tension.  If surfaceTensionTarget is zero, we get
     // the default (isotropic pressure) case.
    
     Tensor ptarget;
     ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
     if ( simParams->surfaceTensionTarget != 0. ) {
       ptarget.xx = ptarget.yy = simParams->langevinPistonTarget - 
         simParams->surfaceTensionTarget / state->lattice.c().z;
     }
 
     strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
       ( controlPressure - ptarget );
 
  
 #ifdef DEBUG_PRESSURE
     iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
 #endif
 
     if ( ! ( step % slowFreq ) )
     {
       BigReal gamma = 1 / simParams->langevinPistonDecay;
       BigReal f1 = exp( -0.5 * dt_long * gamma );
       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
       strainRate *= f1;
       if ( simParams->useFlexibleCell ) {
         // We only use on-diagonal terms (for now)
         if ( simParams->useConstantRatio ) {
           BigReal r = f2 * random->gaussian();
           strainRate.xx += r;
           strainRate.yy += r;
           strainRate.zz += f2 * random->gaussian();
         } else {
           strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
         }
       } else {
         strainRate += f2 * Tensor::identity(random->gaussian());
       }
 #ifdef DEBUG_PRESSURE
       iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
 #endif
     }
 
 #ifdef DEBUG_PRESSURE
     iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
 #endif
     if (simParams->fixCellDims) {
       if (simParams->fixCellDimX) strainRate.xx = 0;
       if (simParams->fixCellDimY) strainRate.yy = 0;
       if (simParams->fixCellDimZ) strainRate.zz = 0;
     }
   }
 
 
 }

void Controller::minimize ( )

protected

Definition at line 594 of file Controller.C.

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

Referenced by algorithm().

                           {
   // iout << "Controller::minimize() called.\n" << endi;
 
   const int minVerbose = simParams->minVerbose;
   const int numberOfSteps = simParams->N;
   int step = simParams->firstTimestep;
   slowFreq = nbondFreq = 1;
   BigReal linegoal = simParams->minLineGoal;  // 1.0e-3
   const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
 
   CALCULATE
 
   int minSeq = 0;
 
   // just move downhill until initial bad contacts go away or 100 steps
   int old_min_huge_count = 2000000000;
   int steps_at_same_min_huge_count = 0;
   for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
     if ( min_huge_count >= old_min_huge_count ) {
       if ( ++steps_at_same_min_huge_count > 10 ) break;
     } else {
       old_min_huge_count = min_huge_count;
       steps_at_same_min_huge_count = 0;
     }
     iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
     broadcast->minimizeCoefficient.publish(minSeq++,1.);
     enqueueCollections(step);
     CALCULATE
   }
   if ( min_huge_count ) {
     iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
   }
   iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
 
   int atStart = 2;
   int errorFactor = 10;
   BigReal old_f_dot_f = min_f_dot_f;
   broadcast->minimizeCoefficient.publish(minSeq++,0.);
   broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
   int newDir = 1;
   min_f_dot_v = min_f_dot_f;
   min_v_dot_v = min_f_dot_f;
   while ( 1 ) {
     // line minimization
     // bracket minimum on line
     minpoint lo,hi,mid,last;
     BigReal x = 0;
     lo.x = x;
     lo.u = min_energy;
     lo.dudx = -1. * min_f_dot_v;
     lo.noGradient = min_huge_count;
     mid = lo;
     last = mid;
     if ( minVerbose ) {
       iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
       if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
       iout << "\n" << endi;
     }
     BigReal tol = fabs( linegoal * min_f_dot_v );
     iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
             fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
     int start_with_huge = last.noGradient;
     min_reduction->require();
     BigReal maxstep = 0.1 / sqrt(min_reduction->item(0));
     x = maxstep; MOVETO(x);
     // bracket minimum on line
     while ( last.u < mid.u ) {
       if ( last.dudx < mid.dudx * (5.5 - x/maxstep)/5. ) {
         x = 6 * maxstep; break;
       }
       lo = mid; mid = last;
       x += maxstep;
       if ( x > 5.5 * maxstep ) break;
       MOVETO(x)
     }
     if ( x > 5.5 * maxstep ) {
       iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" << endi;
       broadcast->minimizeCoefficient.publish(minSeq++,0.);
       broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
       newDir = 1;
       old_f_dot_f = min_f_dot_f;
       min_f_dot_v = min_f_dot_f;
       min_v_dot_v = min_f_dot_f;
       continue;
     }
     hi = last;
 #define PRINT_BRACKET \
     iout << "LINE MINIMIZER BRACKET: DX " \
          << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
         " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
         lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
     PRINT_BRACKET
     // converge on minimum on line
     int itcnt;
     for ( itcnt = 10; itcnt > 0; --itcnt ) {
       int progress = 1;
       // select new position
       if ( mid.noGradient ) {
        if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) {  // subdivide left side
         x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
         MOVETO(x)
         if ( last.u <= mid.u ) {
           hi = mid; mid = last;
         } else {
           lo = last;
         }
        } else {  // subdivide right side
         x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
         MOVETO(x)
         if ( last.u <= mid.u ) {
           lo = mid; mid = last;
         } else {
           hi = last;
         }
        }
       } else {
        if ( mid.dudx > 0. ) {  // subdivide left side
         BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
         BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
         x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
         BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
         if ( xdiv ) x /= xdiv; else x = last.x;
         if ( x > altxhi ) x = altxhi;
         if ( x < altxlo ) x = altxlo;
         if ( x-last.x == 0 ) break;
         MOVETO(x)
         if ( last.u <= mid.u ) {
           hi = mid; mid = last;
         } else {
           if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
           lo = last;
         }
        } else {  // subdivide right side
         BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
         BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
         x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
         BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
         if ( xdiv ) x /= xdiv; else x = last.x;
         if ( x < altxlo ) x = altxlo;
         if ( x > altxhi ) x = altxhi;
         if ( x-last.x == 0 ) break;
         MOVETO(x)
         if ( last.u <= mid.u ) {
           lo = mid; mid = last;
         } else {
           if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
           hi = last;
         }
        }
       }
       PRINT_BRACKET
       if ( ! progress ) {
         MOVETO(mid.x);
         break;
       }
       if ( fabs(last.dudx) < tol ) break;
       if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
       if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
     }
     // new direction
     broadcast->minimizeCoefficient.publish(minSeq++,0.);
     BigReal c = min_f_dot_f / old_f_dot_f;
     c = ( c > 1.5 ? 1.5 : c );
     if ( atStart ) { c = 0; --atStart; }
     if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
       c = 0;
       if ( errorFactor < 100 ) errorFactor += 10;
     }
     if ( c == 0 ) {
       iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
     }
     broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
     newDir = 1;
     old_f_dot_f = min_f_dot_f;
     min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
     min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
   }
 
 }

BigReal Controller::multigatorCalcEnthalpy	(	BigReal	potentialEnergy,
		int	step,
		int	minimize
	)

protected

Definition at line 919 of file Controller.C.

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

Referenced by printEnergies().

                                                                                           {
   Node *node = Node::Object();
   Molecule *molecule = node->molecule;
   SimParameters *simParameters = node->simParameters;
 
   BigReal V = state->lattice.volume();
   BigReal P0 = simParams->multigratorPressureTarget;
   BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
   BigReal NG = controlNumDegFreedom;
   BigReal Nf = numDegFreedom;
   BigReal tauV = simParams->multigratorPressureRelaxationTime;
   BigReal sumZeta = 0.0;
   for (int i=1;i < simParams->multigratorNoseHooverChainLength;i++) {
     sumZeta += multigratorZeta[i];
   }
   BigReal nuOmegaSum = 0.0;
   for (int i=0;i < simParams->multigratorNoseHooverChainLength;i++) {
     nuOmegaSum += multigratorNu[i]*multigratorNu[i]/(2.0*multigratorOmega[i]);
   }
   BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
   BigReal eta = sqrt(kT0*W)*multigratorXi;
 
   BigReal enthalpy = kineticEnergy + potentialEnergy + eta*eta/(2.0*W) + P0*V + nuOmegaSum + kT0*(Nf*multigratorZeta[0] + sumZeta);
 
 //  if (!(step % 100))
     // fprintf(stderr, "enthalpy %lf %lf %lf %lf %lf %lf %lf\n", enthalpy,
     //   kineticEnergy, potentialEnergy, eta*eta/(2.0*W), P0*V, nuOmegaSum, kT0*(Nf*multigratorZeta[0] + sumZeta));
 
   return enthalpy;
 }

void Controller::multigratorPressure	(	int	step,
		int	callNumber
	)

protected

Definition at line 778 of file Controller.C.

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

Referenced by integrate().

                                                              {
   if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
     BigReal P0 = simParams->multigratorPressureTarget;
     BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
     BigReal NG = controlNumDegFreedom;
     BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
     BigReal tau = simParams->multigratorPressureRelaxationTime;
     BigReal s = 0.5*simParams->multigratorPressureFreq*simParams->dt/tau;
     {
       // Compute new scaling factors and send them to Sequencer
       BigReal V = state->lattice.volume();
       BigReal Pinst = trace(controlPressure)/3.0;
       BigReal PGsum = trace(momentumSqrSum);
       //
       multigratorXiT = multigratorXi + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
       BigReal scale = exp(s*NGfac*multigratorXiT);
       BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
       // fprintf(stderr, "%d | T %lf P %lf V %1.3lf\n", step, temperature, Pinst*PRESSUREFACTOR, V);
       Tensor scaleTensor = Tensor::identity(scale);
       Tensor volScaleTensor = Tensor::identity(scale);
       Tensor velScaleTensor = Tensor::identity(velScale);
       state->lattice.rescale(volScaleTensor);
       if (callNumber == 1) {
         broadcast->positionRescaleFactor.publish(step,scaleTensor);
         broadcast->velocityRescaleTensor.publish(step,velScaleTensor);
       } else {
         broadcast->positionRescaleFactor2.publish(step,scaleTensor);
         broadcast->velocityRescaleTensor2.publish(step,velScaleTensor);      
       }
     }
 
     {
       // Wait here for Sequencer to finish scaling and force computation
       reduction->require();
       Tensor virial_normal;
       Tensor virial_nbond;
       Tensor virial_slow;
       Tensor intVirial_normal;
       Tensor intVirial_nbond;
       Tensor intVirial_slow;
       Vector extForce_normal;
       Vector extForce_nbond;
       Vector extForce_slow;
       GET_TENSOR(momentumSqrSum, reduction, REDUCTION_MOMENTUM_SQUARED);
       GET_TENSOR(virial_normal, reduction, REDUCTION_VIRIAL_NORMAL);
       GET_TENSOR(virial_nbond, reduction, REDUCTION_VIRIAL_NBOND);
       GET_TENSOR(virial_slow, reduction, REDUCTION_VIRIAL_SLOW);
       GET_TENSOR(intVirial_normal, reduction, REDUCTION_INT_VIRIAL_NORMAL);
       GET_TENSOR(intVirial_nbond, reduction, REDUCTION_INT_VIRIAL_NBOND);
       GET_TENSOR(intVirial_slow, reduction, REDUCTION_INT_VIRIAL_SLOW);
       GET_VECTOR(extForce_normal, reduction, REDUCTION_EXT_FORCE_NORMAL);
       GET_VECTOR(extForce_nbond, reduction, REDUCTION_EXT_FORCE_NBOND);
       GET_VECTOR(extForce_slow, reduction, REDUCTION_EXT_FORCE_SLOW);
       calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
         intVirial_normal, intVirial_nbond, intVirial_slow,
         extForce_normal, extForce_nbond, extForce_slow);
       if (callNumber == 2) {
         // Update temperature for the Temperature Cycle that is coming next
         BigReal kineticEnergy = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
         temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
       }
     }
 
     {
       // Update pressure integrator
       BigReal V = state->lattice.volume();
       BigReal Pinst = trace(controlPressure)/3.0;
       BigReal PGsum = trace(momentumSqrSum);
       //
       multigratorXi = multigratorXiT + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
     }
 
   }
 }

void Controller::multigratorTemperature	(	int	step,
		int	callNumber
	)

protected

Definition at line 853 of file Controller.C.

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

Referenced by integrate().

                                                                 {
   if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
     BigReal tau = simParams->multigratorTemperatureRelaxationTime;
     BigReal t = 0.5*simParams->multigratorTemperatureFreq*simParams->dt;
     BigReal kT0 = BOLTZMANN * simParams->multigratorTemperatureTarget;
     BigReal Nf = numDegFreedom;
     int n = simParams->multigratorNoseHooverChainLength;
     BigReal T1, T2, v;
     {
       T1 = temperature;
       BigReal kTinst = BOLTZMANN * temperature;
       for (int i=n-1;i >= 0;i--) {
         if (i == 0) {
           BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
           multigratorNuT[0] = exp(-0.5*t*NuOmega)*multigratorNu[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
         } else if (i == n-1) {
           multigratorNuT[n-1] = multigratorNu[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
         } else {
           BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
           multigratorNuT[i] = exp(-0.5*t*NuOmega)*multigratorNu[i] + 
           0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
         }
       }
       BigReal velScale = exp(-t*multigratorNuT[0]/multigratorOmega[0]);
       v = velScale;
       if (callNumber == 1)
         broadcast->velocityRescaleFactor.publish(step,velScale);
       else
         broadcast->velocityRescaleFactor2.publish(step,velScale);
     }
 
     {
       // Wait here for Sequencer to finish scaling and re-calculating kinetic energy
       multigratorReduction->require();
       BigReal kineticEnergy = multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY);
       temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
       T2 = temperature;
       if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
         // If this is pressure cycle, receive new momentum product
         GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
       }
     }
 
     // fprintf(stderr, "%d | T %lf scale %lf T' %lf\n", step, T1, v, T2);
 
     {
       BigReal kTinst = BOLTZMANN * temperature;
       for (int i=0;i < n;i++) {
         if (i == 0) {
           BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
           multigratorNu[0] = exp(-0.5*t*NuOmega)*multigratorNuT[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
         } else if (i == n-1) {
           multigratorNu[n-1] = multigratorNuT[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
         } else {
           BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
           multigratorNu[i] = exp(-0.5*t*NuOmega)*multigratorNuT[i] + 
           0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
         }
         multigratorZeta[i] += t*multigratorNuT[i]/multigratorOmega[i];
       }
     }
 
   }
 }

void Controller::outputExtendedSystem ( int step )

protected

Definition at line 3968 of file Controller.C.

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

Referenced by algorithm(), and integrate().

 {
 
   if ( step >= 0 ) {
 
     // Write out eXtended System Trajectory (XST) file
     if ( simParams->xstFrequency &&
          ((step % simParams->xstFrequency) == 0) )
     {
       if ( ! xstFile.is_open() )
       {
         iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
         NAMD_backup_file(simParams->xstFilename);
         xstFile.open(simParams->xstFilename);
         while (!xstFile) {
           if ( errno == EINTR ) {
             CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
             xstFile.clear();
             xstFile.open(simParams->xstFilename);
             continue;
           }
           char err_msg[257];
           sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
           NAMD_err(err_msg);
         }
         xstFile << "# NAMD extended system trajectory file" << std::endl;
         writeExtendedSystemLabels(xstFile);
       }
       writeExtendedSystemData(step,xstFile);
       xstFile.flush();
     }
 
     // Write out eXtended System Configuration (XSC) files
     //  Output a restart file
     if ( simParams->restartFrequency &&
          ((step % simParams->restartFrequency) == 0) &&
          (step != simParams->firstTimestep) )
     {
       iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
                 << step << "\n" << endi;
       char fname[140];
       const char *bsuffix = ".old";
       strcpy(fname, simParams->restartFilename);
       if ( simParams->restartSave ) {
         char timestepstr[20];
         sprintf(timestepstr,".%d",step);
         strcat(fname, timestepstr);
         bsuffix = ".BAK";
       }
       strcat(fname, ".xsc");
       NAMD_backup_file(fname,bsuffix);
       ofstream_namd xscFile(fname);
       while (!xscFile) {
         if ( errno == EINTR ) {
           CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
           xscFile.clear();
           xscFile.open(fname);
           continue;
         }
         char err_msg[257];
         sprintf(err_msg, "Error opening XSC restart file %s",fname);
         NAMD_err(err_msg);
       } 
       xscFile << "# NAMD extended system configuration restart file" << std::endl;
       writeExtendedSystemLabels(xscFile);
       writeExtendedSystemData(step,xscFile);
       if (!xscFile) {
         char err_msg[257];
         sprintf(err_msg, "Error writing XSC restart file %s",fname);
         NAMD_err(err_msg);
       } 
     }
 
   }
 
   //  Output final coordinates
   if (step == FILE_OUTPUT || step == END_OF_RUN)
   {
     int realstep = ( step == FILE_OUTPUT ?
         simParams->firstTimestep : simParams->N );
     iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
                 << realstep << "\n" << endi;
     static char fname[140];
     strcpy(fname, simParams->outputFilename);
     strcat(fname, ".xsc");
     NAMD_backup_file(fname);
     ofstream_namd xscFile(fname);
     while (!xscFile) {
       if ( errno == EINTR ) {
         CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
         xscFile.clear();
         xscFile.open(fname);
         continue;
       }
       char err_msg[257];
       sprintf(err_msg, "Error opening XSC output file %s",fname);
       NAMD_err(err_msg);
     } 
     xscFile << "# NAMD extended system configuration output file" << std::endl;
     writeExtendedSystemLabels(xscFile);
     writeExtendedSystemData(realstep,xscFile);
     if (!xscFile) {
       char err_msg[257];
       sprintf(err_msg, "Error writing XSC output file %s",fname);
       NAMD_err(err_msg);
     } 
   }
 
   //  Close trajectory file
   if (step == END_OF_RUN) {
     if ( xstFile.is_open() ) {
       xstFile.close();
       iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
     }
   }
 
 }

void Controller::outputFepEnergy ( int step )

protected

Definition at line 3676 of file Controller.C.

Referenced by integrate().

                                          {
   if (simParams->alchOn && simParams->alchFepOn) {
     const int stepInRun = step - simParams->firstTimestep;
     const int alchEquilSteps = simParams->alchEquilSteps;
     const BigReal alchLambda = simParams->alchLambda;
     const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
 
     if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
       FepNo = 0;
       exp_dE_ByRT = 0.0;
       net_dE = 0.0;
     }
     dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
          electEnergy - electEnergySlow - ljEnergy;
     BigReal RT = BOLTZMANN * simParams->alchTemp;
 
     if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
         (step / simParams->alchIDWSFreq) % 2 == 1 )) {
       // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
       FepNo++;
       exp_dE_ByRT += exp(-dE/RT);
       net_dE += dE;
     }
 
     if (simParams->alchOutFreq) {
       if (stepInRun == 0) {
         if (!fepFile.is_open()) {
           NAMD_backup_file(simParams->alchOutFile);
           fepFile.open(simParams->alchOutFile);
           iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
           if(alchEnsembleAvg){
             fepSum = 0.0;
             fepFile << "#            STEP                 Elec                            "
                     << "vdW                    dE           dE_avg         Temp             dG\n"
                     << "#                           l             l+dl      "
                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
           }
           else{
             fepFile << "#            STEP                 Elec                            "
                     << "vdW                    dE         Temp\n"
                     << "#                           l             l+dl      "
                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
           }
         }
         if(!step){
           fepFile << "#NEW FEP WINDOW: "
                   << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
                   << simParams->alchLambda2;
           if ( simParams->alchLambdaIDWS >= 0. ) {
             fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
           }
           fepFile << std::endl;
         }
       }
       if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
         fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
                 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
                 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
       }
       if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
         writeFepEnergyData(step, fepFile);
         fepFile.flush();
       }
       if (alchEnsembleAvg && (step == simParams->N)) {
         fepSum = fepSum + dG;
         fepFile << "#Free energy change for lambda window [ " << alchLambda
                 << " " << simParams->alchLambda2 << " ] is " << dG
                 << " ; net change until now is " << fepSum << std::endl;
         fepFile.flush();
       }
     }
   }
 }

void Controller::outputTiEnergy ( int step )

protected

Definition at line 3750 of file Controller.C.

Referenced by integrate().

                                         {
  if (simParams->alchOn && simParams->alchThermIntOn) {
   const int stepInRun = step - simParams->firstTimestep;
   const int alchEquilSteps = simParams->alchEquilSteps;
   const int alchLambdaFreq = simParams->alchLambdaFreq;
 
   if (stepInRun == 0 || stepInRun == alchEquilSteps) {
     TiNo = 0;
     net_dEdl_bond_1 = 0;
     net_dEdl_bond_2 = 0;
     net_dEdl_elec_1 = 0;
     net_dEdl_elec_2 = 0;
     net_dEdl_lj_1 = 0;
     net_dEdl_lj_2 = 0;
   }
   TiNo++;
   net_dEdl_bond_1 += bondedEnergy_ti_1;
   net_dEdl_bond_2 += bondedEnergy_ti_2;
   net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1
                       + electEnergyPME_ti_1);
   net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2
                       + electEnergyPME_ti_2);
   net_dEdl_lj_1 += ljEnergy_ti_1;
   net_dEdl_lj_2 += ljEnergy_ti_2;
 
   // Don't accumulate block averages (you'll get a divide by zero!) or even 
   // open the TI output if alchOutFreq is zero.
   if (simParams->alchOutFreq) {
     if (stepInRun == 0 || stepInRun == alchEquilSteps 
         || (! ((step - 1) % simParams->alchOutFreq))) {
       // output of instantaneous dU/dl now replaced with running average
       // over last alchOutFreq steps (except for step 0)
       recent_TiNo = 0;
       recent_dEdl_bond_1 = 0;
       recent_dEdl_bond_2 = 0;
       recent_dEdl_elec_1 = 0;
       recent_dEdl_elec_2 = 0;
       recent_dEdl_lj_1 = 0;
       recent_dEdl_lj_2 = 0;
       recent_alchWork = 0;
     }
     recent_TiNo++;
     recent_dEdl_bond_1 += bondedEnergy_ti_1;
     recent_dEdl_bond_2 += bondedEnergy_ti_2;
     recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 
                            + electEnergyPME_ti_1); 
     recent_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 
                            + electEnergyPME_ti_2); 
     recent_dEdl_lj_1 += ljEnergy_ti_1;
     recent_dEdl_lj_2 += ljEnergy_ti_2;
     recent_alchWork += alchWork;
 
     if (stepInRun == 0) {
       if (!tiFile.is_open()) {
         NAMD_backup_file(simParams->alchOutFile);
         tiFile.open(simParams->alchOutFile);
         /* BKR - This has been rather drastically updated to better match 
            stdout. This was necessary for several reasons:
            1) PME global scaling is obsolete (now removed)
            2) scaling of bonded terms was added
            3) alchemical work is now accumulated when switching is active
          */
         iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
         tiFile << "#TITITLE:    TS";
         tiFile << FORMAT("BOND1");
         tiFile << FORMAT("AVGBOND1");
         tiFile << FORMAT("ELECT1");
         tiFile << FORMAT("AVGELECT1");
         tiFile << "     ";
         tiFile << FORMAT("VDW1");
         tiFile << FORMAT("AVGVDW1");
         tiFile << FORMAT("BOND2");
         tiFile << FORMAT("AVGBOND2");
         tiFile << FORMAT("ELECT2");
         tiFile << "     ";
         tiFile << FORMAT("AVGELECT2");
         tiFile << FORMAT("VDW2");
         tiFile << FORMAT("AVGVDW2");
         if (alchLambdaFreq > 0) {
           tiFile << FORMAT("ALCHWORK");
           tiFile << FORMAT("CUMALCHWORK");
         }
         tiFile << std::endl;
       }
 
       if (alchLambdaFreq > 0) {
         tiFile << "#ALCHEMICAL SWITCHING ACTIVE " 
                << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
                << "\n#LAMBDA SCHEDULE: " 
                << "dL: " << simParams->getLambdaDelta() 
                << " Freq: " << alchLambdaFreq;
       }
       else {
         const BigReal alchLambda = simParams->alchLambda;    
         const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
         const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
         const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
         const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
         const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
         const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
         tiFile << "#NEW TI WINDOW: "
                << "LAMBDA " << alchLambda 
                << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
                << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
                << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
                << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
       }
       tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
              << std::endl;
     }
 
     if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
       tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
              << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
              << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
     }
     if ((step%simParams->alchOutFreq) == 0) {
       writeTiEnergyData(step, tiFile);
       tiFile.flush();
     }
   }
  }
 }

void Controller::printDynamicsEnergies ( int step )

protected

Definition at line 3023 of file Controller.C.

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

Referenced by integrate().

                                                {
 
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     SimParameters *simParameters = node->simParameters;
     Lattice &lattice = state->lattice;
 
     compareChecksums(step);
 
     printEnergies(step,0);
 }

void Controller::printEnergies	(	int	step,
		int	minimize
	)

protected

Definition at line 3035 of file Controller.C.

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

 {
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     SimParameters *simParameters = node->simParameters;
     Lattice &lattice = state->lattice;
 
     // Drude model ANISO energy is added into BOND energy
     // and THOLE energy is added into ELECT energy
 
     BigReal bondEnergy;
     BigReal angleEnergy;
     BigReal dihedralEnergy;
     BigReal improperEnergy;
     BigReal crosstermEnergy;
     BigReal boundaryEnergy;
     BigReal miscEnergy;
     BigReal potentialEnergy;
     BigReal flatEnergy;
     BigReal smoothEnergy;
     BigReal work;
 
     Vector momentum;
     Vector angularMomentum;
     BigReal volume = lattice.volume();
 
     bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
     angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
     dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
     improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
     crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
     boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
     miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
 
     if ( minimize || ! ( step % nbondFreq ) )
     {
       electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
       ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
 
       // JLai
       groLJEnergy = reduction->item(REDUCTION_GRO_LJ_ENERGY);
       groGaussEnergy = reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
 
       goNativeEnergy = reduction->item(REDUCTION_GO_NATIVE_ENERGY);
       goNonnativeEnergy = reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
       goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
 
 //fepb
       bondedEnergyDiff_f = reduction->item(REDUCTION_BONDED_ENERGY_F);
       electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
       ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
 
       bondedEnergy_ti_1 = reduction->item(REDUCTION_BONDED_ENERGY_TI_1);
       electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
       ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
       bondedEnergy_ti_2 = reduction->item(REDUCTION_BONDED_ENERGY_TI_2);
       electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
       ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
 //fepe
     }
 
     if ( minimize || ! ( step % slowFreq ) )
     {
       electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
 //fepb
       electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
 
       electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
       electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
       electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
       electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
 //fepe
     }
 
     if (simParameters->LJcorrection && volume) {
 #ifdef MEM_OPT_VERSION
       NAMD_bug("LJcorrection not supported in memory optimized build.");
 #else
       // Apply tail correction to energy.
       BigReal alchLambda = simParameters->getCurrentLambda(step);
       ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
 
       if (simParameters->alchOn) {
         if (simParameters->alchFepOn) {
           BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
           ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
         } else if (simParameters->alchThermIntOn) {
           ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
           ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
         }
       }
 #endif
     }
 
 //fepb BKR - Compute alchemical work if using dynamic lambda.  This is here so
 //           that the cumulative work can be given during a callback.
     if (simParameters->alchLambdaFreq > 0) {
       if (step <= 
           simParameters->firstTimestep + simParameters->alchEquilSteps) {
         cumAlchWork = 0.0;
       }
       alchWork = computeAlchWork(step);
       cumAlchWork += alchWork;
     }
 //fepe
 
     momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
     momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
     momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
     angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
     angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
     angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
 
     // Ported by JLai
     potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
   + improperEnergy + electEnergy + electEnergySlow + ljEnergy
   + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy 
         + groLJEnergy + groGaussEnergy);
     // End of port
     totalEnergy = potentialEnergy + kineticEnergy;
     flatEnergy = (totalEnergy
         + (1.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
     if ( !(step%slowFreq) ) {
       // only adjust based on most accurate energies
       BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
       if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
       if ( step != simParams->firstTimestep ) {
         smooth2_avg *= 0.9375;
         smooth2_avg += 0.0625 * s;
       }
     }
     smoothEnergy = (flatEnergy + smooth2_avg
         - (4.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
 
     // Reset values for accumulated heat and work.
     if (step <= simParameters->firstTimestep &&
         (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
       heat = 0.0;
       totalEnergy0 = totalEnergy;
     }
     if ( simParameters->outputMomenta && ! minimize &&
          ! ( step % simParameters->outputMomenta ) )
     {
       iout << "MOMENTUM: " << step 
            << " P: " << momentum
            << " L: " << angularMomentum
            << "\n" << endi;
     }
 
     if ( simParameters->outputPressure ) {
       pressure_tavg += pressure;
       groupPressure_tavg += groupPressure;
       tavg_count += 1;
       if ( minimize || ! ( step % simParameters->outputPressure ) ) {
         iout << "PRESSURE: " << step << " "
            << PRESSUREFACTOR * pressure << "\n"
            << "GPRESSURE: " << step << " "
            << PRESSUREFACTOR * groupPressure << "\n";
         if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
            << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
            << "GPRESSAVG: " << step << " "
            << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
         iout << endi;
         pressure_tavg = 0;
         groupPressure_tavg = 0;
         tavg_count = 0;
       }
     }
 
     // pressure profile reductions
     if (pressureProfileSlabs) {
       const int freq = simParams->pressureProfileFreq;
       const int arraysize = 3*pressureProfileSlabs;
       
       BigReal *total = new BigReal[arraysize];
       memset(total, 0, arraysize*sizeof(BigReal));
       const int first = simParams->firstTimestep;
 
       if (ppbonded)    ppbonded->getData(first, step, lattice, total);
       if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
       if (ppint)       ppint->getData(first, step, lattice, total);
       for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
       pressureProfileCount++;
 
       if (!(step % freq)) {
         // convert NAMD internal virial to pressure in units of bar 
         BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
    
         iout << "PRESSUREPROFILE: " << step << " ";
         if (step == first) {
           // output pressure profile for this step
           for (int i=0; i<arraysize; i++) {
             iout << total[i] * scalefac << " ";
           }
         } else {
           // output pressure profile averaged over the last count steps.
           scalefac /= pressureProfileCount;
           for (int i=0; i<arraysize; i++) 
             iout << pressureProfileAverage[i]*scalefac << " ";
         }
         iout << "\n" << endi; 
        
         // Clear the average for the next block
         memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
         pressureProfileCount = 0;
       }
       delete [] total;
     }
   
    if ( step != simParams->firstTimestep || stepInFullRun == 0 ) {  // skip repeated first step
     if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
      int benchPhase = stepInFullRun / simParams->firstLdbStep;
      if ( benchPhase > 0 && benchPhase < 10 ) {
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
       if ( simParams->outputEnergies < 60 ) {
         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
       }
 #endif
       iout << iINFO;
       if ( benchPhase < 4 ) iout << "Initial time: ";
       else iout << "Benchmark time: ";
       iout << CkNumPes() << " CPUs ";
       {
         BigReal wallPerStep =
     (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
   BigReal ns = simParams->dt / 1000000.0;
   BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
   BigReal daysPerNano = wallPerStep * days / ns;
   iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
         iout << memusage_MB() << " MB memory\n" << endi;
       }
      }
      startBenchTime = CmiWallTimer();
     }
     ++stepInFullRun;
    }
 
     printTiming(step);
 
     Vector pairVDWForce, pairElectForce;
     if ( simParameters->pairInteractionOn ) {
       GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
       GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
     }
 
     // Compute cumulative nonequilibrium work (including shadow work).
     if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
       work = totalEnergy - totalEnergy0 - heat;
     }
 
     // callback to Tcl with whatever we can
 #ifdef NAMD_TCL
 #define CALLBACKDATA(LABEL,VALUE) \
                 labels << (LABEL) << " "; values << (VALUE) << " ";
 #define CALLBACKLIST(LABEL,VALUE) \
                 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
     if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
       std::ostringstream labels, values;
 #if CMK_BLUEGENEL
       // the normal version below gives a compiler error
       values.precision(16);
 #else
       values << std::setprecision(16);
 #endif
       CALLBACKDATA("TS",step);
       CALLBACKDATA("BOND",bondEnergy);
       CALLBACKDATA("ANGLE",angleEnergy);
       CALLBACKDATA("DIHED",dihedralEnergy);
       CALLBACKDATA("CROSS",crosstermEnergy);
       CALLBACKDATA("IMPRP",improperEnergy);
       CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
       CALLBACKDATA("VDW",ljEnergy);
       CALLBACKDATA("BOUNDARY",boundaryEnergy);
       CALLBACKDATA("MISC",miscEnergy);
       CALLBACKDATA("POTENTIAL",potentialEnergy);
       CALLBACKDATA("KINETIC",kineticEnergy);
       CALLBACKDATA("TOTAL",totalEnergy);
       CALLBACKDATA("TEMP",temperature);
       CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
       CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
       CALLBACKDATA("VOLUME",lattice.volume());
       CALLBACKLIST("CELL_A",lattice.a());
       CALLBACKLIST("CELL_B",lattice.b());
       CALLBACKLIST("CELL_C",lattice.c());
       CALLBACKLIST("CELL_O",lattice.origin());
       labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
                 << lattice.b_p() << " " << lattice.c_p() << "} ";
       if (simParameters->drudeOn) {
         CALLBACKDATA("DRUDEBOND",drudeBondTemp);
       }
       if ( simParameters->pairInteractionOn ) {
         CALLBACKLIST("VDW_FORCE",pairVDWForce);
         CALLBACKLIST("ELECT_FORCE",pairElectForce);
       }
       if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
         CALLBACKLIST("HEAT",heat);
         CALLBACKLIST("WORK",work);
       }
       if (simParameters->alchOn) {
         if (simParameters->alchThermIntOn) {
           CALLBACKLIST("BOND1", bondedEnergy_ti_1);
           CALLBACKLIST("ELEC1", electEnergy_ti_1 + electEnergySlow_ti_1 +
                                 electEnergyPME_ti_1);
           CALLBACKLIST("VDW1", ljEnergy_ti_1);
           CALLBACKLIST("BOND2", bondedEnergy_ti_2);
           CALLBACKLIST("ELEC2", electEnergy_ti_2 + electEnergySlow_ti_2 +
                                 electEnergyPME_ti_2);
           CALLBACKLIST("VDW2", ljEnergy_ti_2);
           if (simParameters->alchLambdaFreq > 0) {
             CALLBACKLIST("CUMALCHWORK", cumAlchWork);
           }
         } else if (simParameters->alchFepOn) {
           CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
                                 improperEnergy + bondedEnergyDiff_f);
           CALLBACKLIST("ELEC2", electEnergy_f + electEnergySlow_f);
           CALLBACKLIST("VDW2", ljEnergy_f);
         } 
       }
 
       labels << '\0';  values << '\0';  // insane but makes Linux work
       state->callback_labelstring = labels.str();
       state->callback_valuestring = values.str();
       // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
     }
 #undef CALLBACKDATA
 #endif
 
     drudeBondTempAvg += drudeBondTemp;
 
     temp_avg += temperature;
     pressure_avg += trace(pressure)/3.;
     groupPressure_avg += trace(groupPressure)/3.;
     avg_count += 1;
 
     if ( simParams->outputPairlists && pairlistWarnings &&
                                 ! (step % simParams->outputPairlists) ) {
       iout << iINFO << pairlistWarnings <<
         " pairlist warnings in past " << simParams->outputPairlists <<
         " steps.\n" << endi;
       pairlistWarnings = 0;
     }
 
     BigReal enthalpy;
     if (simParameters->multigratorOn && ((step % simParameters->outputEnergies) == 0)) {
       enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
     }
 
     // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
     if ( ! minimize &&  step % simParameters->outputEnergies ) return;
     // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
 
     if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
       IMDEnergies energies;
       energies.tstep = step;
       energies.T = temp_avg/avg_count;
       energies.Etot = totalEnergy;
       energies.Epot = potentialEnergy;
       energies.Evdw = ljEnergy;
       energies.Eelec = electEnergy + electEnergySlow;
       energies.Ebond = bondEnergy;
       energies.Eangle = angleEnergy;
       energies.Edihe = dihedralEnergy + crosstermEnergy;
       energies.Eimpr = improperEnergy;
       Node::Object()->imd->gather_energies(&energies);
     }
 
     if ( marginViolations ) {
       iout << iERROR << marginViolations <<
         " margin violations detected since previous energy output.\n" << endi;
     }
     marginViolations = 0;
 
     int precision = simParameters->outputEnergiesPrecision;
     if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
     {
         iout << "ETITLE:      TS";
         iout << FORMAT("BOND", precision);
         iout << FORMAT("ANGLE", precision);
         iout << FORMAT("DIHED", precision);
         if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS", precision);
         iout << FORMAT("IMPRP", precision);
               iout << "     ";
         iout << FORMAT("ELECT", precision);
         iout << FORMAT("VDW", precision);
         iout << FORMAT("BOUNDARY", precision);
         iout << FORMAT("MISC", precision);
         iout << FORMAT("KINETIC", precision);
               iout << "     ";
         iout << FORMAT("TOTAL", precision);
         iout << FORMAT("TEMP", precision);
         iout << FORMAT("POTENTIAL", precision);
         // iout << FORMAT("TOTAL2", precision);
         iout << FORMAT("TOTAL3", precision);
         iout << FORMAT("TEMPAVG", precision);
         if ( volume != 0. ) {
                 iout << "     ";
           iout << FORMAT("PRESSURE", precision);
           iout << FORMAT("GPRESSURE", precision);
           iout << FORMAT("VOLUME", precision);
           iout << FORMAT("PRESSAVG", precision);
           iout << FORMAT("GPRESSAVG", precision);
         }
         if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
           iout << "     ";
           iout << FORMAT("HEAT", precision);
           iout << FORMAT("WORK", precision);
         }
         if (simParameters->drudeOn) {
           iout << "     ";
           iout << FORMAT("DRUDEBOND", precision);
           iout << FORMAT("DRBONDAVG", precision);
         }
         // Ported by JLai
         if (simParameters->goGroPair) {
           iout << "     ";
           iout << FORMAT("GRO_PAIR_LJ", precision);
           iout << FORMAT("GRO_PAIR_GAUSS", precision);
         }
 
         if (simParameters->goForcesOn) {
           iout << "     ";
           iout << FORMAT("NATIVE", precision);
           iout << FORMAT("NONNATIVE", precision);
           //iout << FORMAT("REL_NATIVE", precision);
           //iout << FORMAT("REL_NONNATIVE", precision);
           iout << FORMAT("GOTOTAL", precision);
           //iout << FORMAT("GOAVG", precision);
         }
         // End of port -- JLai
 
         if (simParameters->alchOn) {
           if (simParameters->alchThermIntOn) {
             iout << "\nTITITLE:     TS";
             iout << FORMAT("BOND1", precision);
             iout << FORMAT("ELECT1", precision);
             iout << FORMAT("VDW1", precision);
             iout << FORMAT("BOND2", precision);
             iout << "     ";
             iout << FORMAT("ELECT2", precision);
             iout << FORMAT("VDW2", precision);
             if (simParameters->alchLambdaFreq > 0) {
               iout << FORMAT("LAMBDA", precision);
               iout << FORMAT("ALCHWORK", precision);
               iout << FORMAT("CUMALCHWORK", precision);
             }
           } else if (simParameters->alchFepOn) {
             iout << "\nFEPTITLE:    TS";
             iout << FORMAT("BOND2", precision);
             iout << FORMAT("ELECT2", precision);
             iout << FORMAT("VDW2", precision);
             if (simParameters->alchLambdaFreq > 0) {
               iout << FORMAT("LAMBDA", precision);
             }
           }
         }
 
         iout << "\n\n" << endi;
         
         if (simParameters->qmForcesOn) {
             iout << "QMETITLE:      TS";
             iout << FORMAT("QMID", precision);
             iout << FORMAT("ENERGY", precision);
             if (simParameters->PMEOn) iout << FORMAT("PMECORRENERGY", precision);
             iout << "\n\n" << endi;
         }
         
     }
 
     // N.B.  HP's aCC compiler merges FORMAT calls in the same expression.
     //       Need separate statements because data returned in static array.
     iout << ETITLE(step);
     iout << FORMAT(bondEnergy, precision);
     iout << FORMAT(angleEnergy, precision);
     if ( simParameters->mergeCrossterms ) {
       iout << FORMAT(dihedralEnergy+crosstermEnergy, precision);
     } else {
       iout << FORMAT(dihedralEnergy, precision);
       iout << FORMAT(crosstermEnergy, precision);
     }
     iout << FORMAT(improperEnergy, precision);
     iout << "     ";
     iout << FORMAT(electEnergy+electEnergySlow, precision);
     iout << FORMAT(ljEnergy, precision);
     iout << FORMAT(boundaryEnergy, precision);
     iout << FORMAT(miscEnergy, precision);
     iout << FORMAT(kineticEnergy, precision);
     iout << "     ";
     iout << FORMAT(totalEnergy, precision);
     iout << FORMAT(temperature, precision);
     iout << FORMAT(potentialEnergy, precision);
     // iout << FORMAT(flatEnergy, precision);
     iout << FORMAT(smoothEnergy, precision);
     iout << FORMAT(temp_avg/avg_count, precision);
     if ( volume != 0. )
     {
         iout << "     ";
         iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3., precision);
         iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3., precision);
         iout << FORMAT(volume, precision);
         iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count, precision);
         iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count, precision);
     }
     if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
           iout << "     ";
           iout << FORMAT(heat, precision);
           iout << FORMAT(work, precision);
     }
     if (simParameters->drudeOn) {
         iout << "     ";
         iout << FORMAT(drudeBondTemp, precision);
         iout << FORMAT(drudeBondTempAvg/avg_count, precision);
     }
     // Ported by JLai
     if (simParameters->goGroPair) {
       iout << "     ";
       iout << FORMAT(groLJEnergy, precision);
       iout << FORMAT(groGaussEnergy, precision);
     }
 
     if (simParameters->goForcesOn) {
       iout << "     ";
       iout << FORMAT(goNativeEnergy, precision);
       iout << FORMAT(goNonnativeEnergy, precision);
       //iout << FORMAT(relgoNativeEnergy, precision);
       //iout << FORMAT(relgoNonnativeEnergy, precision);
       iout << FORMAT(goTotalEnergy, precision);
       //iout << FORMAT("not implemented", precision);
     } // End of port -- JLai
 
     if (simParameters->alchOn) { 
       if (simParameters->alchThermIntOn) {
         iout << "\n";
         iout << TITITLE(step);
         iout << FORMAT(bondedEnergy_ti_1, precision);
         iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 + 
                        electEnergyPME_ti_1, precision);
         iout << FORMAT(ljEnergy_ti_1, precision);
         iout << FORMAT(bondedEnergy_ti_2, precision);
         iout << "     ";
         iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 +
                        electEnergyPME_ti_2, precision);
         iout << FORMAT(ljEnergy_ti_2, precision);
         if (simParameters->alchLambdaFreq > 0) {
           iout << FORMAT(simParameters->getCurrentLambda(step), precision);
           iout << FORMAT(alchWork, precision);
           iout << FORMAT(cumAlchWork, precision);
         }
       } else if (simParameters->alchFepOn) {
         iout << "\n";
         iout << FEPTITLE2(step);
         iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy 
                        + improperEnergy + bondedEnergyDiff_f, precision);
         iout << FORMAT(electEnergy_f + electEnergySlow_f, precision);
         iout << FORMAT(ljEnergy_f, precision);
         if (simParameters->alchLambdaFreq > 0) {
           iout << FORMAT(simParameters->getCurrentLambda(step), precision);
         }
       }
     }
 
     iout << "\n\n" << endi;
 
 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
      char webout[80];
      sprintf(webout,"%d %d %d %d",(int)totalEnergy,
              (int)(potentialEnergy),
              (int)kineticEnergy,(int)temperature);
      CApplicationDepositNode0Data(webout);
 #endif
 
     if (simParameters->pairInteractionOn) {
       iout << "PAIR INTERACTION:";
       iout << " STEP: " << step;
       iout << " VDW_FORCE: ";
       iout << FORMAT(pairVDWForce.x, precision);
       iout << FORMAT(pairVDWForce.y, precision);
       iout << FORMAT(pairVDWForce.z, precision);
       iout << " ELECT_FORCE: ";
       iout << FORMAT(pairElectForce.x, precision);
       iout << FORMAT(pairElectForce.y, precision);
       iout << FORMAT(pairElectForce.z, precision);
       iout << "\n" << endi;
     }
     drudeBondTempAvg = 0;
     temp_avg = 0;
     pressure_avg = 0;
     groupPressure_avg = 0;
     avg_count = 0;
 
     if ( fflush_count ) {
       --fflush_count;
       fflush(stdout);
     }
 }

void Controller::printFepMessage ( int step )

protected

Definition at line 1274 of file Controller.C.

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

Referenced by integrate().

 {
   if (simParams->alchOn && simParams->alchFepOn 
       && !simParams->alchLambdaFreq) {
     const BigReal alchLambda = simParams->alchLambda;
     const BigReal alchLambda2 = simParams->alchLambda2;
     const BigReal alchLambdaIDWS = simParams->alchLambdaIDWS;
     const BigReal alchTemp = simParams->alchTemp;
     const int alchEquilSteps = simParams->alchEquilSteps;
     iout << "FEP: RESETTING FOR NEW FEP WINDOW "
          << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2;
     if ( alchLambdaIDWS >= 0. ) {
       iout << " LAMBDA_IDWS " << alchLambdaIDWS;
     }
     iout << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
          << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
          << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp 
          << " K FOR FEP CALCULATION\n" << endi;
   }
 } 

void Controller::printMinimizeEnergies ( int step )

protected

Definition at line 3004 of file Controller.C.

References compareChecksums(), RequireReduction::item(), min_energy, min_f_dot_f, 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.

                                                {
 
     rescaleaccelMD(step,1);
     receivePressure(step,1);
 
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     compareChecksums(step,1);
 
     printEnergies(step,1);
 
     min_energy = totalEnergy;
     min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
     min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
     min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
     min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
 
 }

void Controller::printTiMessage ( int step )

protected

Definition at line 1294 of file Controller.C.

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

Referenced by integrate().

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

void Controller::printTiming ( int step )

protected

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

                                      {
 
     if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
     {
       const double endWTime = CmiWallTimer() - firstWTime;
       const double endCTime = CmiTimer() - firstCTime;
 
       // fflush about once per minute
       if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
 
       const double elapsedW = 
         (endWTime - startWTime) / simParams->outputTiming;
       const double elapsedC = 
         (endCTime - startCTime) / simParams->outputTiming;
 
       const double remainingW = elapsedW * (simParams->N - step);
       const double remainingW_hours = remainingW / 3600;
 
       startWTime = endWTime;
       startCTime = endCTime;
 
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
       if ( simParams->outputEnergies < 60 &&
            step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
         iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
       }
 #endif
       if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
         CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
                   ", %g hours remaining, %f MB of memory in use.\n",
                   step, endCTime, elapsedC, endWTime, elapsedW,
                   remainingW_hours, memusage_MB());
         if ( fflush_count ) { --fflush_count; fflush(stdout); }
       }
     }
 }

void Controller::reassignVelocities ( int step )

protected

Definition at line 1308 of file Controller.C.

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

Referenced by integrate().

 {
   const int reassignFreq = simParams->reassignFreq;
   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
     BigReal newTemp = simParams->reassignTemp;
     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
     if ( simParams->reassignIncr > 0.0 ) {
       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
         newTemp = simParams->reassignHold;
     } else {
       if ( newTemp < simParams->reassignHold )
         newTemp = simParams->reassignHold;
     }
     iout << "REASSIGNING VELOCITIES AT STEP " << step
          << " TO " << newTemp << " KELVIN.\n" << endi;
   }
 }

void Controller::rebalanceLoad ( int step )

protected

Definition at line 4126 of file Controller.C.

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

Referenced by integrate().

 {
   if ( ! ldbSteps ) { 
     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
   }
   if ( ! --ldbSteps ) {
     startBenchTime -= CmiWallTimer();
         Node::Object()->outputPatchComputeMaps("before_ldb", step);
     LdbCoordinator::Object()->rebalance(this);  
         startBenchTime += CmiWallTimer();
     fflush_count = 3;
   }
 }

void Controller::receivePressure	(	int	step,
		int	minimize = `0`
	)

protected

Definition at line 1420 of file Controller.C.

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

Referenced by integrate(), and printMinimizeEnergies().

 {
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     SimParameters *simParameters = node->simParameters;
     Lattice &lattice = state->lattice;
 
     reduction->require();
 
     Tensor virial_normal;
     Tensor virial_nbond;
     Tensor virial_slow;
 #ifdef ALTVIRIAL
     Tensor altVirial_normal;
     Tensor altVirial_nbond;
     Tensor altVirial_slow;
 #endif
     Tensor intVirial_normal;
     Tensor intVirial_nbond;
     Tensor intVirial_slow;
     Vector extForce_normal;
     Vector extForce_nbond;
     Vector extForce_slow;
 
 #if 1
     numDegFreedom = molecule->num_deg_freedom();
     int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
     int numFixedGroups = molecule->num_fixed_groups();
     int numFixedAtoms = molecule->num_fixed_atoms();
 #endif
 #if 0
     int numAtoms = molecule->numAtoms;
     numDegFreedom = 3 * numAtoms;
     int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
     int numFixedAtoms =
         ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
     int numLonepairs = molecule->numLonepairs;
     int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
     if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
     if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
     if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
     if ( ! ( numFixedAtoms || molecule->numConstraints
         || simParameters->comMove || simParameters->langevinOn ) ) {
       numDegFreedom -= 3;
       numGroupDegFreedom -= 3;
     }
     if (simParameters->pairInteractionOn) {
       // this doesn't attempt to deal with fixed atoms or constraints
       numDegFreedom = 3 * molecule->numFepInitial;
     }
     int numRigidBonds = molecule->numRigidBonds;
     int numFixedRigidBonds =
         ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
 
     // numLonepairs is subtracted here because all lonepairs have a rigid bond
     // to oxygen, but all of the LP degrees of freedom are dealt with above
     numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
 #endif
 
     kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
     kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
 
     BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
         reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
     BigReal groupKineticEnergyCentered = kineticEnergyCentered -
         reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
 
     BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
                                         / ( numDegFreedom * BOLTZMANN );
     BigReal atomTempCentered = 2.0 * kineticEnergyCentered
                                         / ( numDegFreedom * BOLTZMANN );
     BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
                                         / ( numGroupDegFreedom * BOLTZMANN );
     BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
                                         / ( numGroupDegFreedom * BOLTZMANN );
 
     /*  test code for comparing different temperatures  
     iout << "TEMPTEST: " << step << " " << 
         atomTempHalfstep << " " <<
         atomTempCentered << " " <<
         groupTempHalfstep << " " <<
         groupTempCentered << "\n" << endi;
   iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
     numGroupDegFreedom << "\n" << endi;
      */
 
     GET_TENSOR(momentumSqrSum,reduction,REDUCTION_MOMENTUM_SQUARED);
 
     GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
     GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
     GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
 
 #ifdef ALTVIRIAL
     GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
     GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
     GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
 #endif
 
     GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
     GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
     GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
 
     GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
     GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
     GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
     // APH NOTE: These four lines are now done in calcPressure()
     // Vector extPosition = lattice.origin();
     // virial_normal -= outer(extForce_normal,extPosition);
     // virial_nbond -= outer(extForce_nbond,extPosition);
     // virial_slow -= outer(extForce_slow,extPosition);
 
     kineticEnergy = kineticEnergyCentered;
     temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );
 
     if (simParameters->drudeOn) {
       BigReal drudeComKE
         = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
       BigReal drudeBondKE
         = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
       int g_bond = 3 * molecule->numDrudeAtoms;
       int g_com = numDegFreedom - g_bond;
       temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
       drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
     }
 
     // Calculate number of degrees of freedom (controlNumDegFreedom)
     if ( simParameters->useGroupPressure )
     {
       controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
       if ( ! ( numFixedAtoms || molecule->numConstraints
   || simParameters->comMove || simParameters->langevinOn ) ) {
         controlNumDegFreedom -= 1;
       }
     }
     else
     {
       controlNumDegFreedom = numDegFreedom / 3;
     }
     if (simParameters->fixCellDims) {
       if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
       if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
       if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
     }
 
     // Calculate pressure tensors using virials
     calcPressure(step, minimize,
       virial_normal, virial_nbond, virial_slow,
       intVirial_normal, intVirial_nbond, intVirial_slow,
       extForce_normal, extForce_nbond, extForce_slow);
 
 #ifdef DEBUG_PRESSURE
     iout << iINFO << "Control pressure = " << controlPressure <<
       " = " << controlPressure_normal << " + " <<
       controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
 #endif
     if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
         iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
              << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
              << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
              << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
              << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
              << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
              << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
              << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
              << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
              << endi;
    }
 }

void Controller::recvCheckpointAck ( checkpoint & cp )

protected

Definition at line 4117 of file Controller.C.

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

Referenced by Node::recvCheckpointAck().

                                                  {  // initiating replica
   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
     state->lattice = cp.lattice;
     *(ControllerState*)this = cp.state;
   }
   CkpvAccess(_qd)->process();
 }

void Controller::recvCheckpointReq	(	const char *	key,
		int	task,
		checkpoint &	cp
	)

protected

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

                                                                             {  // responding replica
   switch ( task ) {
     case SCRIPT_CHECKPOINT_STORE:
       if ( ! checkpoints.count(key) ) {
         checkpoints[key] = new checkpoint;
       }
       *checkpoints[key] = cp;
       break;
     case SCRIPT_CHECKPOINT_LOAD:
       if ( ! checkpoints.count(key) ) {
         NAMD_die("Unable to load checkpoint, requested key was never stored.");
       }
       cp = *checkpoints[key];
       break;
     case SCRIPT_CHECKPOINT_SWAP:
       if ( ! checkpoints.count(key) ) {
         NAMD_die("Unable to swap checkpoint, requested key was never stored.");
       }
       std::swap(cp,*checkpoints[key]);
       break;
     case SCRIPT_CHECKPOINT_FREE:
       if ( ! checkpoints.count(key) ) {
         NAMD_die("Unable to free checkpoint, requested key was never stored.");
       }
       delete checkpoints[key];
       checkpoints.erase(key);
       break;
   }
 }

void Controller::rescaleaccelMD	(	int	step,
		int	minimize = `0`
	)

protected

Definition at line 1839 of file Controller.C.

Referenced by integrate(), and printMinimizeEnergies().

 {
     if ( !simParams->accelMDOn ) return;
 
     amd_reduction->require();
 
     // copy all to submit_reduction
     for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
       submit_reduction->item(i) += amd_reduction->item(i);
     }
     submit_reduction->submit();
 
     if (step == simParams->firstTimestep) {
         accelMDdVAverage = 0;
         accelMDdV = 0;
     }
 //    if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
     if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
 
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     Lattice &lattice = state->lattice;
 
     const BigReal accelMDE = simParams->accelMDE;
     const BigReal accelMDalpha = simParams->accelMDalpha;
     const BigReal accelMDTE = simParams->accelMDTE;
     const BigReal accelMDTalpha = simParams->accelMDTalpha;
     const int accelMDOutFreq = simParams->accelMDOutFreq;
 
     //GaMD
     static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
     static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
     static BigReal k0P, kP, EP;
     static BigReal k0D, kD, ED;
     static int V_n = 1, iEusedD, iEusedP;
     static char warnD, warnP;
     static int restfreq;
 
     const int firststep = (simParams->accelMDFirstStep > simParams->firstTimestep) ? simParams->accelMDFirstStep : simParams->firstTimestep;
     const int ntcmdprep = (simParams->accelMDGcMDPrepSteps > 0) ? simParams->accelMDFirstStep + simParams->accelMDGcMDPrepSteps : 0;
     const int ntcmd     = simParams->accelMDFirstStep + simParams->accelMDGcMDSteps;
     const int ntebprep  = ntcmd + simParams->accelMDGEquiPrepSteps;
     const int nteb      = ntcmd + simParams->accelMDGEquiSteps;
     const int ntave     = simParams->accelMDGStatWindow;
 
     BigReal volume;
     BigReal bondEnergy;
     BigReal angleEnergy;
     BigReal dihedralEnergy;
     BigReal improperEnergy;
     BigReal crosstermEnergy;
     BigReal boundaryEnergy;
     BigReal miscEnergy;
     BigReal amd_electEnergy;
     BigReal amd_ljEnergy;
     BigReal amd_ljEnergy_Corr = 0.;
     BigReal amd_electEnergySlow;
     BigReal amd_groLJEnergy;
     BigReal amd_groGaussEnergy;
     BigReal amd_goTotalEnergy;
     BigReal potentialEnergy;
     BigReal factor_dihe = 1;
     BigReal factor_tot = 1;
     BigReal testV;
     BigReal dV = 0.;
     Vector  accelMDfactor;
     Tensor vir; //auto initialized to 0
     Tensor vir_dihe;
     Tensor vir_normal;
     Tensor vir_nbond;
     Tensor vir_slow;
 
     volume = lattice.volume();
 
     bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
     angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
     dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
     improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
     crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
     boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
     miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);
 
     GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
     GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
     GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
     GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
 
     if ( !( step % nbondFreq ) ) {
       amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
       amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
       amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
       amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
       BigReal goNativeEnergy = amd_reduction->item(REDUCTION_GO_NATIVE_ENERGY);
       BigReal goNonnativeEnergy = amd_reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
       amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
     } else {
       amd_electEnergy = electEnergy;
       amd_ljEnergy = ljEnergy;
       amd_groLJEnergy = groLJEnergy;
       amd_groGaussEnergy = groGaussEnergy;
       amd_goTotalEnergy = goTotalEnergy;
     }
 
     if( simParams->LJcorrection && volume ) {
         // Obtain tail correction (copied from printEnergies())
         // This value is only printed for sanity check
         // accelMD currently does not 'boost' tail correction
         amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
     }
 
     if ( !( step % slowFreq ) ) {
       amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
     } else {
       amd_electEnergySlow = electEnergySlow;
     }
 
     potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
         improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
         crosstermEnergy + boundaryEnergy + miscEnergy +
         amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
 
     //GaMD
     if(simParams->accelMDG){
        // if first time running accelMDG module
        if(step == firststep) {
           iEusedD = iEusedP = simParams->accelMDGiE;
           warnD   = warnP   = '\0';
 
           //restart file reading
           if(simParams->accelMDGRestart){
              FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
              char line[256];
              int dihe_n=0, tot_n=0;
              if(!rest){
                 sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
                 NAMD_die(line);
              }
 
              while(fgets(line, 256, rest)){
                 if(line[0] == 'D'){
                    dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
                                    &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
                 }
                 else if(line[0] == 'T'){
                    tot_n  = sscanf(line+1, " %d %la %la %la %la %la %la %la",
                                    &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
                 }
              }
 
              fclose(rest);
 
              V_n++;
 
              //check dihe settings
              if(simParams->accelMDdihe || simParams->accelMDdual){
                 if(dihe_n !=8)
                    NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
                 k0D = kD * (VmaxD - VminD);
                 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
                    << " Vmax " << VmaxD << " Vmin " << VminD 
                    << " Vavg " << VavgD << " sigmaV " << sigmaVD
                    << " E " << ED << " k " << kD
                    << "\n" << endi;
              }
 
              //check tot settings
              if(!simParams->accelMDdihe || simParams->accelMDdual){
                 if(tot_n !=8)
                    NAMD_die("Invalid total potential energy format in accelMDG restart file");
                 k0P = kP * (VmaxP - VminP);
                 iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
                    << " Vmax " << VmaxP << " Vmin " << VminP 
                    << " Vavg " << VavgP << " sigmaV " << sigmaVP
                    << " E " << EP << " k " << kP
                    << "\n" << endi;
             }
 
             iEusedD = (ED == VmaxD) ? 1 : 2;
             iEusedP = (EP == VmaxP) ? 1 : 2;
           }
           //local restfreq follows NAMD restartfreq (default: 500)
           restfreq = simParams->restartFrequency ? simParams->restartFrequency : 500;
        }
 
        //for dihedral potential
        if(simParams->accelMDdihe || simParams->accelMDdual){
           testV = dihedralEnergy + crosstermEnergy;
 
           //write restart file every restartfreq steps
           if(step > firststep && step % restfreq == 0)
              write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
                    true, false);
           //write restart file at the end of the simulation
           if( simParams->accelMDLastStep > 0 ){
               if( step == simParams->accelMDLastStep )
                   write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
                           true, true);
           }
           else if(step == simParams->N)
               write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD, 
                       true, true);
 
           //conventional MD
           if(step < ntcmd){
              //very first step
              if(step == firststep && !simParams->accelMDGRestart){
                 //initialize stat
                 VmaxD = VminD = VavgD = testV;
                 M2D = sigmaVD = 0.;
              }
              //first step after cmdprep
              else if(step == ntcmdprep && ntcmdprep != 0){
                 //reset stat
                 VmaxD = VminD = VavgD = testV;
                 M2D = sigmaVD = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
              }
              //every ntave steps
              else if(ntave > 0 && step % ntave == 0){
                 //update Vmax Vmin
                 if(testV > VmaxD) VmaxD = testV;
                 if(testV < VminD) VminD = testV;
                 //reset avg and std
                 VavgD = testV;
                 M2D = sigmaVD = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
              }
              //normal steps
              else
                 calc_accelMDG_mean_std(testV, V_n,
                       &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
 
              //last cmd step
              if(step == ntcmd - 1){
                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D, 
                       VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
              }
           }
           //equilibration
           else if(step < nteb){
              calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
                    &dV, &factor_dihe, &vir);
 
              //first step after cmd
              if(step == ntcmd && simParams->accelMDGresetVaftercmd){
                 //reset stat
                 VmaxD = VminD = VavgD = testV;
                 M2D = sigmaVD = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
              }
              //every ntave steps
              else if(ntave > 0 && step % ntave == 0){
                 //update Vmax Vmin
                 if(testV > VmaxD) VmaxD = testV;
                 if(testV < VminD) VminD = testV;
                 //reset avg and std
                 VavgD = testV;
                 M2D = sigmaVD = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
              }
              else
                 calc_accelMDG_mean_std(testV, V_n, 
                       &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
 
              //steps after ebprep
              if(step >= ntebprep)
                 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
                     calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0D,
                           VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
           }
           //production
           else{
              calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
                    &dV, &factor_dihe, &vir);
           }
 
        }
        //for total potential
        if(!simParams->accelMDdihe || simParams->accelMDdual){
           Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
           testV = potentialEnergy;
           if(simParams->accelMDdual){
              testV -= dihedralEnergy + crosstermEnergy;
              vir_tot -= vir_dihe;
           }
 
           //write restart file every restartfreq steps
           if(step > firststep && step % restfreq == 0)
              write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
                    !simParams->accelMDdual, false);
           //write restart file at the end of simulation
           if( simParams->accelMDLastStep > 0 ){
               if( step == simParams->accelMDLastStep )
                   write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
                           !simParams->accelMDdual, true);
           }
           else if(step == simParams->N)
              write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP, 
                    !simParams->accelMDdual, true);
 
           //conventional MD
           if(step < ntcmd){
              //very first step
              if(step == firststep && !simParams->accelMDGRestart){
                 //initialize stat
                 VmaxP = VminP = VavgP = testV;
                 M2P = sigmaVP = 0.;
              }
              //first step after cmdprep
              else if(step == ntcmdprep && ntcmdprep != 0){
                 //reset stat
                 VmaxP = VminP = VavgP = testV;
                 M2P = sigmaVP = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
              }
              //every ntave steps
              else if(ntave > 0 && step % ntave == 0){
                 //update Vmax Vmin
                 if(testV > VmaxP) VmaxP = testV;
                 if(testV < VminP) VminP = testV;
                 //reset avg and std
                 VavgP = testV;
                 M2P = sigmaVP = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
              }
              //normal steps
              else
                 calc_accelMDG_mean_std(testV, V_n,
                       &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
              //last cmd step
              if(step == ntcmd - 1){
                 calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P, 
                       VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
              }
           }
           //equilibration
           else if(step < nteb){
              calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
                    &dV, &factor_tot, &vir);
 
              //first step after cmd
              if(step == ntcmd && simParams->accelMDGresetVaftercmd){
                 //reset stat
                 VmaxP = VminP = VavgP = testV;
                 M2P = sigmaVP = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
              }
              //every ntave steps
              else if(ntave > 0 && step % ntave == 0){
                 //update Vmax Vmin
                 if(testV > VmaxP) VmaxP = testV;
                 if(testV < VminP) VminP = testV;
                 //reset avg and std
                 VavgP = testV;
                 M2P = sigmaVP = 0.;
                 iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
              }
              else
                 calc_accelMDG_mean_std(testV, V_n,
                       &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
 
              //steps after ebprep
              if(step >= ntebprep)
                 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
                     calc_accelMDG_E_k(simParams->accelMDGiE, step, simParams->accelMDGSigma0P,
                           VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
           }
           //production
           else{
              calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
                    &dV, &factor_tot, &vir);
           }
 
        }
        accelMDdVAverage += dV;
 
        //first step after ntcmdprep
        if((ntcmdprep > 0 && step == ntcmdprep) || 
           (step < nteb && ntave > 0 && step % ntave == 0) ||
           (simParams->accelMDGresetVaftercmd && step == ntcmd)){
           V_n = 1;
        }
 
        if(step < nteb)
            V_n++;
 
     }
     //aMD
     else{
         if (simParams->accelMDdihe) {
 
             testV = dihedralEnergy + crosstermEnergy;
             if ( testV < accelMDE ) {
                 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
                 factor_dihe *= factor_dihe;
                 vir = vir_dihe * (factor_dihe - 1.0);
                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
                 accelMDdVAverage += dV;
             }  
 
         } else if (simParams->accelMDdual) {
 
             testV = dihedralEnergy + crosstermEnergy;
             if ( testV < accelMDE ) {
                 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
                 factor_dihe *= factor_dihe;
                 vir = vir_dihe * (factor_dihe - 1.0) ;
                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
             }
 
             testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
             if ( testV < accelMDTE ) {
                 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
                 factor_tot *= factor_tot;
                 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
                 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
             }
             accelMDdVAverage += dV;
 
         } else {
 
             testV = potentialEnergy;
             if ( testV < accelMDE ) {
                 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
                 factor_tot *= factor_tot;
                 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
                 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
                 accelMDdVAverage += dV;
             }
         } 
     }
 
     accelMDfactor[0]=factor_dihe;
     accelMDfactor[1]=factor_tot;
     accelMDfactor[2]=1;
     broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
     virial_amd = vir;
     accelMDdV = dV;
 
     if ( factor_tot < 0.001 ) {
        iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
             << "\n" << endi;
     }
 
     if ( ! (step % accelMDOutFreq) ) {
         if ( !(step == simParams->firstTimestep) ) {
              accelMDdVAverage = accelMDdVAverage/accelMDOutFreq; 
         }
         iout << "ACCELERATED MD: STEP " << step
              << " dV "   << dV
              << " dVAVG " << accelMDdVAverage 
              << " BOND " << bondEnergy
              << " ANGLE " << angleEnergy
              << " DIHED " << dihedralEnergy+crosstermEnergy
              << " IMPRP " << improperEnergy
              << " ELECT " << amd_electEnergy+amd_electEnergySlow
              << " VDW "  << amd_ljEnergy
              << " POTENTIAL "  << potentialEnergy;
         if(amd_ljEnergy_Corr)
             iout << " LJcorr " << amd_ljEnergy_Corr;
         iout << "\n" << endi;
         if(simParams->accelMDG){
             if(simParams->accelMDdihe || simParams->accelMDdual)
                 iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD 
                     << " Vmax " << VmaxD << " Vmin " << VminD 
                     << " Vavg " << VavgD << " sigmaV " << sigmaVD
                     << " E " << ED << " k0 " << k0D << " k " << kD
                     << "\n" << endi;
             if(warnD & 1)
                 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
             if(warnD & 2)
                 iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
             if(!simParams->accelMDdihe || simParams->accelMDdual)
                 iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP 
                     << " Vmax " << VmaxP << " Vmin " << VminP 
                     << " Vavg " << VavgP << " sigmaV " << sigmaVP
                     << " E " << EP << " k0 " << k0P << " k " << kP
                     << "\n" << endi;
             if(warnP & 1)
                 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
             if(warnP & 2)
                 iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
                     << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
             warnD = warnP = '\0';
         }
 
         accelMDdVAverage = 0;
 
         if (simParams->accelMDDebugOn) {
            Tensor p_normal;
            Tensor p_nbond;
            Tensor p_slow;
            Tensor p;
            if ( volume != 0. )  {
                  p_normal = vir_normal/volume;
                  p_nbond  = vir_nbond/volume;
                  p_slow = vir_slow/volume;
                  p = vir/volume;
            }    
            iout << " accelMD Scaling Factor: " << accelMDfactor << "\n" 
                 << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
                 << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
                 << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
                 << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n" 
                 << endi;
         }
    }
 }

void Controller::rescaleVelocities ( int step )

protected

Definition at line 1233 of file Controller.C.

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

Referenced by integrate().

 {
   const int rescaleFreq = simParams->rescaleFreq;
   if ( rescaleFreq > 0 ) {
     rescaleVelocities_sumTemps += temperature;  ++rescaleVelocities_numTemps;
     if ( rescaleVelocities_numTemps == rescaleFreq ) {
       BigReal avgTemp = rescaleVelocities_sumTemps / rescaleVelocities_numTemps;
       BigReal rescaleTemp = simParams->rescaleTemp;
       if ( simParams->adaptTempOn && simParams->adaptTempRescale && (step > simParams->adaptTempFirstStep ) && 
                 (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
         rescaleTemp = adaptTempT;
       }
       BigReal factor = sqrt(rescaleTemp/avgTemp);
       broadcast->velocityRescaleFactor.publish(step,factor);
       //iout << "RESCALING VELOCITIES AT STEP " << step
       //     << " FROM AVERAGE TEMPERATURE OF " << avgTemp
       //     << " TO " << rescaleTemp << " KELVIN.\n" << endi;
       rescaleVelocities_sumTemps = 0;  rescaleVelocities_numTemps = 0;
     }
   }
 }

void Controller::resumeAfterTraceBarrier ( int step )

Definition at line 4157 of file Controller.C.

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

Referenced by Node::resumeAfterTraceBarrier().

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

void Controller::run ( void )

Definition at line 268 of file Controller.C.

References awaken(), CTRL_STK_SZ, and DebugM.

Referenced by NamdState::runController().

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

double Controller::stochRescaleCoefficient ( )

protected

Calculate new coefficient for stochastic velocity rescaling and update heat.

Definition at line 1364 of file Controller.C.

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

Referenced by stochRescaleVelocities().

                                            {
   const double stochRescaleTemp = simParams->stochRescaleTemp;
   double coefficient = 1;
   if ( temperature > 0 ) {
     double R1 = random->gaussian();
     // double gammaShape = 0.5*(numDegFreedom - 1);
     // R2sum is the sum of (numDegFreedom - 1) squared normal variables,
     // which is chi-squared distributed.
     // This is in turn a special case of the Gamma distribution,
     // which converges to a normal distribution in the limit of a
     // large shape parameter.
     // double R2sum = 2*(gammaShape+sqrt(gammaShape)*random->gaussian())+R1*R1;
     double R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
     double tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
 
     coefficient = sqrt(stochRescaleTimefactor +
         (1 - stochRescaleTimefactor)*tempfactor*R2sum +
         2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*
           stochRescaleTimefactor)); 
   }
   heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient-1);
   return coefficient;
 }

void Controller::stochRescaleVelocities ( int step )

protected

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

Generate and broadcast the scale factor for stochastic velocity rescaling.

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

Parameters

step	the current timestep

Definition at line 1348 of file Controller.C.

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

Referenced by integrate().

 {
   if ( simParams->stochRescaleOn ) {
     ++stochRescale_count;
     if ( stochRescale_count == simParams->stochRescaleFreq ) {
       double coefficient = stochRescaleCoefficient();
       broadcast->stochRescaleCoefficient.publish(step,coefficient);
       stochRescale_count = 0;
     }
   }
 }

void Controller::tcoupleVelocities ( int step )

protected

Definition at line 1326 of file Controller.C.

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

Referenced by integrate().

 {
   if ( simParams->tCoupleOn )
   {
     const BigReal tCoupleTemp = simParams->tCoupleTemp;
     BigReal coefficient = 1.;
     if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
     broadcast->tcoupleCoefficient.publish(step,coefficient);
   }
 }

void Controller::terminate ( void )

protected

Definition at line 4178 of file Controller.C.

References BackEnd::awaken().

Referenced by algorithm().

                                {
   BackEnd::awaken();
   CthFree(thread);
   CthSuspend();
 }

void Controller::traceBarrier	(	int	turnOnTrace,
		int	step
	)

protected

Definition at line 4150 of file Controller.C.

Referenced by integrate().

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

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

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

Referenced by rescaleaccelMD().

                                                                                                                                                              {
    FILE *rest;
    char timestepstr[20];
    char *filename;
    int baselen;
    char *restart_name;
    const char *bsuffix;
 
    if(lasttime || simParams->restartFrequency == 0){
            filename = simParams->outputFilename;
            bsuffix = ".BAK";
    }
    else{
            filename = simParams->restartFilename;
            bsuffix = ".old";
    }
 
    baselen = strlen(filename);
    restart_name = new char[baselen+26];
 
    strcpy(restart_name, filename);
    if ( simParams->restartSave && !lasttime) {
       sprintf(timestepstr,".%d",step_n);
       strcat(restart_name, timestepstr);
       bsuffix = ".BAK";
    }
    strcat(restart_name, ".gamd");
 
    if(write_topic){
       NAMD_backup_file(restart_name,bsuffix);
 
       rest = fopen(restart_name, "w");
       if(rest){
          fprintf(rest, "# NAMD accelMDG restart file\n"
                "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
                "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
                type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
          fclose(rest);
          iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
       }
       else
          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
    }
    else{
       rest = fopen(restart_name, "a");
       if(rest){
          fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
                           type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
          fclose(rest);
       }
       else
          iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
    }
 
    delete[] restart_name;
 }

void Controller::writeExtendedSystemData	(	int	step,
		ofstream_namd &	file
	)

protected

Definition at line 3643 of file Controller.C.

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

Referenced by outputExtendedSystem().

                                                                       {
   Lattice &lattice = state->lattice;
   file.precision(12);
   file << step;
     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
   if ( simParams->langevinPistonOn ) {
     Vector strainRate = diagonal(langevinPiston_strainRate);
     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
     file << " " << strainRate.x;
     file << " " << strainRate.y;
     file << " " << strainRate.z;
     file << " " << strainRate2.x;
     file << " " << strainRate2.y;
     file << " " << strainRate2.z;
   }
   file << std::endl;
 }

void Controller::writeExtendedSystemLabels ( ofstream_namd & file )

protected

Definition at line 3630 of file Controller.C.

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

Referenced by outputExtendedSystem().

                                                               {
   Lattice &lattice = state->lattice;
   file << "#$LABELS step";
   if ( lattice.a_p() ) file << " a_x a_y a_z";
   if ( lattice.b_p() ) file << " b_x b_y b_z";
   if ( lattice.c_p() ) file << " c_x c_y c_z";
   file << " o_x o_y o_z";
   if ( simParams->langevinPistonOn ) {
     file << " s_x s_y s_z s_u s_v s_w";
   }
   file << std::endl;
 }

void Controller::writeFepEnergyData	(	int	step,
		ofstream_namd &	file
	)

protected

Definition at line 3908 of file Controller.C.

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

Referenced by outputFepEnergy().

                                                                  {
   BigReal eeng = electEnergy + electEnergySlow;
   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
   BigReal RT = BOLTZMANN * simParams->alchTemp;
 
   const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
   const int stepInRun = step - simParams->firstTimestep;
   const BigReal alchLambda = simParams->alchLambda;
 
   if(stepInRun){
     if ( simParams->alchLambdaIDWS >= 0. &&
         (step / simParams->alchIDWSFreq) % 2 == 0 ) {
       // IDWS is active and we are on a "backward" timestep
       fepFile << FEPTITLE_BACK(step);
     } else {
       // "forward" timestep
       fepFile << FEPTITLE(step);
     }
     fepFile << FORMAT(eeng);
     fepFile << FORMAT(eeng_f);
     fepFile << FORMAT(ljEnergy);
     fepFile << FORMAT(ljEnergy_f);
     fepFile << FORMAT(dE);
     if(alchEnsembleAvg){
       BigReal dE_avg = net_dE / FepNo;
       fepFile << FORMAT(dE_avg);
     }
     fepFile << FORMAT(temperature);
     if(alchEnsembleAvg){
       dG = -(RT * log(exp_dE_ByRT / FepNo));
       fepFile << FORMAT(dG);
     }
     fepFile << std::endl;
   }
 }

void Controller::writeTiEnergyData	(	int	step,
		ofstream_namd &	file
	)

protected

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

                                                                 {
   tiFile << TITITLE(step);
   tiFile << FORMAT(recent_dEdl_bond_1 / recent_TiNo);
   tiFile << FORMAT(net_dEdl_bond_1 / TiNo);
   tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
   tiFile << FORMAT(net_dEdl_elec_1/TiNo);
   tiFile << "     ";
   tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
   tiFile << FORMAT(net_dEdl_lj_1/TiNo);
   tiFile << FORMAT(recent_dEdl_bond_2 / recent_TiNo);
   tiFile << FORMAT(net_dEdl_bond_2 / TiNo);
   tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
   tiFile << "     ";
   tiFile << FORMAT(net_dEdl_elec_2/TiNo);
   tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
   tiFile << FORMAT(net_dEdl_lj_2/TiNo);
   if (simParams->alchLambdaFreq > 0) {
     tiFile << FORMAT(recent_alchWork / recent_TiNo);
     tiFile << FORMAT(cumAlchWork);
   }
   tiFile << std::endl;
 }

Friends And Related Function Documentation

friend class CheckpointMsg

friend

Definition at line 55 of file Controller.h.

friend class Node

friend

Definition at line 54 of file Controller.h.

friend class ScriptTcl

friend

Definition at line 53 of file Controller.h.

Member Data Documentation

BigReal Controller::accelMDdV

Definition at line 50 of file Controller.h.

Referenced by rescaleaccelMD().

BigReal Controller::accelMDdVAverage

protected

Definition at line 301 of file Controller.h.

Referenced by rescaleaccelMD().

Bool Controller::adaptTempAutoDt

protected

Definition at line 324 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal Controller::adaptTempBetaMax

protected

Definition at line 318 of file Controller.h.

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

BigReal Controller::adaptTempBetaMin

protected

Definition at line 317 of file Controller.h.

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

BigReal* Controller::adaptTempBetaN

protected

Definition at line 313 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

int Controller::adaptTempBin

protected

Definition at line 319 of file Controller.h.

Referenced by adaptTempUpdate().

int Controller::adaptTempBins

protected

Definition at line 320 of file Controller.h.

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

BigReal Controller::adaptTempCg

protected

Definition at line 322 of file Controller.h.

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

BigReal Controller::adaptTempDBeta

protected

Definition at line 321 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal Controller::adaptTempDt

protected

Definition at line 323 of file Controller.h.

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

BigReal Controller::adaptTempDTave

protected

Definition at line 315 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal Controller::adaptTempDTavenum

protected

Definition at line 316 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal Controller::adaptTempDtMax

protected

Definition at line 326 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal Controller::adaptTempDtMin

protected

Definition at line 325 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

BigReal* Controller::adaptTempPotEnergyAve

protected

Definition at line 310 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyAveDen

protected

Definition at line 308 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyAveNum

protected

Definition at line 307 of file Controller.h.

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

int* Controller::adaptTempPotEnergySamples

protected

Definition at line 312 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyVar

protected

Definition at line 311 of file Controller.h.

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

BigReal* Controller::adaptTempPotEnergyVarNum

protected

Definition at line 309 of file Controller.h.

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

ofstream_namd Controller::adaptTempRestartFile

protected

Definition at line 327 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempWriteRestart().

BigReal Controller::adaptTempT

protected

Definition at line 314 of file Controller.h.

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

BigReal Controller::alchWork

protected

Definition at line 151 of file Controller.h.

Referenced by outputTiEnergy(), and printEnergies().

RequireReduction* Controller::amd_reduction

protected

Definition at line 238 of file Controller.h.

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

int Controller::avg_count

protected

Definition at line 84 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::bondedEnergy_ti_1

protected

Definition at line 126 of file Controller.h.

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

BigReal Controller::bondedEnergy_ti_2

protected

Definition at line 127 of file Controller.h.

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

BigReal Controller::bondedEnergyDiff_f

protected

Definition at line 113 of file Controller.h.

Referenced by outputFepEnergy(), and printEnergies().

ControllerBroadcasts* Controller::broadcast

protected

Definition at line 251 of file Controller.h.

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

Lattice Controller::checkpoint_lattice

protected

Definition at line 269 of file Controller.h.

Referenced by algorithm().

ControllerState Controller::checkpoint_state

protected

Definition at line 270 of file Controller.h.

Referenced by algorithm().

int Controller::checkpoint_stored

protected

Definition at line 268 of file Controller.h.

Referenced by algorithm(), and Controller().

int Controller::checkpoint_task

protected

Definition at line 277 of file Controller.h.

Referenced by recvCheckpointAck().

std::map<std::string,checkpoint*> Controller::checkpoints

protected

Definition at line 276 of file Controller.h.

Referenced by algorithm(), and recvCheckpointReq().

CollectionMaster* const Controller::collection

protected

Definition at line 249 of file Controller.h.

Referenced by enqueueCollections().

int Controller::computeChecksum

protected

Definition at line 89 of file Controller.h.

Referenced by compareChecksums().

int Controller::controlNumDegFreedom

protected

Definition at line 169 of file Controller.h.

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

Tensor Controller::controlPressure

protected

Definition at line 170 of file Controller.h.

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

Tensor Controller::controlPressure_nbond

protected

Definition at line 77 of file Controller.h.

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

Tensor Controller::controlPressure_normal

protected

Definition at line 76 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::controlPressure_slow

protected

Definition at line 78 of file Controller.h.

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

BigReal Controller::cumAlchWork

protected

Definition at line 140 of file Controller.h.

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

BigReal Controller::dE

protected

Definition at line 119 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::dG

protected

Definition at line 121 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::drudeBondTemp

protected

Definition at line 155 of file Controller.h.

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

BigReal Controller::drudeBondTempAvg

protected

Definition at line 156 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::electEnergy

protected

Definition at line 104 of file Controller.h.

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

BigReal Controller::electEnergy_f

protected

Definition at line 114 of file Controller.h.

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

BigReal Controller::electEnergy_ti_1

protected

Definition at line 128 of file Controller.h.

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

BigReal Controller::electEnergy_ti_2

protected

Definition at line 131 of file Controller.h.

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

BigReal Controller::electEnergyPME_ti_1

protected

Definition at line 141 of file Controller.h.

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

BigReal Controller::electEnergyPME_ti_2

protected

Definition at line 142 of file Controller.h.

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

BigReal Controller::electEnergySlow

protected

Definition at line 105 of file Controller.h.

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

BigReal Controller::electEnergySlow_f

protected

Definition at line 115 of file Controller.h.

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

BigReal Controller::electEnergySlow_ti_1

protected

Definition at line 129 of file Controller.h.

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

BigReal Controller::electEnergySlow_ti_2

protected

Definition at line 132 of file Controller.h.

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

BigReal Controller::exp_dE_ByRT

protected

Definition at line 118 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

ofstream_namd Controller::fepFile

protected

Definition at line 258 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

int Controller::FepNo

protected

Definition at line 122 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::fepSum

protected

Definition at line 124 of file Controller.h.

Referenced by outputFepEnergy().

int Controller::fflush_count

protected

Definition at line 222 of file Controller.h.

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

BigReal Controller::goNativeEnergy

protected

Definition at line 109 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::goNonnativeEnergy

protected

Definition at line 110 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::goTotalEnergy

protected

Definition at line 111 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::groGaussEnergy

protected

Definition at line 108 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

BigReal Controller::groLJEnergy

protected

Definition at line 107 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

Tensor Controller::groupPressure

protected

Definition at line 168 of file Controller.h.

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

BigReal Controller::groupPressure_avg

protected

Definition at line 83 of file Controller.h.

Referenced by Controller(), and printEnergies().

Tensor Controller::groupPressure_nbond

protected

Definition at line 74 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_normal

protected

Definition at line 73 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_slow

protected

Definition at line 75 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::groupPressure_tavg

protected

Definition at line 86 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::heat

protected

heat exchanged with the thermostat since firstTimestep

Definition at line 162 of file Controller.h.

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

BigReal Controller::kineticEnergy

protected

Definition at line 158 of file Controller.h.

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

BigReal Controller::kineticEnergyCentered

protected

Definition at line 160 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

BigReal Controller::kineticEnergyHalfstep

protected

Definition at line 159 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

Tensor Controller::langevinPiston_origStrainRate

protected

Definition at line 204 of file Controller.h.

Referenced by Controller().

int Controller::ldbSteps

protected

Definition at line 220 of file Controller.h.

Referenced by rebalanceLoad().

BigReal Controller::ljEnergy

protected

Definition at line 106 of file Controller.h.

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

BigReal Controller::ljEnergy_f

protected

Definition at line 116 of file Controller.h.

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

BigReal Controller::ljEnergy_f_left

protected

Definition at line 117 of file Controller.h.

BigReal Controller::ljEnergy_ti_1

protected

Definition at line 130 of file Controller.h.

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

BigReal Controller::ljEnergy_ti_2

protected

Definition at line 133 of file Controller.h.

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

int Controller::marginViolations

protected

Definition at line 90 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

BigReal Controller::min_energy

protected

Definition at line 94 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

BigReal Controller::min_f_dot_f

protected

Definition at line 95 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

BigReal Controller::min_f_dot_v

protected

Definition at line 96 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

int Controller::min_huge_count

protected

Definition at line 98 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

RequireReduction* Controller::min_reduction

protected

Definition at line 60 of file Controller.h.

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

BigReal Controller::min_v_dot_v

protected

Definition at line 97 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

Tensor Controller::momentumSqrSum

protected

Definition at line 211 of file Controller.h.

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

std::vector<BigReal> Controller::multigratorNu

protected

Definition at line 213 of file Controller.h.

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

std::vector<BigReal> Controller::multigratorNuT

protected

Definition at line 214 of file Controller.h.

Referenced by Controller(), and multigratorTemperature().

std::vector<BigReal> Controller::multigratorOmega

protected

Definition at line 215 of file Controller.h.

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

RequireReduction* Controller::multigratorReduction

protected

Definition at line 217 of file Controller.h.

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

BigReal Controller::multigratorXi

protected

Definition at line 209 of file Controller.h.

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

BigReal Controller::multigratorXiT

protected

Definition at line 210 of file Controller.h.

Referenced by multigratorPressure().

std::vector<BigReal> Controller::multigratorZeta

protected

Definition at line 216 of file Controller.h.

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

int Controller::nbondFreq

protected

Definition at line 79 of file Controller.h.

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

BigReal Controller::net_dE

protected

Definition at line 120 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

BigReal Controller::net_dEdl_bond_1

protected

Definition at line 134 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_bond_2

protected

Definition at line 135 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_elec_1

protected

Definition at line 136 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_elec_2

protected

Definition at line 137 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_lj_1

protected

Definition at line 138 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::net_dEdl_lj_2

protected

Definition at line 139 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int64_t Controller::numDegFreedom

protected

Definition at line 101 of file Controller.h.

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

Lattice Controller::origLattice

protected

Definition at line 281 of file Controller.h.

Referenced by Controller().

int Controller::pairlistWarnings

protected

Definition at line 91 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

Tensor Controller::positionRescaleFactor

protected

Definition at line 206 of file Controller.h.

Referenced by langevinPiston1().

PressureProfileReduction* Controller::ppbonded

protected

Definition at line 242 of file Controller.h.

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

PressureProfileReduction* Controller::ppint

protected

Definition at line 244 of file Controller.h.

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

PressureProfileReduction* Controller::ppnonbonded

protected

Definition at line 243 of file Controller.h.

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

Tensor Controller::pressure

protected

Definition at line 167 of file Controller.h.

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

Tensor Controller::pressure_amd

protected

Definition at line 71 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

BigReal Controller::pressure_avg

protected

Definition at line 82 of file Controller.h.

Referenced by Controller(), and printEnergies().

Tensor Controller::pressure_nbond

protected

Definition at line 69 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_normal

protected

Definition at line 68 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_slow

protected

Definition at line 70 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

Tensor Controller::pressure_tavg

protected

Definition at line 85 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal* Controller::pressureProfileAverage

protected

Definition at line 247 of file Controller.h.

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

int Controller::pressureProfileCount

protected

Definition at line 246 of file Controller.h.

Referenced by Controller(), and printEnergies().

int Controller::pressureProfileSlabs

protected

Definition at line 245 of file Controller.h.

Referenced by Controller(), and printEnergies().

Random* Controller::random

protected

Definition at line 234 of file Controller.h.

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

BigReal Controller::recent_alchWork

protected

Definition at line 150 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_bond_1

protected

Definition at line 144 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_bond_2

protected

Definition at line 145 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_elec_1

protected

Definition at line 146 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_elec_2

protected

Definition at line 147 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_lj_1

protected

Definition at line 148 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::recent_dEdl_lj_2

protected

Definition at line 149 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::recent_TiNo

protected

Definition at line 152 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

RequireReduction* Controller::reduction

protected

Definition at line 237 of file Controller.h.

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

int Controller::rescaleVelocities_numTemps

protected

Definition at line 175 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

BigReal Controller::rescaleVelocities_sumTemps

protected

Definition at line 174 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

SimParameters* const Controller::simParams

protected

int Controller::slowFreq

protected

Definition at line 80 of file Controller.h.

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

BigReal Controller::smooth2_avg2

protected

Definition at line 166 of file Controller.h.

NamdState* const Controller::state

protected

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

Referenced by printEnergies().

int Controller::stochRescale_count

protected

Count time steps until next stochastic velocity rescaling.

Definition at line 192 of file Controller.h.

Referenced by Controller(), and stochRescaleVelocities().

BigReal Controller::stochRescaleTimefactor

protected

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

Definition at line 195 of file Controller.h.

Referenced by Controller(), and stochRescaleCoefficient().

Tensor Controller::strainRate_old

protected

Definition at line 205 of file Controller.h.

Referenced by langevinPiston1().

SubmitReduction* Controller::submit_reduction

protected

Definition at line 239 of file Controller.h.

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

int Controller::tavg_count

protected

Definition at line 87 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::temp_avg

protected

Definition at line 81 of file Controller.h.

Referenced by Controller(), and printEnergies().

BigReal Controller::temperature

protected

Definition at line 161 of file Controller.h.

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

ofstream_namd Controller::tiFile

protected

Definition at line 262 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

int Controller::TiNo

protected

Definition at line 143 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

BigReal Controller::totalEnergy

protected

Definition at line 103 of file Controller.h.

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

BigReal Controller::totalEnergy0

protected

totalEnergy at firstTimestep

Definition at line 164 of file Controller.h.

Referenced by Controller(), and printEnergies().

Tensor Controller::virial_amd

protected

Definition at line 72 of file Controller.h.

Referenced by calcPressure(), and rescaleaccelMD().

ofstream_namd Controller::xstFile

protected

Definition at line 252 of file Controller.h.

Referenced by outputExtendedSystem().

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

Classes

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation