#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 (

Classes

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation