#include <Controller.h>

Inheritance diagram for Controller:

Classes
struct	checkpoint

struct	CthThreadWrapper

Public Member Functions
	Controller (NamdState *s)

virtual	~Controller (void)

void	run (void)

void	awaken (void)

void	resumeAfterTraceBarrier (int)

BigReal	getTIderivative (void) const

void	resetMovingAverage ()

void	stochRescaleVelocities (int)

double	stochRescaleCoefficient ()

Public Attributes
BigReal	accelMDdV

int	stochRescale_count

BigReal	stochRescaleTimefactor

Protected Types
enum	mc_axis_pick { MC_X = 0, MC_Y, MC_Z, MC_AXIS_TOTAL }

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)

BigReal	getTotalPotentialEnergy (int step)

void	printTiming (int)

void	printMinimizeEnergies (int)

void	printDynamicsEnergies (int)

void	printEnergies (int step, int minimize)

void	printFepMessage (int)

void	printTiMessage (int)

void	enqueueCollections (int)

void	correctMomentum (int step)

void	rescaleVelocities (int)

void	reassignVelocities (int)

void	tcoupleVelocities (int)

void	berendsenPressure (int)

void	langevinPiston1 (int)

void	langevinPiston2 (int)

void	monteCarloPressure_prepare (int)

void	monteCarloPressure_accept (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	suspend (void)

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	ljEnergySlow

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

Tensor	langevinPiston_origStrainRate

Tensor	strainRate_old

Tensor	positionRescaleFactor

int	mc_trial [MC_AXIS_TOTAL]

int	mc_accept [MC_AXIS_TOTAL]

int	mc_totalTry

int	mc_totalAccept

int	mc_picked_axis

BigReal	mc_totalEnergyOld

Lattice	mc_oldLattice

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 *	reductionBasic

RequireReduction *	reductionGpuResident

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

RunningAverage	perfstats

MovingAverage	totalEnergyAverage

MovingAverage	temperatureAverage

MovingAverage	pressureAverage

MovingAverage	groupPressureAverage

MovingAverage	pressureAverage_xx

MovingAverage	pressureAverage_yx

MovingAverage	pressureAverage_yy

MovingAverage	pressureAverage_zx

MovingAverage	pressureAverage_zy

MovingAverage	pressureAverage_zz

MovingAverage	groupPressureAverage_xx

MovingAverage	groupPressureAverage_xy

MovingAverage	groupPressureAverage_xz

MovingAverage	groupPressureAverage_yx

MovingAverage	groupPressureAverage_yy

MovingAverage	groupPressureAverage_yz

MovingAverage	groupPressureAverage_zx

MovingAverage	groupPressureAverage_zy

MovingAverage	groupPressureAverage_zz

Protected Attributes inherited from ControllerState
Tensor	langevinPiston_strainRate

Tensor	berendsenPressure_avg

Vector	monteCarloMaxVolume

int	berendsenPressure_count

BigReal	smooth2_avg

Friends
class	ScriptTcl

class	Node

class	CheckpointMsg

Detailed Description

Definition at line 105 of file Controller.h.

Member Enumeration Documentation

◆ mc_axis_pick

enum Controller::mc_axis_pick

protected

Enumerator
MC_X
MC_Y
MC_Z
MC_AXIS_TOTAL

Definition at line 302 of file Controller.h.

                       {
       MC_X = 0, // for x-axis, constant ratio, or isotropic fluctuation
       MC_Y, // for y-axis
       MC_Z, // for z-axis or constant area fluctuations
       MC_AXIS_TOTAL
     };

Constructor & Destructor Documentation

◆ Controller()

Controller::Controller ( NamdState * s )

Definition at line 174 of file Controller.C.

References SimParameters::accelMDOn, amd_reduction, avg_count, AVGXY, ControllerState::berendsenPressure_avg, ControllerState::berendsenPressure_count, BOLTZMANN, bondedEnergy_ti_1, bondedEnergy_ti_2, broadcast, checkpoint_stored, SimParameters::CUDASOAintegrateMode, cumAlchWork, drudeBondTemp, drudeBondTempAvg, SimParameters::dt, electEnergy_ti_1, electEnergy_ti_2, electEnergySlow_ti_1, electEnergySlow_ti_2, groupPressure_avg, groupPressure_tavg, heat, Tensor::identity(), langevinPiston_origStrainRate, ControllerState::langevinPiston_strainRate, ljEnergy_ti_1, ljEnergy_ti_2, mc_accept, MC_AXIS_TOTAL, mc_totalAccept, mc_totalTry, mc_trial, min_reduction, Node::molecule, ControllerState::monteCarloMaxVolume, SimParameters::monteCarloMaxVolume, SimParameters::monteCarloPressureOn, MULTIGRATOR_REDUCTION_MAX_RESERVED, SimParameters::multigratorNoseHooverChainLength, multigratorNu, multigratorNuT, multigratorOmega, SimParameters::multigratorOn, multigratorReduction, SimParameters::multigratorTemperatureRelaxationTime, SimParameters::multigratorTemperatureTarget, multigratorXi, multigratorZeta, PatchData::nodeBroadcast, Molecule::num_deg_freedom(), 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, reductionBasic, reductionGpuResident, REDUCTIONS_AMD, REDUCTIONS_BASIC, REDUCTIONS_GPURESIDENT, 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(), Vector::x, XXXBIGREAL, Vector::y, and Vector::z.

                                    :
         computeChecksum(0), marginViolations(0), pairlistWarnings(0),
         simParams(Node::Object()->simParameters),
         state(s),
         collection(CollectionMaster::Object()),
         startCTime(0),
         firstCTime(CmiTimer()),
         startWTime(0),
         firstWTime(namdWallTimer()),
         startBenchTime(0),
         stepInFullRun(0),
         ldbSteps(0),
         fflush_count(3)
 {
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
     if (simParams->CUDASOAintegrateMode) {
         CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
         PatchData* patchData = cpdata.ckLocalBranch();
         broadcast = new ControllerBroadcasts(0 /* ldObjPtr */, patchData->nodeBroadcast);
     } else 
 #endif  // defined(NAMD_CUDA) || defined(NAMD_HIP)
     {
         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
        reductionBasic = ReductionMgr::Object()->willRequire(REDUCTIONS_AMD);
     } else {
        amd_reduction = NULL;
        submit_reduction = NULL;
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
        if (simParams->CUDASOAintegrateMode) {
          reductionGpuResident = ReductionMgr::Object()->willRequire(REDUCTIONS_GPURESIDENT);
        }
 #endif  // defined(NAMD_CUDA) || defined(NAMD_HIP)
        reductionBasic = 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;
 
     // Monte Carlo pressure control
     if (simParams->monteCarloPressureOn) {
       monteCarloMaxVolume.x = simParams->monteCarloMaxVolume.x;
       monteCarloMaxVolume.y = simParams->monteCarloMaxVolume.y;
       monteCarloMaxVolume.z = simParams->monteCarloMaxVolume.z;
       mc_totalTry = mc_totalAccept = 0;
       for(int ax = 0; ax < MC_AXIS_TOTAL; ++ax) {
         mc_trial[ax] = 0;
         mc_accept[ax] = 0;
       }
     }
 
     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;
 
     bondedEnergy_ti_1 = 0;
     bondedEnergy_ti_2 = 0;
     ljEnergy_ti_1 = 0;
     ljEnergy_ti_2 = 0;
     electEnergy_ti_1 = 0;
     electEnergy_ti_2 = 0;
     electEnergySlow_ti_1 = 0;
     electEnergySlow_ti_2 = 0;
     TIderivative = 0;
 }

◆ ~Controller()

Controller::~Controller ( void )

virtual

Definition at line 331 of file Controller.C.

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

 {
     delete broadcast;
     if (reductionBasic) {
       delete reductionBasic;
     }
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
     if (reductionGpuResident) {
       delete reductionGpuResident;
     }
 #endif
     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

◆ adaptTempInit()

void Controller::adaptTempInit ( int step )

protected

Definition at line 2840 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");
     }
 }

◆ adaptTempUpdate()

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

protected

Definition at line 2966 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, minimize(), 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;
    }
    
 }

◆ adaptTempWriteRestart()

void Controller::adaptTempWriteRestart ( int step )

protected

Definition at line 2940 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(); 
     }
 }    

◆ algorithm()

void Controller::algorithm ( void )

protectedvirtual

Definition at line 379 of file Controller.C.

References BackEnd::awaken(), broadcast, checkpoint_lattice, checkpoint_state, checkpoint_stored, checkpoints, END_OF_RUN, endi(), enqueueCollections(), EVAL_MEASURE, FILE_OUTPUT, SimParameters::firstTimestep, FORCE_OUTPUT, SimpleBroadcastObject< T >::get(), SimParameters::initialTemp, integrate(), iout, Controller::checkpoint::lattice, 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) {
 #ifdef NODEGROUP_FORCE_REGISTER
     CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
     PatchData *patchData = cpdata.ckLocalBranch();
 #endif
     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();
 }

◆ awaken()

void Controller::awaken ( void )

Definition at line 371 of file Controller.C.

References Controller::CthThreadWrapper::thread.

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

                             {
   CthAwaken(threadWrapper->thread);
 }

◆ berendsenPressure()

void Controller::berendsenPressure ( int step )

protected

Definition at line 1370 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);
     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)
     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;
   }
 }

◆ calc_accelMDG_E_k()

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 2225 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);
 }

◆ calc_accelMDG_force_factor()

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 2254 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.;
    }
 }

◆ calc_accelMDG_mean_std()

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 2205 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);
 }

◆ calcPressure()

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

References Lattice::a(), SimParameters::accelMDOn, AVGXY, Lattice::b(), Lattice::c(), controlPressure, controlPressure_nbond, controlPressure_normal, controlPressure_slow, SimParameters::CUDASOAintegrate, Tensor::diagonal(), SimParameters::getCurrentLambda(), Molecule::getVirialTailCorr(), groupPressure, groupPressure_nbond, groupPressure_normal, groupPressure_slow, Tensor::identity(), SimParameters::isMultiTimeStepping(), SimParameters::LJcorrection, SimParameters::LJcorrectionAlt, minimize(), NAMD_bug(), nbondFreq, Node::Object(), Lattice::origin(), outer(), pressure, pressure_amd, pressure_nbond, pressure_normal, pressure_slow, simParams, slowFreq, state, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, SimParameters::useGroupPressure, virial_amd, Lattice::volume(), Vector::x, Tensor::xx, Tensor::xy, Tensor::xz, Vector::y, Tensor::yx, Tensor::yy, Tensor::yz, Vector::z, Tensor::zx, Tensor::zy, and Tensor::zz.

Referenced by multigratorPressure(), and receivePressure().

                                                                                             {
 
   Tensor virial_normal = virial_normal_in;
   Tensor virial_nbond = virial_nbond_in;
   Tensor virial_slow = virial_slow_in;
 
 #if 0
   Tensor tmp = virial_normal;
   fprintf(stderr, "VIRIAL_NORMAL  %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);
   
   tmp = virial_nbond;
   fprintf(stderr, "VIRIAL_NBOND  %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);
   
   tmp = virial_slow;
   fprintf(stderr, "VIRIAL_SLOW  %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);
      
   tmp = intVirial_normal;
   fprintf(stderr, "INT_VIRIAL_NORMAL  %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);
   
   tmp = intVirial_nbond;
   fprintf(stderr, "INT_VIRIAL_NBOND  %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);
   
   tmp = intVirial_slow;
   fprintf(stderr, "INT_VIRIAL_SLOW  %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);
   
   Vector ext = extForce_normal;
   fprintf(stderr, "EXT_FORCE_NORMAL  %1.2lf %1.2lf %1.2lf\n",
      ext.x, ext.y, ext.z);
   
   ext = extForce_nbond;
   fprintf(stderr, "EXT_FORCE_NBOND  %1.2lf %1.2lf %1.2lf\n",
      ext.x, ext.y, ext.z);
   
   ext = extForce_slow;
   fprintf(stderr, "EXT_FORCE_SLOW  %1.2lf %1.2lf %1.2lf\n",
      ext.x, ext.y, ext.z);
 
   //fprintf(stderr, "State Lattice for timestep %d\n", step);
   fprintf(stderr, "LATTICE  %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", 
     state->lattice.a().x, state->lattice.a().y, state->lattice.a().z,
     state->lattice.b().x, state->lattice.b().y, state->lattice.b().z,
     state->lattice.c().x, state->lattice.c().y, state->lattice.c().z);
     
   fprintf(stderr, "VOLUME  %3.10lf\n", state->lattice.volume());
   
 #endif
 
 
   Node *node = Node::Object();
   Molecule *molecule = node->molecule;
   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 ((simParams->LJcorrection || simParams->LJcorrectionAlt) && volume) {
 #ifdef MEM_OPT_VERSION
       NAMD_bug("LJcorrection not supported in memory optimized build.");
 #else
       // Apply tail correction to pressure
       BigReal alchLambda = simParams->getCurrentLambda(step);
       virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
 #endif
     }
 
 
     // kinetic energy component included in virials
     if(simParams->CUDASOAintegrate && step != 0 && !simParams->isMultiTimeStepping()){
       pressure_normal = virial_normal / volume;
       // groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
       if (simParams->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 = pressure - intVirial_normal / volume;
       // groupPressure = groupPressure_normal + groupPressure_nbond +
       //       groupPressure_slow;
     }else{
       pressure_normal = virial_normal / volume;
 
       groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
       if (simParams->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 ( simParams->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 ( simParams->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 ( 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(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.);
   }
 }

◆ compareChecksums()

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

protected

Definition at line 3264 of file Controller.C.

Referenced by getTotalPotentialEnergy(), printDynamicsEnergies(), and printMinimizeEnergies().

                                                          {
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     bool cudaIntegrator = simParams->CUDASOAintegrate;
     RequireReduction* reduction = getCurrentReduction();
 
     // Some basic correctness checking
     BigReal checksum, checksum_b;
     char errmsg[256];
 
     checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
     if ( ((int)checksum) != molecule->numAtoms && step != 0) {
       sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
               (int)checksum, molecule->numAtoms);
       NAMD_bug(errmsg);
     }
     // do we ever need this if we're on a GPU-based code?
     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_ONEFOURNBTHOLE_CHECKSUM);
     if ( checksum_b && (((int)checksum) != molecule->numCalcOneFourNbTholes) ) {
       sprintf(errmsg, "Bad global 1-4 NbThole count! (%d vs %d)\n",
               (int)checksum, molecule->numCalcOneFourNbTholes);
       if ( forgiving && (((int)checksum) < molecule->numCalcOneFourNbTholes) )
         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 << ") on step " << step;
         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");
     }
 }

◆ computeAlchWork()

BigReal Controller::computeAlchWork ( const int step )

protected

Definition at line 4677 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
   );
 }

◆ correctMomentum()

void Controller::correctMomentum ( int step )

protected

Definition at line 1673 of file Controller.C.

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

Referenced by integrate().

                                          {
     RequireReduction* reduction = getCurrentReduction();
 
     Vector momentum;
     BigReal mass;
     momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
     momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
     momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
     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);
     // What does this mean???
     broadcast->momentumCorrection.publish(step+slowFreq,momentum);
 }

◆ cycleBarrier()

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

protected

Definition at line 4967 of file Controller.C.

References broadcast, and namdWallTimer().

Referenced by integrate().

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

◆ enqueueCollections()

void Controller::enqueueCollections ( int timestep )

protected

Definition at line 4424 of file Controller.C.

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

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

 {
   int prec;
   int dcdIndex;
 
   std::tie(prec, dcdIndex) = Output::coordinateNeeded(timestep);
 #ifndef MEM_OPT_VERSION
   if ( prec == 4 ){
     collection->enqueuePositionsDcdSelection(timestep,state->lattice);
   }
   else
 #endif
     if ( prec ){
     collection->enqueuePositions(timestep,state->lattice);
   }
   if ( Output::velocityNeeded(timestep) )
     collection->enqueueVelocities(timestep);
   if ( Output::forceNeeded(timestep) )
     collection->enqueueForces(timestep);
 }

◆ getTIderivative()

BigReal Controller::getTIderivative ( void ) const

inline

Definition at line 121 of file Controller.h.

Referenced by colvarproxy_namd::get_dE_dlambda().

121 { return TIderivative; }

◆ getTotalPotentialEnergy()

BigReal Controller::getTotalPotentialEnergy ( int step )

protected

Definition at line 3572 of file Controller.C.

References compareChecksums(), SimParameters::CUDASOAintegrate, electEnergy, electEnergySlow, SimParameters::getCurrentLambda(), Molecule::getEnergyTailCorr(), RequireReduction::item(), SimParameters::LJcorrection, ljEnergy, ljEnergySlow, Node::molecule, NAMD_bug(), nbondFreq, Node::Object(), REDUCTION_ANGLE_ENERGY, REDUCTION_BC_ENERGY, REDUCTION_BOND_ENERGY, REDUCTION_CROSSTERM_ENERGY, REDUCTION_DIHEDRAL_ENERGY, REDUCTION_ELECT_ENERGY, REDUCTION_ELECT_ENERGY_SLOW, REDUCTION_IMPROPER_ENERGY, REDUCTION_LJ_ENERGY, REDUCTION_LJ_ENERGY_SLOW, REDUCTION_MISC_ENERGY, simParams, slowFreq, state, and Lattice::volume().

Referenced by monteCarloPressure_accept(), and monteCarloPressure_prepare().

 {
     compareChecksums(step, 0);
     
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     Lattice &lattice = state->lattice;
     BigReal volume = lattice.volume();
     // store total energy
     BigReal totalPotentialEnergy = 0.0;
     BigReal bondEnergy, angleEnergy, dihedralEnergy;
     BigReal improperEnergy, crosstermEnergy;
     BigReal electEnergy, ljEnergy, electEnergySlow, ljEnergySlow;
     BigReal miscEnergy, bcEnergy;
     bool cudaIntegrator = simParams->CUDASOAintegrate;
     RequireReduction* reduction = getCurrentReduction();
 
     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);
     bcEnergy = reduction->item(REDUCTION_BC_ENERGY);
     miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
 
     totalPotentialEnergy += bondEnergy + angleEnergy + dihedralEnergy +
                   improperEnergy + crosstermEnergy + miscEnergy
                   + bcEnergy;
 
     if (! ( step % nbondFreq ))
     {
       electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
       ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
       totalPotentialEnergy += electEnergy + ljEnergy;
     }
 
     if (! ( step % slowFreq ))
     {
       electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
       ljEnergySlow = reduction->item(REDUCTION_LJ_ENERGY_SLOW);
       totalPotentialEnergy += electEnergySlow;
       totalPotentialEnergy += ljEnergySlow;
     }
 
     if (simParams->LJcorrection && volume) {
 #ifdef MEM_OPT_VERSION
       NAMD_bug("LJcorrection not supported in memory optimized build.");
 #else
       // Apply tail correction to energy.
       BigReal alchLambda = simParams->getCurrentLambda(step);
       totalPotentialEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
 #endif
     }
 
   #if 0
     printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, Bond: %f, Angle: %f, Dih: %f, IMP: %f, Cross: %f\n", 
         step, CkMyPe(), bondEnergy, angleEnergy, dihedralEnergy, improperEnergy, crosstermEnergy);
     printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, LJ: %f, Real: %f, Recip: %f\n", step, CkMyPe(), 
         ljEnergy, electEnergy, electEnergySlow);
     printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, MISC: %f, BC: %f, TOTAL: %f\n", step, CkMyPe(), 
         miscEnergy, bcEnergy, totalPotentialEnergy);
   #endif
 
     return totalPotentialEnergy;
 }

◆ integrate()

void Controller::integrate ( int scriptTask )

protected

Definition at line 531 of file Controller.C.

References adaptTempUpdate(), berendsenPressure(), correctMomentum(), SimParameters::CUDASOAintegrate, cycleBarrier(), enqueueCollections(), eventEndOfTimeStep, 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(!simParams->CUDASOAintegrate){
      
      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
     }
   }else CthSuspend();
     // signal(SIGINT, oldhandler);
 }

◆ langevinPiston1()

void Controller::langevinPiston1 ( int step )

protected

Definition at line 1405 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;
     }
 
   }
 
 }

◆ langevinPiston2()

void Controller::langevinPiston2 ( int step )

protected

Definition at line 1560 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;
     }
   }
 
 
 }

◆ minimize()

void Controller::minimize ( )

protected

Definition at line 780 of file Controller.C.

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

Referenced by adaptTempUpdate(), algorithm(), calcPressure(), printEnergies(), receivePressure(), and rescaleaccelMD().

                           {
   // 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;
   }
 
 }

◆ monteCarloPressure_accept()

void Controller::monteCarloPressure_accept ( int step )

protected

Calculate the MC acceptance criteria for volume change and determin if this volume change is accepted or not.

Definition at line 1137 of file Controller.C.

References BOLTZMANN, broadcast, Lattice::c(), endi(), getTotalPotentialEnergy(), SimParameters::GPUresidentSingleProcessMode, iout, iWARN(), mc_accept, mc_oldLattice, mc_picked_axis, mc_totalAccept, mc_totalEnergyOld, mc_totalTry, mc_trial, MC_X, MC_Y, MC_Z, Node::molecule, SimParameters::monteCarloAcceptanceRate, SimParameters::monteCarloAdjustmentFreq, ControllerBroadcasts::monteCarloBarostatAcceptance, ControllerState::monteCarloMaxVolume, SimParameters::monteCarloPressureFreq, SimParameters::monteCarloPressureOn, SimParameters::monteCarloPressureTarget, SimParameters::monteCarloTemp, Molecule::numMolecules, Node::Object(), SimParameters::outputEnergies, SimpleBroadcastObject< T >::publish(), random, RequireReduction::require(), simParams, state, SimParameters::surfaceTensionTarget, Random::uniform(), SimParameters::useConstantArea, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Vector::x, Vector::y, and Vector::z.

 {
   if (simParams->monteCarloPressureOn) {
     if ( !(step % simParams->monteCarloPressureFreq) ) {    
       Node *node = Node::Object();
       Molecule *molecule = node->molecule;
       // need to consider the fix atom case!
       int numAtom = molecule->numMolecules;
       BigReal volumeOld = mc_oldLattice.volume();
       BigReal volumeNew = state->lattice.volume();
       BigReal areaNew = (volumeNew / state->lattice.c().z);
       BigReal areaOld = (volumeOld / mc_oldLattice.c().z);
       BigReal kbT = BOLTZMANN * simParams->monteCarloTemp;
       BigReal pressureTarget = simParams->monteCarloPressureTarget;
       BigReal surfTensionTarget = simParams->surfaceTensionTarget;
 
       // When using the GPU resident code path in single process mode, the reduction
       // data needs to persist so it can be consumed by the printStep function
       const bool clear_reduction_data = !simParams->GPUresidentSingleProcessMode;
       RequireReduction* reduction = getCurrentReduction();
       reduction->require(clear_reduction_data);
       
       BigReal totalEnergyOld = mc_totalEnergyOld;
       BigReal totalEnergyNew = getTotalPotentialEnergy(step);
       BigReal weight = (totalEnergyNew - totalEnergyOld);
       weight += (pressureTarget * (volumeNew - volumeOld));
       weight -= (numAtom * kbT * log(volumeNew / volumeOld));
       weight -= (surfTensionTarget * (areaNew - areaOld));
       BigReal totalWeight = exp(-weight / kbT);
       int accepted = (random->uniform() < totalWeight);
       //accepted = 0;
       #if 0
       printf("MC-data (accept): step: %d, Pe: %d, numAtom: %d\n", step, CkMyPe(), numAtom);
       printf("                  oldE: %f, newE: %f, deltaE: %f\n", totalEnergyOld, totalEnergyNew, totalEnergyNew - totalEnergyOld);
       printf("                  oldV: %f, newV: %f, deltaV: %f\n", volumeOld, volumeNew, volumeNew - volumeOld);
       printf("                  surfTen: %f, deltaA: %f \n", surfTensionTarget, areaNew - areaOld);
       printf("                  weight: %f, accepted: %d\n\n", totalWeight, accepted);
       #endif
       if(accepted) {
         // what should I do if its accepted?
         ++mc_totalAccept;
         ++mc_accept[mc_picked_axis];
       } else {
         // what should I do if its rejected?
         state->lattice = mc_oldLattice;
       }
       // update the maximum scale no matter if it's accepted or not
       {
         if (mc_trial[mc_picked_axis] == simParams->monteCarloAdjustmentFreq) {
           BigReal factor = 1.0;
           BigReal accRate = (BigReal) mc_accept[mc_picked_axis] / (BigReal) mc_trial[mc_picked_axis];
           // safe check to make sure there is any accepted move
           if (accRate > 0.0) {
             factor = accRate / simParams->monteCarloAcceptanceRate;
           } else {
             factor = 0.5; // if no move was accepted, we scale by half
           }
           // reset the trial for picked axis
           mc_trial[mc_picked_axis] = mc_accept[mc_picked_axis] = 0;
           BigReal maxVolume = state->lattice.volume() * 0.3;
           switch (mc_picked_axis) {
             case MC_X:
               monteCarloMaxVolume.x *= factor;
               if (monteCarloMaxVolume.x < 0.001) {
                 iout << iWARN << "Small volume change in MC barostat!\n" << endi;
                 iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
                 monteCarloMaxVolume.x = 0.001;
               } else if (monteCarloMaxVolume.x > maxVolume ) {
                 iout << iWARN << "Large volume change in MC barostat!\n" << endi;
                 iout << iWARN << "Maximum volume change is set to " <<  maxVolume << " A^3.\n" << endi;
                 monteCarloMaxVolume.x = maxVolume;
               }
               break;
 
             case MC_Y:
               monteCarloMaxVolume.y *= factor;
               if (monteCarloMaxVolume.y < 0.001) {
                 iout << iWARN << "Small volume change in MC barostat!\n" << endi;
                 iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
                 monteCarloMaxVolume.y = 0.001;
               } else if (monteCarloMaxVolume.y > maxVolume ) {
                 iout << iWARN << "Large volume change in MC barostat!\n" << endi;
                 iout << iWARN << "Maximum volume change is set to " <<  maxVolume << " A^3.\n" << endi;
                 monteCarloMaxVolume.y = maxVolume;
               }
               break;
 
             case MC_Z:
               monteCarloMaxVolume.z *= factor;
               if (monteCarloMaxVolume.z < 0.001) {
                 iout << iWARN << "Small volume change in MC barostat!\n" << endi;
                 iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
                 monteCarloMaxVolume.z = 0.001;
               } else if (monteCarloMaxVolume.z > maxVolume ) {
                 iout << iWARN << "Large volume change in MC barostat!\n" << endi;
                 iout << iWARN << "Maximum volume change is set to " <<  maxVolume << " A^3.\n" << endi;
                 monteCarloMaxVolume.z = maxVolume;
               }
               break;
           }
         }
       }
 
       if ( !(step % simParams->outputEnergies) ) {
         BigReal acceptPerc = 100.0 * ((BigReal) mc_totalAccept /
                             (BigReal) mc_totalTry);
         iout << "MC_BAROSTAT STEP: " << step << " TRIALS: " << mc_totalTry << 
           " ACCEPTED: " << mc_totalAccept << " %ACCEPTANCE: " << acceptPerc;
         if (simParams->useFlexibleCell) {
           if (simParams->useConstantArea) {
             iout << " MAX_VOLUME_Z: " << monteCarloMaxVolume.z << "\n\n";
           } else if (simParams->useConstantRatio) {
             iout << " MAX_VOLUME_XY: " << monteCarloMaxVolume.x << 
               " MAX_VOLUME_Z: " << monteCarloMaxVolume.z << "\n\n";
           } else {
             iout << " MAX_VOLUME_X: " << monteCarloMaxVolume.x
                  << " MAX_VOLUME_Y: " << monteCarloMaxVolume.y
                  << " MAX_VOLUME_Z: " << monteCarloMaxVolume.z
                  << "\n\n";
           }
         } else {
           iout << " MAX_VOLUME_XYZ: " << monteCarloMaxVolume.x << "\n\n";
         }
       }
 
       broadcast->monteCarloBarostatAcceptance.publish(step, accepted);
     }
   }
 }

◆ monteCarloPressure_prepare()

void Controller::monteCarloPressure_prepare ( int step )

protected

Perform a random walk in volume and calculate the rescale factor for lattice and atom coordinates.

Definition at line 1267 of file Controller.C.

References broadcast, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, getTotalPotentialEnergy(), SimParameters::GPUresidentSingleProcessMode, Tensor::identity(), mc_oldLattice, mc_picked_axis, mc_totalEnergyOld, mc_totalTry, mc_trial, MC_X, MC_Y, MC_Z, ControllerState::monteCarloMaxVolume, SimParameters::monteCarloPressureFreq, SimParameters::monteCarloPressureOn, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), random, RequireReduction::require(), Lattice::rescale(), simParams, state, Random::uniform(), SimParameters::useConstantArea, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz.

 {
   if (simParams->monteCarloPressureOn) {
     if ( !(step % simParams->monteCarloPressureFreq) ) {   
       // Wait here for Sequencer to finish energy/force computation
       // When using the GPU resident code path in single process mode, the reduction
       // data needs to persist so it can be consumed by the MCPressure accept function
       const bool clear_reduction_data = !simParams->GPUresidentSingleProcessMode;
       RequireReduction* reduction = getCurrentReduction();
       reduction->require(clear_reduction_data);
       mc_oldLattice = state->lattice;
       mc_totalEnergyOld = getTotalPotentialEnergy(step);
       Tensor factor = Tensor::identity(1.0);
       // pick a random deltaV and calc scaling factor
       {
         BigReal oldVolume = mc_oldLattice.volume();
         bool flexible = simParams->useFlexibleCell;
         bool constArea = simParams->useConstantArea;
         bool constRatio = simParams->useConstantRatio;
         bool constZ, constY, constX;
         if (simParams->fixCellDims) {
           constX = simParams->fixCellDimX;
           constY = simParams->fixCellDimY;
           constZ = simParams->fixCellDimZ;
         } else {
           constZ = constY = constX = false;
         }
 
         if(flexible) {
           //Anisotropic fluctuation
           if(constArea) {
             // always pick z
             mc_picked_axis = MC_Z;
             BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
             factor.zz = (oldVolume + dV) / oldVolume;
           } else if (constRatio) {
             while(true) {
               // decide to scale x-y or scale z
               BigReal probAxis = 2.0 * random->uniform();
               if (probAxis < 1.0) {
                 //picking x or y. We use MC_X for constant ratio
                 mc_picked_axis = MC_X;
                 BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
                 BigReal scale = sqrt((oldVolume + dV)/oldVolume);
                 factor.xx = factor.yy = scale;
                 break;
               } else if (probAxis < 2.0 && !constZ) {
                 //picking z
                 mc_picked_axis = MC_Z;
                 BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
                 BigReal scale = (oldVolume + dV)/oldVolume;
                 factor.zz = scale;
                 break;
               }
             }
           } else {
             while(true) {
               // decide to scale x, y or z
               BigReal probAxis = 3.0 * random->uniform();
               if (probAxis < 1.0 && !constX) {
                 // picking x
                 mc_picked_axis = MC_X;
                 BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
                 factor.xx = (oldVolume + dV) / oldVolume;
                 break;
               } else if (probAxis < 2.0 && !constY) {
                 // picking y
                 mc_picked_axis = MC_Y;
                 BigReal dV = monteCarloMaxVolume.y * (2.0 * random->uniform() - 1.0);
                 factor.yy = (oldVolume + dV) / oldVolume;
                 break;
               } else if (probAxis < 3.0 && !constZ) {
                 // picking z
                 mc_picked_axis = MC_Z;
                 BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
                 factor.zz = (oldVolume + dV) / oldVolume;
                 break;
               }
             }
           }
         } else {
           //Isotropic fluctuation.  We use MC_X for isotropic fluctuation
           mc_picked_axis = MC_X;
           BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
           BigReal scale = pow((oldVolume + dV)/oldVolume, 1.0/3.0);
           factor.xx = factor.yy = factor.zz = scale;
         } 
         // increament the MC trial
         ++mc_totalTry;
         ++mc_trial[mc_picked_axis];
       }
 
       broadcast->positionRescaleFactor.publish(step,factor);
       state->lattice.rescale(factor);
     }
   }
 }

◆ multigatorCalcEnthalpy()

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

protected

Definition at line 1106 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;
 }

◆ multigratorPressure()

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

protected

Definition at line 964 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_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
       RequireReduction* reduction = getCurrentReduction();
       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 );
     }
 
   }
 }

◆ multigratorTemperature()

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

protected

Definition at line 1040 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];
       }
     }
 
   }
 }

◆ outputExtendedSystem()

void Controller::outputExtendedSystem ( int step )

protected

Definition at line 4795 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(), NAMD_FILENAME_BUFFER_SIZE, 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[NAMD_FILENAME_BUFFER_SIZE];
       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[NAMD_FILENAME_BUFFER_SIZE];
     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;
     }
   }
 
 }

◆ outputFepEnergy()

void Controller::outputFepEnergy ( int step )

protected

Definition at line 4446 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 ) 
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
         &&
         // some people rely on the ensemble averages (dE_avg and dG) in fep output for calculating the free energy changes.
         // in the GPU version, dE is only available every outputEnergies steps.
         // so I check if the step is multiples of outputEnergies here.
         // I am still not sure whether we should check the stepInRun or step.
         // JH: dE is also calculated at alchOutFreq
         (step % simParams->computeEnergies == 0 || step % simParams->alchOutFreq == 0)
 #endif
     ) {
       // 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();
       }
     }
   }
 }

◆ outputTiEnergy()

void Controller::outputTiEnergy ( int step )

protected

Definition at line 4530 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;
   }
   // The energies are only computed if they are required to output in CUDA build.
   // In the case when step % simParams->computeEnergies != 0 the energies are 0.
   // All output terms in writeTiEnergyData use TiNo and recent_TiNo as denominators,
   // so the TiNo and recent_TiNo should NOT be increased.
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
   if (step % simParams->computeEnergies == 0) {
 #endif
     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;
     // Compute global dE / dLambda for lambda-dynamics
     computeTIderivative(step);
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
   }
 #endif
 
   // 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;
     }
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
     if (step % simParams->computeEnergies == 0) {
 #endif
       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 defined(NAMD_CUDA) || defined(NAMD_HIP)
     }
 #endif
 
     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();
     }
   }
  }
 }

◆ printDynamicsEnergies()

void Controller::printDynamicsEnergies ( int step )

protected

Definition at line 3567 of file Controller.C.

References compareChecksums(), and printEnergies().

Referenced by integrate().

                                                {
     compareChecksums(step);
     printEnergies(step,0);
 }

◆ printEnergies()

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

protected

Definition at line 3639 of file Controller.C.

Referenced by printDynamicsEnergies(), and printMinimizeEnergies().

 {
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     Lattice &lattice = state->lattice;
     bool cudaIntegrator = simParams->CUDASOAintegrate;
     RequireReduction* reduction = getCurrentReduction();
 
     // 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);
       ljEnergySlow = reduction->item(REDUCTION_LJ_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 ((simParams->LJcorrection || simParams->LJcorrectionAlt) && volume) {
 #ifdef MEM_OPT_VERSION
       NAMD_bug("LJcorrection not supported in memory optimized build.");
 #else
       // Apply tail correction to energy.
       BigReal alchLambda = simParams->getCurrentLambda(step);
       ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
 
       if (simParams->alchOn) {
         if (simParams->alchFepOn) {
           BigReal alchLambda2 = simParams->getCurrentLambda2(step);
           ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
         } else if (simParams->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 (simParams->alchLambdaFreq > 0) {
       if (step <= 
           simParams->firstTimestep + simParams->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 + ljEnergySlow
   + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy 
         + groLJEnergy + groGaussEnergy);
     // End of port
     totalEnergy = potentialEnergy + kineticEnergy;
     if ( ! cudaIntegrator) {
       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 <= simParams->firstTimestep &&
         (simParams->stochRescaleOn && simParams->stochRescaleHeat)) {
       heat = 0.0;
       totalEnergy0 = totalEnergy;
     }
     if ( simParams->outputMomenta && ! minimize &&
          ! ( step % simParams->outputMomenta ) )
     {
       iout << "MOMENTUM: " << step 
            << " P: " << momentum
            << " L: " << angularMomentum
            << "\n" << endi;
     }
 
     if ( simParams->outputPressure ) {
       if ( ! cudaIntegrator ) {
         pressure_tavg += pressure;
         groupPressure_tavg += groupPressure;
       }
       tavg_count += 1;
       if ( minimize || ! ( step % simParams->outputPressure ) ) {
         if (cudaIntegrator) {
           // tensors are valid for outputPressure step
           // update averages, then copy into tensor for output
           pressureAverage_xx.addSample(pressure.xx);
           pressureAverage_yx.addSample(pressure.yx);
           pressureAverage_yy.addSample(pressure.yy);
           pressureAverage_zx.addSample(pressure.zx);
           pressureAverage_zy.addSample(pressure.zy);
           pressureAverage_zz.addSample(pressure.zz);
           groupPressureAverage_xx.addSample(groupPressure.xx);
           groupPressureAverage_xy.addSample(groupPressure.xy);
           groupPressureAverage_xz.addSample(groupPressure.xz);
           groupPressureAverage_yx.addSample(groupPressure.yx);
           groupPressureAverage_yy.addSample(groupPressure.yy);
           groupPressureAverage_yz.addSample(groupPressure.yz);
           groupPressureAverage_zx.addSample(groupPressure.zx);
           groupPressureAverage_zy.addSample(groupPressure.zy);
           groupPressureAverage_zz.addSample(groupPressure.zz);
           pressure_tavg.xx = pressureAverage_xx.average();
           pressure_tavg.xy = pressureAverage_yx.average();
           pressure_tavg.xz = pressureAverage_zx.average();
           pressure_tavg.yx = pressureAverage_yx.average();
           pressure_tavg.yy = pressureAverage_yy.average();
           pressure_tavg.yz = pressureAverage_zy.average();
           pressure_tavg.zx = pressureAverage_zx.average();
           pressure_tavg.zy = pressureAverage_zy.average();
           pressure_tavg.zz = pressureAverage_zz.average();
           groupPressure_tavg.xx = groupPressureAverage_xx.average();
           groupPressure_tavg.xy = groupPressureAverage_xy.average();
           groupPressure_tavg.xz = groupPressureAverage_xz.average();
           groupPressure_tavg.yx = groupPressureAverage_yx.average();
           groupPressure_tavg.yy = groupPressureAverage_yy.average();
           groupPressure_tavg.yz = groupPressureAverage_yz.average();
           groupPressure_tavg.zx = groupPressureAverage_zx.average();
           groupPressure_tavg.zy = groupPressureAverage_zy.average();
           groupPressure_tavg.zz = groupPressureAverage_zz.average();
           iout << "PRESSURE: " << step << " "
              << PRESSUREFACTOR * pressure << "\n"
              << "GPRESSURE: " << step << " "
              << PRESSUREFACTOR * groupPressure << "\n";
           // tavg_count makes output behaves the same as for standard case
           // but is no longer part of the average
           if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
              << (PRESSUREFACTOR) * pressure_tavg << "\n"
              << "GPRESSAVG: " << step << " "
              << (PRESSUREFACTOR) * groupPressure_tavg << "\n";
           iout << endi;
           tavg_count = 0;
         }
         else {
           iout << "PRESSURE: " << step << " "
              << PRESSUREFACTOR * pressure << "\n"
              << "GPRESSURE: " << step << " "
              << PRESSUREFACTOR * groupPressure << "\n";
           if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
              << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
              << "GPRESSAVG: " << step << " "
              << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
           iout << endi;
           pressure_tavg = 0;
           groupPressure_tavg = 0;
           tavg_count = 0;
         }
       }
     }
 
     // pressure profile reductions
     if (pressureProfileSlabs) {
       const int freq = simParams->pressureProfileFreq;
       const int arraysize = 3*pressureProfileSlabs;
       
       BigReal *total = new BigReal[arraysize];
       memset(total, 0, arraysize*sizeof(BigReal));
       const int first = simParams->firstTimestep;
 
       if (ppbonded)    ppbonded->getData(first, step, lattice, total);
       if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
       if (ppint)       ppint->getData(first, step, lattice, total);
       for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
       pressureProfileCount++;
 
       if (!(step % freq)) {
         // convert NAMD internal virial to pressure in units of bar 
         BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
    
         iout << "PRESSUREPROFILE: " << step << " ";
         if (step == first) {
           // output pressure profile for this step
           for (int i=0; i<arraysize; i++) {
             iout << total[i] * scalefac << " ";
           }
         } else {
           // output pressure profile averaged over the last count steps.
           scalefac /= pressureProfileCount;
           for (int i=0; i<arraysize; i++) 
             iout << pressureProfileAverage[i]*scalefac << " ";
         }
         iout << "\n" << endi; 
        
         // Clear the average for the next block
         memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
         pressureProfileCount = 0;
       }
       delete [] total;
     }
   
    if ( step != simParams->firstTimestep || stepInFullRun == 0 ) {  // skip repeated first step
     if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
      int benchPhase = stepInFullRun / simParams->firstLdbStep;
      if ( benchPhase > 0 && benchPhase < 10 ) {
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
       if ( simParams->computeEnergies < 60 ) {
         iout << iWARN << "Energy evaluation is expensive, increase computeEnergies to improve performance.\n";
         iout << iWARN << "If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n";
       }
 #endif
       iout << iINFO;
       if ( benchPhase < 4 ) iout << "Initial time: ";
       else iout << "Benchmark time: ";
       iout << CkNumPes() << " CPUs ";
 
       {
         BigReal wallPerStep =
           (namdWallTimer() - startBenchTime) / simParams->firstLdbStep;
         BigReal ns = simParams->dt / 1000000.0;
         BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
 
         if (simParams->nsPerDayOn) {
           BigReal nsPerDay = ns / (wallPerStep * days);
           iout << wallPerStep << " s/step " << nsPerDay << " ns/day ";
         }
         else {
           BigReal daysPerNano = wallPerStep * days / ns;
           iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
         }
         iout << memusage_MB() << " MB memory\n" << endi;
       }
 
      }
      startBenchTime = namdWallTimer();
     }
     ++stepInFullRun;
    }
 
     printTiming(step);
 
     Vector pairVDWForce, pairElectForce;
     if ( simParams->pairInteractionOn ){
       GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
       GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
     }
 
     // Compute cumulative nonequilibrium work (including shadow work).
     if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
       work = totalEnergy - totalEnergy0 - heat;
     }
 
     // callback to Tcl with whatever we can
 #ifdef NAMD_TCL
 #define CALLBACKDATA(LABEL,VALUE) \
                 labels << (LABEL) << " "; values << (VALUE) << " ";
 #define CALLBACKLIST(LABEL,VALUE) \
                 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
     if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
 
       std::ostringstream labels, values;
 #if CMK_BLUEGENEL
       // the normal version below gives a compiler error
       values.precision(16);
 #else
       values << std::setprecision(16);
 #endif
       CALLBACKDATA("TS",step);
       CALLBACKDATA("BOND",bondEnergy);
       CALLBACKDATA("ANGLE",angleEnergy);
       CALLBACKDATA("DIHED",dihedralEnergy);
       CALLBACKDATA("CROSS",crosstermEnergy);
       CALLBACKDATA("IMPRP",improperEnergy);
       CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
       CALLBACKDATA("VDW",ljEnergy+ljEnergySlow);
       CALLBACKDATA("BOUNDARY",boundaryEnergy);
       CALLBACKDATA("MISC",miscEnergy);
       CALLBACKDATA("POTENTIAL",potentialEnergy);
       CALLBACKDATA("KINETIC",kineticEnergy);
       CALLBACKDATA("TOTAL",totalEnergy);
       CALLBACKDATA("TEMP",temperature);
       CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
       CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
       CALLBACKDATA("VOLUME",lattice.volume());
       CALLBACKLIST("CELL_A",lattice.a());
       CALLBACKLIST("CELL_B",lattice.b());
       CALLBACKLIST("CELL_C",lattice.c());
       CALLBACKLIST("CELL_O",lattice.origin());
       labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
                 << lattice.b_p() << " " << lattice.c_p() << "} ";
       if (simParams->drudeOn) {
         CALLBACKDATA("DRUDEBOND",drudeBondTemp);
       }
       if ( simParams->pairInteractionOn ) {
         CALLBACKLIST("VDW_FORCE",pairVDWForce);
         CALLBACKLIST("ELECT_FORCE",pairElectForce);
       }
       if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
         CALLBACKLIST("HEAT",heat);
         CALLBACKLIST("WORK",work);
       }
       if (simParams->alchOn) {
         if (simParams->alchThermIntOn) {
           CALLBACKLIST("BOND1", bondedEnergy_ti_1);
           CALLBACKLIST("ELEC1", electEnergy_ti_1 + electEnergySlow_ti_1 +
                                 electEnergyPME_ti_1);
           CALLBACKLIST("VDW1", ljEnergy_ti_1);
           CALLBACKLIST("BOND2", bondedEnergy_ti_2);
           CALLBACKLIST("ELEC2", electEnergy_ti_2 + electEnergySlow_ti_2 +
                                 electEnergyPME_ti_2);
           CALLBACKLIST("VDW2", ljEnergy_ti_2);
           if (simParams->alchLambdaFreq > 0) {
             CALLBACKLIST("CUMALCHWORK", cumAlchWork);
           }
         } else if (simParams->alchFepOn) {
           CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
                                 improperEnergy + bondedEnergyDiff_f);
           CALLBACKLIST("ELEC2", electEnergy_f + electEnergySlow_f);
           CALLBACKLIST("VDW2", ljEnergy_f);
         } 
       }
 
       labels << '\0';  values << '\0';  // insane but makes Linux work
       state->callback_labelstring = labels.str();
       state->callback_valuestring = values.str();
       // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
     }
 #undef CALLBACKDATA
 #endif
 
     if ( ! cudaIntegrator) {
       drudeBondTempAvg += drudeBondTemp;
       temp_avg += temperature;
       pressure_avg += trace(pressure)/3.;
       groupPressure_avg += trace(groupPressure)/3.;
       avg_count += 1;
     }
 
     if ( simParams->outputPairlists && pairlistWarnings &&
                                 ! (step % simParams->outputPairlists) ) {
       iout << iINFO << pairlistWarnings <<
         " pairlist warnings in past " << simParams->outputPairlists <<
         " steps.\n" << endi;
       pairlistWarnings = 0;
     }
 
     BigReal enthalpy;
     if (simParams->multigratorOn && ((step % simParams->computeEnergies) == 0)) {
       enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
     }
 
     if (simParams->IMDon) {
       IMDOutput *imd = NULL;
 #ifdef NODEGROUP_FORCE_REGISTER
       if (simParams->CUDASOAintegrate) {
         CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
         imd = cpdata.ckLocalBranch()->imd;
       } else
 #endif
       {
         imd = Node::Object()->imd;
       }
       if ( !(step % simParams->IMDfreq) &&
           ((simParams->IMDversion == IMDversion_t::IMDv3)
           && (simParams->IMDsendsettings.time_switch == 1)
           && (step != simParams->firstTimestep)) ) {
 
         IMDTime time = {simParams->dt * 0.001, simParams->dt * 0.001 * step, step}; // convert to ps
 
         if (imd != NULL) imd->gather_time(&time);
       }
 
       if (!(step % simParams->IMDfreq) && ((simParams->IMDversion == IMDversion_t::IMDv2) ||
           ((simParams->IMDversion == IMDversion_t::IMDv3) && (simParams->IMDsendsettings.energies_switch == 1) && (step != simParams->firstTimestep)))) {
         IMDEnergies energies;
         energies.tstep = step;
         energies.T = temp_avg/avg_count;
         energies.Etot = totalEnergy;
         energies.Epot = potentialEnergy;
         energies.Evdw = ljEnergy + ljEnergySlow;
         energies.Eelec = electEnergy + electEnergySlow;
         energies.Ebond = bondEnergy;
         energies.Eangle = angleEnergy;
         energies.Edihe = dihedralEnergy + crosstermEnergy;
         energies.Eimpr = improperEnergy;
         if (imd != NULL) imd->gather_energies(&energies);
       }
 
       if (!simParams->CUDASOAintegrate &&
           (simParams->IMDversion == IMDversion_t::IMDv3) &&
           !(step % simParams->IMDfreq) &&
           ((simParams->IMDsendsettings.time_switch == 1) || (simParams->IMDsendsettings.energies_switch == 1)) &&
           (step != simParams->firstTimestep)) {
         broadcast->IMDTimeEnergyBarrier.publish(step, 1);
       }
     }
     // XXX
     // Important note: in CPU-only/GPU-offload modes
     // printEnergies is called EVERY STEP.
     // Reduction averages are summed above and then later are output 
     // divided by count. Sums and counter are reset after printing.
     // XXX
     // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
     if ( ! minimize &&  step % simParams->outputEnergies ) return;
 
     if (cudaIntegrator) {
       totalEnergyAverage.addSample(totalEnergy);
       temperatureAverage.addSample(temperature);
       pressureAverage.addSample((1./3)*trace(pressure));
       groupPressureAverage.addSample((1./3)*trace(groupPressure));
     }
 
     // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
 
     if ( marginViolations ) {
       iout << iERROR << marginViolations <<
         " margin violations detected since previous energy output.\n" << endi;
     }
     marginViolations = 0;
 
     int precision = simParams->outputEnergiesPrecision;
     if ( (step % (10 * (minimize?1:simParams->outputEnergies) ) ) == 0 )
     {
         iout << "ETITLE:      TS";
         iout << FORMAT("BOND", precision);
         iout << FORMAT("ANGLE", precision);
         iout << FORMAT("DIHED", precision);
         if ( ! simParams->mergeCrossterms ) iout << FORMAT("CROSS", precision);
         iout << FORMAT("IMPRP", precision);
         iout << "     ";
         iout << FORMAT("ELECT", precision);
         iout << FORMAT("VDW", precision);
         iout << FORMAT("BOUNDARY", precision);
         iout << FORMAT("MISC", precision);
         iout << FORMAT("KINETIC", precision);
         iout << "     ";
         iout << FORMAT("TOTAL", precision);
         iout << FORMAT("TEMP", precision);
         iout << FORMAT("POTENTIAL", precision);
         if (cudaIntegrator) {
           iout << FORMAT("TOTALAVG", precision);
         }
         else {
           // iout << FORMAT("TOTAL2", precision);
           iout << FORMAT("TOTAL3", precision);
         }
         iout << FORMAT("TEMPAVG", precision);
         if ( volume != 0. ) {
           iout << "     ";
           iout << FORMAT("PRESSURE", precision);
           iout << FORMAT("GPRESSURE", precision);
           iout << FORMAT("VOLUME", precision);
           iout << FORMAT("PRESSAVG", precision);
           iout << FORMAT("GPRESSAVG", precision);
         }
         if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
           iout << "     ";
           iout << FORMAT("HEAT", precision);
           iout << FORMAT("WORK", precision);
         }
         if (simParams->drudeOn) {
           iout << "     ";
           iout << FORMAT("DRUDEBOND", precision);
           iout << FORMAT("DRBONDAVG", precision);
         }
         // Ported by JLai
         if (simParams->goGroPair) {
           iout << "     ";
           iout << FORMAT("GRO_PAIR_LJ", precision);
           iout << FORMAT("GRO_PAIR_GAUSS", precision);
         }
 
         if (simParams->goForcesOn) {
           iout << "     ";
           iout << FORMAT("NATIVE", precision);
           iout << FORMAT("NONNATIVE", precision);
           //iout << FORMAT("REL_NATIVE", precision);
           //iout << FORMAT("REL_NONNATIVE", precision);
           iout << FORMAT("GOTOTAL", precision);
           //iout << FORMAT("GOAVG", precision);
         }
         // End of port -- JLai
 
         if (simParams->alchOn) {
           if (simParams->alchThermIntOn) {
             iout << "\nTITITLE:     TS";
             iout << FORMAT("BOND1", precision);
             iout << FORMAT("ELECT1", precision);
             iout << FORMAT("VDW1", precision);
             iout << FORMAT("BOND2", precision);
             iout << "     ";
             iout << FORMAT("ELECT2", precision);
             iout << FORMAT("VDW2", precision);
             if (simParams->alchLambdaFreq > 0) {
               iout << FORMAT("LAMBDA", precision);
               iout << FORMAT("ALCHWORK", precision);
               iout << FORMAT("CUMALCHWORK", precision);
             }
           } else if (simParams->alchFepOn) {
             iout << "\nFEPTITLE:    TS";
             iout << FORMAT("BOND2", precision);
             iout << FORMAT("ELECT2", precision);
             iout << FORMAT("VDW2", precision);
             if (simParams->alchLambdaFreq > 0) {
               iout << FORMAT("LAMBDA", precision);
             }
           }
         }
 
         iout << "\n\n" << endi;
         
         if (simParams->qmForcesOn) {
             iout << "QMETITLE:      TS";
             iout << FORMAT("QMID", precision);
             iout << FORMAT("ENERGY", precision);
             if (simParams->PMEOn) iout << FORMAT("PMECORRENERGY", precision);
             iout << "\n\n" << endi;
         }
         
     }
 
     // N.B.  HP's aCC compiler merges FORMAT calls in the same expression.
     //       Need separate statements because data returned in static array.
     iout << ETITLE(step);
     iout << FORMAT(bondEnergy, precision);
     iout << FORMAT(angleEnergy, precision);
     if ( simParams->mergeCrossterms ) {
       iout << FORMAT(dihedralEnergy+crosstermEnergy, precision);
     } else {
       iout << FORMAT(dihedralEnergy, precision);
       iout << FORMAT(crosstermEnergy, precision);
     }
     iout << FORMAT(improperEnergy, precision);
     iout << "     ";
     iout << FORMAT(electEnergy+electEnergySlow, precision);
     iout << FORMAT(ljEnergy+ljEnergySlow, precision);
     iout << FORMAT(boundaryEnergy, precision);
     iout << FORMAT(miscEnergy, precision);
     iout << FORMAT(kineticEnergy, precision);
     iout << "     ";
     iout << FORMAT(totalEnergy, precision);
     iout << FORMAT(temperature, precision);
     iout << FORMAT(potentialEnergy, precision);
     if (cudaIntegrator) {
       iout << FORMAT(totalEnergyAverage.average(), precision);
       iout << FORMAT(temperatureAverage.average(), precision);
     }
     else {
       // iout << FORMAT(flatEnergy, precision);
       iout << FORMAT(smoothEnergy, precision);
       iout << FORMAT(temp_avg/avg_count, precision);
     }
     if ( volume != 0. )
     {
         iout << "     ";
         iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3., precision);
         iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3., precision);
         iout << FORMAT(volume, precision);
         if (cudaIntegrator) {
           iout << FORMAT(pressureAverage.average()*PRESSUREFACTOR, precision);
           iout << FORMAT(groupPressureAverage.average()*PRESSUREFACTOR, precision);
         }
         else {
           iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count, precision);
           iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count, precision);
         }
     }
     if (simParams->stochRescaleOn && simParams->stochRescaleHeat) {
           iout << "     ";
           iout << FORMAT(heat, precision);
           iout << FORMAT(work, precision);
     }
     if (simParams->drudeOn) {
         iout << "     ";
         iout << FORMAT(drudeBondTemp, precision);
         iout << FORMAT(drudeBondTempAvg/avg_count, precision);
     }
     // Ported by JLai
     if (simParams->goGroPair) {
       iout << "     ";
       iout << FORMAT(groLJEnergy, precision);
       iout << FORMAT(groGaussEnergy, precision);
     }
 
     if (simParams->goForcesOn) {
       iout << "     ";
       iout << FORMAT(goNativeEnergy, precision);
       iout << FORMAT(goNonnativeEnergy, precision);
       //iout << FORMAT(relgoNativeEnergy, precision);
       //iout << FORMAT(relgoNonnativeEnergy, precision);
       iout << FORMAT(goTotalEnergy, precision);
       //iout << FORMAT("not implemented", precision);
     } // End of port -- JLai
 
     if (simParams->alchOn) { 
       if (simParams->alchThermIntOn) {
         iout << "\n";
         iout << TITITLE(step);
         iout << FORMAT(bondedEnergy_ti_1, precision);
         iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 + 
                        electEnergyPME_ti_1, precision);
         iout << FORMAT(ljEnergy_ti_1, precision);
         iout << FORMAT(bondedEnergy_ti_2, precision);
         iout << "     ";
         iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 +
                        electEnergyPME_ti_2, precision);
         iout << FORMAT(ljEnergy_ti_2, precision);
         if (simParams->alchLambdaFreq > 0) {
           iout << FORMAT(simParams->getCurrentLambda(step), precision);
           iout << FORMAT(alchWork, precision);
           iout << FORMAT(cumAlchWork, precision);
         }
       } else if (simParams->alchFepOn) {
         iout << "\n";
         iout << FEPTITLE2(step);
         iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy 
                        + improperEnergy + bondedEnergyDiff_f, precision);
         iout << FORMAT(electEnergy_f + electEnergySlow_f, precision);
         iout << FORMAT(ljEnergy_f, precision);
         if (simParams->alchLambdaFreq > 0) {
           iout << FORMAT(simParams->getCurrentLambda(step), precision);
         }
       }
     }
 
     iout << "\n\n" << endi;
 
 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
      char webout[80];
      sprintf(webout,"%d %d %d %d",(int)totalEnergy,
              (int)(potentialEnergy),
              (int)kineticEnergy,(int)temperature);
      CApplicationDepositNode0Data(webout);
 #endif
 
     if (simParams->pairInteractionOn) {
       iout << "PAIR INTERACTION:";
       iout << " STEP: " << step;
       iout << " VDW_FORCE: ";
       iout << FORMAT(pairVDWForce.x, precision);
       iout << FORMAT(pairVDWForce.y, precision);
       iout << FORMAT(pairVDWForce.z, precision);
       iout << " ELECT_FORCE: ";
       iout << FORMAT(pairElectForce.x, precision);
       iout << FORMAT(pairElectForce.y, precision);
       iout << FORMAT(pairElectForce.z, precision);
       iout << "\n" << endi;
     }
     drudeBondTempAvg = 0;
     temp_avg = 0;
     pressure_avg = 0;
     groupPressure_avg = 0;
     avg_count = 0;
 
     if ( fflush_count ) {
       --fflush_count;
       fflush(stdout);
     }
 }

◆ printFepMessage()

void Controller::printFepMessage ( int step )

protected

Definition at line 1694 of file Controller.C.

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

Referenced by integrate().

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

◆ printMinimizeEnergies()

void Controller::printMinimizeEnergies ( int step )

protected

Definition at line 3537 of file Controller.C.

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

                                                {
 
     rescaleaccelMD(step,1);
     receivePressure(step,1);
     RequireReduction* reduction = getCurrentReduction();
 
     Node *node = Node::Object();
     Molecule *molecule = node->molecule;
     compareChecksums(step,1);
 
     printEnergies(step,1);
 
     min_energy = totalEnergy;
     min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
     min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
     min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
     min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
 
 #ifdef DEBUG_MINIMIZE
     printf("%s, line %d\n", __FILE__, __LINE__);
     printf("Step %d:\n", step);
     printf("   min_energy = %f\n", min_energy);
     printf("   min_f_dot_f = %f\n", min_f_dot_f);
     printf("   min_f_dot_v = %f\n", min_f_dot_v);
     printf("   min_v_dot_v = %f\n", min_v_dot_v);
     printf("   min_huge_count = %d\n", min_huge_count);
     printf("\n");
 #endif
 }

◆ printTiMessage()

void Controller::printTiMessage ( int step )

protected

Definition at line 1714 of file Controller.C.

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

Referenced by integrate().

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

◆ printTiming()

void Controller::printTiming ( int step )

protected

Definition at line 3467 of file Controller.C.

References RunningAverage::add_sample(), RunningAverage::average(), SimParameters::benchmarkTime, SimParameters::computeEnergies, SimParameters::dt, endi(), fflush_count, SimParameters::firstTimestep, iout, iWARN(), memusage_MB(), SimParameters::N, NAMD_quit(), namdWallTimer(), SimParameters::nsPerDayOn, SimParameters::outputPerformance, SimParameters::outputTiming, perfstats, simParams, and RunningAverage::standard_deviation().

Referenced by printEnergies().

                                      {
 
     if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
     {
       const double endWTime = namdWallTimer() - firstWTime;
       const double endCTime = CmiTimer() - firstCTime;
 
       // fflush about once per minute
       if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
 
       const double elapsedW = 
         (endWTime - startWTime) / simParams->outputTiming;
       const double elapsedC = 
         (endCTime - startCTime) / simParams->outputTiming;
 
       const double remainingW = elapsedW * (simParams->N - step);
       const double remainingW_hours = remainingW / 3600;
 
       startWTime = endWTime;
       startCTime = endCTime;
 
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
       if ( simParams->computeEnergies < 60 &&
            step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
         iout << iWARN << "Energy evaluation is expensive, increase computeEnergies to improve performance.\n" << endi;
         iout << iWARN << "If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n" << endi;
       }
 #endif
       if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
         if (simParams->nsPerDayOn) {
           BigReal ns = simParams->dt / 1000000.0;
           BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
           BigReal nsPerDay = ns / (elapsedW * days);
           CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
                     ", %g ns/days"
                     ", %g hours remaining, %f MB of memory in use.\n",
                     step, endCTime, elapsedC, endWTime, elapsedW,
                     nsPerDay,
                     remainingW_hours, memusage_MB());
         }
         else {
           CmiPrintf("TIMING: %d  CPU: %g, %g/step  Wall: %g, %g/step"
                     ", %g hours remaining, %f MB of memory in use.\n",
                     step, endCTime, elapsedC, endWTime, elapsedW,
                     remainingW_hours, memusage_MB());
         }
         if (simParams->outputPerformance) {
           // calculate running average and sample variance of
           // the measured time per step (elapsedW)
           perfstats.add_sample(elapsedW);
           double t_avg = perfstats.average();
           double t_std = perfstats.standard_deviation();
           // convert t_avg to average number of nanoseconds simulated per day
           double ns_per_day = (simParams->dt / t_avg) * (60 * 60 * 24 * 1e-6);
           CmiPrintf("PERFORMANCE: %d  averaging %g ns/day, %g sec/step with"
               " standard deviation %g\n", step, ns_per_day, t_avg, t_std);
         }
         if ( fflush_count ) { --fflush_count; fflush(stdout); }
         double benchmarkTime = (double) simParams->benchmarkTime;
         if (benchmarkTime > 0 && benchmarkTime <= endWTime) {
           // terminate NAMD benchmark
           char s[100];
           sprintf(s, "Exceeded benchmark time limit %g seconds",
               benchmarkTime);
           NAMD_quit(s);
         }
       }
     }
 }

◆ reassignVelocities()

void Controller::reassignVelocities ( int step )

protected

Definition at line 1728 of file Controller.C.

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

Referenced by integrate().

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

◆ rebalanceLoad()

void Controller::rebalanceLoad ( int step )

protected

Definition at line 4953 of file Controller.C.

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

Referenced by integrate().

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

◆ receivePressure()

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

protected

Definition at line 1840 of file Controller.C.

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

Referenced by integrate(), and printMinimizeEnergies().

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

◆ recvCheckpointAck()

void Controller::recvCheckpointAck ( checkpoint & cp )

protected

Definition at line 4944 of file Controller.C.

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

Referenced by Node::recvCheckpointAck().

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

◆ recvCheckpointReq()

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

protected

Definition at line 4914 of file Controller.C.

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

Referenced by Node::recvCheckpointReq().

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

◆ rescaleaccelMD()

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

protected

Definition at line 2328 of file Controller.C.

Referenced by integrate(), and printMinimizeEnergies().

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

◆ rescaleVelocities()

void Controller::rescaleVelocities ( int step )

protected

Definition at line 1651 of file Controller.C.

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

Referenced by integrate().

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

◆ resetMovingAverage()

void Controller::resetMovingAverage ( )

Definition at line 656 of file Controller.C.

References groupPressureAverage, SimParameters::movingAverageWindowSize, pressureAverage, MovingAverage::reset(), simParams, temperatureAverage, and totalEnergyAverage.

                                     {
   const int size = simParams->movingAverageWindowSize;
   totalEnergyAverage.reset(size);
   temperatureAverage.reset(size);
   pressureAverage.reset(size);
   groupPressureAverage.reset(size);
 }

◆ resumeAfterTraceBarrier()

void Controller::resumeAfterTraceBarrier ( int step )

Definition at line 4984 of file Controller.C.

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

Referenced by Node::resumeAfterTraceBarrier().

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

◆ run()

void Controller::run ( void )

Definition at line 358 of file Controller.C.

References awaken(), CTRL_STK_SZ, DebugM, and Controller::CthThreadWrapper::thread.

Referenced by NamdState::runController().

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

◆ stochRescaleCoefficient()

double Controller::stochRescaleCoefficient ( )

Calculate new coefficient for stochastic velocity rescaling and update heat.

Definition at line 1784 of file Controller.C.

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

Referenced by stochRescaleVelocities().

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

◆ stochRescaleVelocities()

void Controller::stochRescaleVelocities ( int step )

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

Generate and broadcast the scale factor for stochastic velocity rescaling.

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

Parameters

step	the current timestep

Definition at line 1768 of file Controller.C.

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

Referenced by integrate().

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

◆ suspend()

void Controller::suspend ( void )

protected

Definition at line 375 of file Controller.C.

                              {
   CthSuspend();
 }

◆ tcoupleVelocities()

void Controller::tcoupleVelocities ( int step )

protected

Definition at line 1746 of file Controller.C.

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

Referenced by integrate().

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

◆ terminate()

void Controller::terminate ( void )

protected

Definition at line 5005 of file Controller.C.

References BackEnd::awaken(), and Controller::CthThreadWrapper::thread.

Referenced by algorithm().

                                {
   BackEnd::awaken();
   CthFree(threadWrapper->thread);
   delete threadWrapper;
   CthSuspend();
 }

◆ traceBarrier()

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

protected

Definition at line 4977 of file Controller.C.

References namdWallTimer().

Referenced by integrate().

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

◆ write_accelMDG_rest_file()

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

protected

Definition at line 2270 of file Controller.C.

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

Referenced by rescaleaccelMD().

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

◆ writeExtendedSystemData()

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

protected

Definition at line 4379 of file Controller.C.

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

Referenced by outputExtendedSystem().

                                                                       {
   Lattice &lattice = state->lattice;
   file.precision(12);
   file << step;
     if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
     if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
     if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
     file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
   if ( simParams->langevinPistonOn ) {
     Vector strainRate = diagonal(langevinPiston_strainRate);
     Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
     file << " " << strainRate.x;
     file << " " << strainRate.y;
     file << " " << strainRate.z;
     file << " " << strainRate2.x;
     file << " " << strainRate2.y;
     file << " " << strainRate2.z;
   }
   if (simParams->monteCarloPressureOn) {
     // We could print the elements as they are, but it would be meaningless,
     // because for constant ratio and isotropic fluctuations, we only operate
     // on x element. 
     if (!simParams->useFlexibleCell) {
       // isotropic fluctuations: use x element for v_x, v_y, and v_z
       file << " " << monteCarloMaxVolume.x;
       file << " " << monteCarloMaxVolume.x;
       file << " " << monteCarloMaxVolume.x;
     } else {
       // we have flexible cell
       if (simParams->useConstantRatio) {
         // use x element for v_x and v_y
         file << " " << monteCarloMaxVolume.x;
         file << " " << monteCarloMaxVolume.x;
         file << " " << monteCarloMaxVolume.z;
       } else {
         // constant area or just flexible in each axis dimension
         file << " " << monteCarloMaxVolume.x;
         file << " " << monteCarloMaxVolume.y;
         file << " " << monteCarloMaxVolume.z;
       }
     }
   }
   file << std::endl;
 }

◆ writeExtendedSystemLabels()

void Controller::writeExtendedSystemLabels ( ofstream_namd & file )

protected

Definition at line 4363 of file Controller.C.

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

Referenced by outputExtendedSystem().

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

◆ writeFepEnergyData()

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

protected

Definition at line 4735 of file Controller.C.

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

Referenced by outputFepEnergy().

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

◆ writeTiEnergyData()

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

protected

Definition at line 4771 of file Controller.C.

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

Referenced by outputTiEnergy().

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

Friends And Related Function Documentation

◆ CheckpointMsg

friend class CheckpointMsg

friend

Definition at line 135 of file Controller.h.

◆ Node

friend class Node

friend

Definition at line 134 of file Controller.h.

◆ ScriptTcl

friend class ScriptTcl

friend

Definition at line 133 of file Controller.h.

Member Data Documentation

◆ accelMDdV

BigReal Controller::accelMDdV

Definition at line 116 of file Controller.h.

Referenced by rescaleaccelMD(), and colvarproxy_namd::update_accelMD_info().

◆ accelMDdVAverage

BigReal Controller::accelMDdVAverage

protected

Definition at line 415 of file Controller.h.

Referenced by rescaleaccelMD().

◆ adaptTempAutoDt

Bool Controller::adaptTempAutoDt

protected

Definition at line 438 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

◆ adaptTempBetaMax

BigReal Controller::adaptTempBetaMax

protected

Definition at line 432 of file Controller.h.

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

◆ adaptTempBetaMin

BigReal Controller::adaptTempBetaMin

protected

Definition at line 431 of file Controller.h.

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

◆ adaptTempBetaN

BigReal* Controller::adaptTempBetaN

protected

Definition at line 427 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

◆ adaptTempBin

int Controller::adaptTempBin

protected

Definition at line 433 of file Controller.h.

Referenced by adaptTempUpdate().

◆ adaptTempBins

int Controller::adaptTempBins

protected

Definition at line 434 of file Controller.h.

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

◆ adaptTempCg

BigReal Controller::adaptTempCg

protected

Definition at line 436 of file Controller.h.

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

◆ adaptTempDBeta

BigReal Controller::adaptTempDBeta

protected

Definition at line 435 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

◆ adaptTempDt

BigReal Controller::adaptTempDt

protected

Definition at line 437 of file Controller.h.

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

◆ adaptTempDTave

BigReal Controller::adaptTempDTave

protected

Definition at line 429 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

◆ adaptTempDTavenum

BigReal Controller::adaptTempDTavenum

protected

Definition at line 430 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

◆ adaptTempDtMax

BigReal Controller::adaptTempDtMax

protected

Definition at line 440 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

◆ adaptTempDtMin

BigReal Controller::adaptTempDtMin

protected

Definition at line 439 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempUpdate().

◆ adaptTempPotEnergyAve

BigReal* Controller::adaptTempPotEnergyAve

protected

Definition at line 424 of file Controller.h.

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

◆ adaptTempPotEnergyAveDen

BigReal* Controller::adaptTempPotEnergyAveDen

protected

Definition at line 422 of file Controller.h.

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

◆ adaptTempPotEnergyAveNum

BigReal* Controller::adaptTempPotEnergyAveNum

protected

Definition at line 421 of file Controller.h.

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

◆ adaptTempPotEnergySamples

int* Controller::adaptTempPotEnergySamples

protected

Definition at line 426 of file Controller.h.

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

◆ adaptTempPotEnergyVar

BigReal* Controller::adaptTempPotEnergyVar

protected

Definition at line 425 of file Controller.h.

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

◆ adaptTempPotEnergyVarNum

BigReal* Controller::adaptTempPotEnergyVarNum

protected

Definition at line 423 of file Controller.h.

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

◆ adaptTempRestartFile

ofstream_namd Controller::adaptTempRestartFile

protected

Definition at line 441 of file Controller.h.

Referenced by adaptTempInit(), and adaptTempWriteRestart().

◆ adaptTempT

BigReal Controller::adaptTempT

protected

Definition at line 428 of file Controller.h.

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

◆ alchWork

BigReal Controller::alchWork

protected

Definition at line 234 of file Controller.h.

Referenced by outputTiEnergy(), and printEnergies().

◆ amd_reduction

RequireReduction* Controller::amd_reduction

protected

Definition at line 354 of file Controller.h.

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

◆ avg_count

int Controller::avg_count

protected

Definition at line 164 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ bondedEnergy_ti_1

BigReal Controller::bondedEnergy_ti_1

protected

Definition at line 209 of file Controller.h.

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

◆ bondedEnergy_ti_2

BigReal Controller::bondedEnergy_ti_2

protected

Definition at line 210 of file Controller.h.

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

◆ bondedEnergyDiff_f

BigReal Controller::bondedEnergyDiff_f

protected

Definition at line 196 of file Controller.h.

Referenced by outputFepEnergy(), and printEnergies().

◆ broadcast

ControllerBroadcasts* Controller::broadcast

protected

◆ checkpoint_lattice

Lattice Controller::checkpoint_lattice

protected

Definition at line 383 of file Controller.h.

Referenced by algorithm().

◆ checkpoint_state

ControllerState Controller::checkpoint_state

protected

Definition at line 384 of file Controller.h.

Referenced by algorithm().

◆ checkpoint_stored

int Controller::checkpoint_stored

protected

Definition at line 382 of file Controller.h.

Referenced by algorithm(), and Controller().

◆ checkpoint_task

int Controller::checkpoint_task

protected

Definition at line 391 of file Controller.h.

Referenced by recvCheckpointAck().

◆ checkpoints

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

protected

Definition at line 390 of file Controller.h.

Referenced by algorithm(), and recvCheckpointReq().

◆ collection

CollectionMaster* const Controller::collection

protected

Definition at line 364 of file Controller.h.

Referenced by enqueueCollections().

◆ computeChecksum

int Controller::computeChecksum

protected

Definition at line 169 of file Controller.h.

Referenced by compareChecksums().

◆ controlNumDegFreedom

int Controller::controlNumDegFreedom

protected

Definition at line 252 of file Controller.h.

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

◆ controlPressure

Tensor Controller::controlPressure

protected

Definition at line 253 of file Controller.h.

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

◆ controlPressure_nbond

Tensor Controller::controlPressure_nbond

protected

Definition at line 157 of file Controller.h.

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

◆ controlPressure_normal

Tensor Controller::controlPressure_normal

protected

Definition at line 156 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ controlPressure_slow

Tensor Controller::controlPressure_slow

protected

Definition at line 158 of file Controller.h.

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

◆ cumAlchWork

BigReal Controller::cumAlchWork

protected

Definition at line 223 of file Controller.h.

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

◆ dE

BigReal Controller::dE

protected

Definition at line 202 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

◆ dG

BigReal Controller::dG

protected

Definition at line 204 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

◆ drudeBondTemp

BigReal Controller::drudeBondTemp

protected

Definition at line 238 of file Controller.h.

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

◆ drudeBondTempAvg

BigReal Controller::drudeBondTempAvg

protected

Definition at line 239 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ electEnergy

BigReal Controller::electEnergy

protected

Definition at line 186 of file Controller.h.

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

◆ electEnergy_f

BigReal Controller::electEnergy_f

protected

Definition at line 197 of file Controller.h.

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

◆ electEnergy_ti_1

BigReal Controller::electEnergy_ti_1

protected

Definition at line 211 of file Controller.h.

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

◆ electEnergy_ti_2

BigReal Controller::electEnergy_ti_2

protected

Definition at line 214 of file Controller.h.

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

◆ electEnergyPME_ti_1

BigReal Controller::electEnergyPME_ti_1

protected

Definition at line 224 of file Controller.h.

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

◆ electEnergyPME_ti_2

BigReal Controller::electEnergyPME_ti_2

protected

Definition at line 225 of file Controller.h.

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

◆ electEnergySlow

BigReal Controller::electEnergySlow

protected

Definition at line 187 of file Controller.h.

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

◆ electEnergySlow_f

BigReal Controller::electEnergySlow_f

protected

Definition at line 198 of file Controller.h.

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

◆ electEnergySlow_ti_1

BigReal Controller::electEnergySlow_ti_1

protected

Definition at line 212 of file Controller.h.

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

◆ electEnergySlow_ti_2

BigReal Controller::electEnergySlow_ti_2

protected

Definition at line 215 of file Controller.h.

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

◆ exp_dE_ByRT

BigReal Controller::exp_dE_ByRT

protected

Definition at line 201 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

◆ fepFile

ofstream_namd Controller::fepFile

protected

Definition at line 372 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

◆ FepNo

int Controller::FepNo

protected

Definition at line 205 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

◆ fepSum

BigReal Controller::fepSum

protected

Definition at line 207 of file Controller.h.

Referenced by outputFepEnergy().

◆ fflush_count

int Controller::fflush_count

protected

Definition at line 329 of file Controller.h.

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

◆ goNativeEnergy

BigReal Controller::goNativeEnergy

protected

Definition at line 192 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

◆ goNonnativeEnergy

BigReal Controller::goNonnativeEnergy

protected

Definition at line 193 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

◆ goTotalEnergy

BigReal Controller::goTotalEnergy

protected

Definition at line 194 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

◆ groGaussEnergy

BigReal Controller::groGaussEnergy

protected

Definition at line 191 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

◆ groLJEnergy

BigReal Controller::groLJEnergy

protected

Definition at line 190 of file Controller.h.

Referenced by printEnergies(), and rescaleaccelMD().

◆ groupPressure

Tensor Controller::groupPressure

protected

Definition at line 251 of file Controller.h.

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

◆ groupPressure_avg

BigReal Controller::groupPressure_avg

protected

Definition at line 163 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ groupPressure_nbond

Tensor Controller::groupPressure_nbond

protected

Definition at line 154 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ groupPressure_normal

Tensor Controller::groupPressure_normal

protected

Definition at line 153 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ groupPressure_slow

Tensor Controller::groupPressure_slow

protected

Definition at line 155 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ groupPressure_tavg

Tensor Controller::groupPressure_tavg

protected

Definition at line 166 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ groupPressureAverage

MovingAverage Controller::groupPressureAverage

protected

Definition at line 450 of file Controller.h.

Referenced by printEnergies(), and resetMovingAverage().

◆ groupPressureAverage_xx

MovingAverage Controller::groupPressureAverage_xx

protected

Definition at line 459 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_xy

MovingAverage Controller::groupPressureAverage_xy

protected

Definition at line 460 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_xz

MovingAverage Controller::groupPressureAverage_xz

protected

Definition at line 461 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_yx

MovingAverage Controller::groupPressureAverage_yx

protected

Definition at line 462 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_yy

MovingAverage Controller::groupPressureAverage_yy

protected

Definition at line 463 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_yz

MovingAverage Controller::groupPressureAverage_yz

protected

Definition at line 464 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_zx

MovingAverage Controller::groupPressureAverage_zx

protected

Definition at line 465 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_zy

MovingAverage Controller::groupPressureAverage_zy

protected

Definition at line 466 of file Controller.h.

Referenced by printEnergies().

◆ groupPressureAverage_zz

MovingAverage Controller::groupPressureAverage_zz

protected

Definition at line 467 of file Controller.h.

Referenced by printEnergies().

◆ heat

BigReal Controller::heat

protected

heat exchanged with the thermostat since firstTimestep

Definition at line 245 of file Controller.h.

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

◆ kineticEnergy

BigReal Controller::kineticEnergy

protected

Definition at line 241 of file Controller.h.

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

◆ kineticEnergyCentered

BigReal Controller::kineticEnergyCentered

protected

Definition at line 243 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

◆ kineticEnergyHalfstep

BigReal Controller::kineticEnergyHalfstep

protected

Definition at line 242 of file Controller.h.

Referenced by printEnergies(), and receivePressure().

◆ langevinPiston_origStrainRate

Tensor Controller::langevinPiston_origStrainRate

protected

Definition at line 289 of file Controller.h.

Referenced by Controller().

◆ ldbSteps

int Controller::ldbSteps

protected

Definition at line 327 of file Controller.h.

Referenced by rebalanceLoad().

◆ ljEnergy

BigReal Controller::ljEnergy

protected

Definition at line 188 of file Controller.h.

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

◆ ljEnergy_f

BigReal Controller::ljEnergy_f

protected

Definition at line 199 of file Controller.h.

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

◆ ljEnergy_f_left

BigReal Controller::ljEnergy_f_left

protected

Definition at line 200 of file Controller.h.

◆ ljEnergy_ti_1

BigReal Controller::ljEnergy_ti_1

protected

Definition at line 213 of file Controller.h.

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

◆ ljEnergy_ti_2

BigReal Controller::ljEnergy_ti_2

protected

Definition at line 216 of file Controller.h.

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

◆ ljEnergySlow

BigReal Controller::ljEnergySlow

protected

Definition at line 189 of file Controller.h.

Referenced by getTotalPotentialEnergy(), and printEnergies().

◆ marginViolations

int Controller::marginViolations

protected

Definition at line 170 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

◆ mc_accept

int Controller::mc_accept[MC_AXIS_TOTAL]

protected

Definition at line 309 of file Controller.h.

Referenced by Controller(), and monteCarloPressure_accept().

◆ mc_oldLattice

Lattice Controller::mc_oldLattice

protected

Definition at line 313 of file Controller.h.

Referenced by monteCarloPressure_accept(), and monteCarloPressure_prepare().

◆ mc_picked_axis

int Controller::mc_picked_axis

protected

Definition at line 311 of file Controller.h.

Referenced by monteCarloPressure_accept(), and monteCarloPressure_prepare().

◆ mc_totalAccept

int Controller::mc_totalAccept

protected

Definition at line 310 of file Controller.h.

Referenced by Controller(), and monteCarloPressure_accept().

◆ mc_totalEnergyOld

BigReal Controller::mc_totalEnergyOld

protected

Definition at line 312 of file Controller.h.

Referenced by monteCarloPressure_accept(), and monteCarloPressure_prepare().

◆ mc_totalTry

int Controller::mc_totalTry

protected

Definition at line 310 of file Controller.h.

Referenced by Controller(), monteCarloPressure_accept(), and monteCarloPressure_prepare().

◆ mc_trial

int Controller::mc_trial[MC_AXIS_TOTAL]

protected

Definition at line 309 of file Controller.h.

Referenced by Controller(), monteCarloPressure_accept(), and monteCarloPressure_prepare().

◆ min_energy

BigReal Controller::min_energy

protected

Definition at line 176 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

◆ min_f_dot_f

BigReal Controller::min_f_dot_f

protected

Definition at line 177 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

◆ min_f_dot_v

BigReal Controller::min_f_dot_v

protected

Definition at line 178 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

◆ min_huge_count

int Controller::min_huge_count

protected

Definition at line 180 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

◆ min_reduction

RequireReduction* Controller::min_reduction

protected

Definition at line 140 of file Controller.h.

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

◆ min_v_dot_v

BigReal Controller::min_v_dot_v

protected

Definition at line 179 of file Controller.h.

Referenced by minimize(), and printMinimizeEnergies().

◆ momentumSqrSum

Tensor Controller::momentumSqrSum

protected

Definition at line 318 of file Controller.h.

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

◆ multigratorNu

std::vector<BigReal> Controller::multigratorNu

protected

Definition at line 320 of file Controller.h.

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

◆ multigratorNuT

std::vector<BigReal> Controller::multigratorNuT

protected

Definition at line 321 of file Controller.h.

Referenced by Controller(), and multigratorTemperature().

◆ multigratorOmega

std::vector<BigReal> Controller::multigratorOmega

protected

Definition at line 322 of file Controller.h.

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

◆ multigratorReduction

RequireReduction* Controller::multigratorReduction

protected

Definition at line 324 of file Controller.h.

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

◆ multigratorXi

BigReal Controller::multigratorXi

protected

Definition at line 316 of file Controller.h.

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

◆ multigratorXiT

BigReal Controller::multigratorXiT

protected

Definition at line 317 of file Controller.h.

Referenced by multigratorPressure().

◆ multigratorZeta

std::vector<BigReal> Controller::multigratorZeta

protected

Definition at line 323 of file Controller.h.

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

◆ nbondFreq

int Controller::nbondFreq

protected

Definition at line 159 of file Controller.h.

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

◆ net_dE

BigReal Controller::net_dE

protected

Definition at line 203 of file Controller.h.

Referenced by outputFepEnergy(), and writeFepEnergyData().

◆ net_dEdl_bond_1

BigReal Controller::net_dEdl_bond_1

protected

Definition at line 217 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ net_dEdl_bond_2

BigReal Controller::net_dEdl_bond_2

protected

Definition at line 218 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ net_dEdl_elec_1

BigReal Controller::net_dEdl_elec_1

protected

Definition at line 219 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ net_dEdl_elec_2

BigReal Controller::net_dEdl_elec_2

protected

Definition at line 220 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ net_dEdl_lj_1

BigReal Controller::net_dEdl_lj_1

protected

Definition at line 221 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ net_dEdl_lj_2

BigReal Controller::net_dEdl_lj_2

protected

Definition at line 222 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ numDegFreedom

int64_t Controller::numDegFreedom

protected

Definition at line 183 of file Controller.h.

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

◆ origLattice

Lattice Controller::origLattice

protected

Definition at line 395 of file Controller.h.

Referenced by Controller().

◆ pairlistWarnings

int Controller::pairlistWarnings

protected

Definition at line 171 of file Controller.h.

Referenced by compareChecksums(), and printEnergies().

◆ perfstats

RunningAverage Controller::perfstats

protected

Definition at line 444 of file Controller.h.

Referenced by printTiming().

◆ positionRescaleFactor

Tensor Controller::positionRescaleFactor

protected

Definition at line 291 of file Controller.h.

Referenced by langevinPiston1().

◆ ppbonded

PressureProfileReduction* Controller::ppbonded

protected

Definition at line 358 of file Controller.h.

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

◆ ppint

PressureProfileReduction* Controller::ppint

protected

Definition at line 360 of file Controller.h.

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

◆ ppnonbonded

PressureProfileReduction* Controller::ppnonbonded

protected

Definition at line 359 of file Controller.h.

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

◆ pressure

Tensor Controller::pressure

protected

Definition at line 250 of file Controller.h.

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

◆ pressure_amd

Tensor Controller::pressure_amd

protected

Definition at line 151 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ pressure_avg

BigReal Controller::pressure_avg

protected

Definition at line 162 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ pressure_nbond

Tensor Controller::pressure_nbond

protected

Definition at line 149 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ pressure_normal

Tensor Controller::pressure_normal

protected

Definition at line 148 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ pressure_slow

Tensor Controller::pressure_slow

protected

Definition at line 150 of file Controller.h.

Referenced by calcPressure(), and receivePressure().

◆ pressure_tavg

Tensor Controller::pressure_tavg

protected

Definition at line 165 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ pressureAverage

MovingAverage Controller::pressureAverage

protected

Definition at line 449 of file Controller.h.

Referenced by printEnergies(), and resetMovingAverage().

◆ pressureAverage_xx

MovingAverage Controller::pressureAverage_xx

protected

Definition at line 453 of file Controller.h.

Referenced by printEnergies().

◆ pressureAverage_yx

MovingAverage Controller::pressureAverage_yx

protected

Definition at line 454 of file Controller.h.

Referenced by printEnergies().

◆ pressureAverage_yy

MovingAverage Controller::pressureAverage_yy

protected

Definition at line 455 of file Controller.h.

Referenced by printEnergies().

◆ pressureAverage_zx

MovingAverage Controller::pressureAverage_zx

protected

Definition at line 456 of file Controller.h.

Referenced by printEnergies().

◆ pressureAverage_zy

MovingAverage Controller::pressureAverage_zy

protected

Definition at line 457 of file Controller.h.

Referenced by printEnergies().

◆ pressureAverage_zz

MovingAverage Controller::pressureAverage_zz

protected

Definition at line 458 of file Controller.h.

Referenced by printEnergies().

◆ pressureProfileAverage

BigReal* Controller::pressureProfileAverage

protected

Definition at line 363 of file Controller.h.

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

◆ pressureProfileCount

int Controller::pressureProfileCount

protected

Definition at line 362 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ pressureProfileSlabs

int Controller::pressureProfileSlabs

protected

Definition at line 361 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ random

Random* Controller::random

protected

Definition at line 341 of file Controller.h.

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

◆ recent_alchWork

BigReal Controller::recent_alchWork

protected

Definition at line 233 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ recent_dEdl_bond_1

BigReal Controller::recent_dEdl_bond_1

protected

Definition at line 227 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ recent_dEdl_bond_2

BigReal Controller::recent_dEdl_bond_2

protected

Definition at line 228 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ recent_dEdl_elec_1

BigReal Controller::recent_dEdl_elec_1

protected

Definition at line 229 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ recent_dEdl_elec_2

BigReal Controller::recent_dEdl_elec_2

protected

Definition at line 230 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ recent_dEdl_lj_1

BigReal Controller::recent_dEdl_lj_1

protected

Definition at line 231 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ recent_dEdl_lj_2

BigReal Controller::recent_dEdl_lj_2

protected

Definition at line 232 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ recent_TiNo

int Controller::recent_TiNo

protected

Definition at line 235 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ reductionBasic

RequireReduction* Controller::reductionBasic

protected

Definition at line 344 of file Controller.h.

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

◆ reductionGpuResident

RequireReduction* Controller::reductionGpuResident

protected

Definition at line 352 of file Controller.h.

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

◆ rescaleVelocities_numTemps

int Controller::rescaleVelocities_numTemps

protected

Definition at line 258 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

◆ rescaleVelocities_sumTemps

BigReal Controller::rescaleVelocities_sumTemps

protected

Definition at line 257 of file Controller.h.

Referenced by Controller(), and rescaleVelocities().

◆ simParams

SimParameters* const Controller::simParams

protected

Definition at line 342 of file Controller.h.

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

◆ slowFreq

int Controller::slowFreq

protected

Definition at line 160 of file Controller.h.

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

◆ smooth2_avg2

BigReal Controller::smooth2_avg2

protected

Definition at line 249 of file Controller.h.

◆ state

NamdState* const Controller::state

protected

◆ stepInFullRun

int Controller::stepInFullRun

protected

Definition at line 184 of file Controller.h.

Referenced by printEnergies().

◆ stochRescale_count

int Controller::stochRescale_count

Count time steps until next stochastic velocity rescaling.

Definition at line 276 of file Controller.h.

Referenced by Controller(), and stochRescaleVelocities().

◆ stochRescaleTimefactor

BigReal Controller::stochRescaleTimefactor

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

Definition at line 279 of file Controller.h.

Referenced by Controller(), and stochRescaleCoefficient().

◆ strainRate_old

Tensor Controller::strainRate_old

protected

Definition at line 290 of file Controller.h.

Referenced by langevinPiston1().

◆ submit_reduction

SubmitReduction* Controller::submit_reduction

protected

Definition at line 355 of file Controller.h.

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

◆ tavg_count

int Controller::tavg_count

protected

Definition at line 167 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ temp_avg

BigReal Controller::temp_avg

protected

Definition at line 161 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ temperature

BigReal Controller::temperature

protected

Definition at line 244 of file Controller.h.

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

◆ temperatureAverage

MovingAverage Controller::temperatureAverage

protected

Definition at line 448 of file Controller.h.

Referenced by printEnergies(), and resetMovingAverage().

◆ tiFile

ofstream_namd Controller::tiFile

protected

Definition at line 376 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ TiNo

int Controller::TiNo

protected

Definition at line 226 of file Controller.h.

Referenced by outputTiEnergy(), and writeTiEnergyData().

◆ totalEnergy

BigReal Controller::totalEnergy

protected

Definition at line 185 of file Controller.h.

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

◆ totalEnergy0

BigReal Controller::totalEnergy0

protected

totalEnergy at firstTimestep

Definition at line 247 of file Controller.h.

Referenced by Controller(), and printEnergies().

◆ totalEnergyAverage

MovingAverage Controller::totalEnergyAverage

protected

Definition at line 447 of file Controller.h.

Referenced by printEnergies(), and resetMovingAverage().

◆ virial_amd

Tensor Controller::virial_amd

protected

Definition at line 152 of file Controller.h.

Referenced by calcPressure(), and rescaleaccelMD().

◆ xstFile

ofstream_namd Controller::xstFile

protected

Definition at line 366 of file Controller.h.

Referenced by outputExtendedSystem().

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

Classes

Public Member Functions

Public Attributes

Protected Types

Protected Member Functions

Protected Attributes

Friends

Detailed Description

Member Enumeration Documentation

◆ mc_axis_pick

Constructor & Destructor Documentation

◆ Controller()

◆ ~Controller()

Member Function Documentation

◆ adaptTempInit()

◆ adaptTempUpdate()

◆ adaptTempWriteRestart()

◆ algorithm()

◆ awaken()

◆ berendsenPressure()

◆ calc_accelMDG_E_k()

◆ calc_accelMDG_force_factor()

◆ calc_accelMDG_mean_std()

◆ calcPressure()

◆ compareChecksums()

◆ computeAlchWork()

◆ correctMomentum()

◆ cycleBarrier()

◆ enqueueCollections()

◆ getTIderivative()

◆ getTotalPotentialEnergy()

◆ integrate()

◆ langevinPiston1()

◆ langevinPiston2()

◆ minimize()

◆ monteCarloPressure_accept()

◆ monteCarloPressure_prepare()

◆ multigatorCalcEnthalpy()

◆ multigratorPressure()

◆ multigratorTemperature()

◆ outputExtendedSystem()

◆ outputFepEnergy()

◆ outputTiEnergy()

◆ printDynamicsEnergies()

◆ printEnergies()

◆ printFepMessage()

◆ printMinimizeEnergies()

◆ printTiMessage()

◆ printTiming()

◆ reassignVelocities()

◆ rebalanceLoad()

◆ receivePressure()

◆ recvCheckpointAck()

◆ recvCheckpointReq()

◆ rescaleaccelMD()

◆ rescaleVelocities()

◆ resetMovingAverage()

◆ resumeAfterTraceBarrier()

◆ run()

◆ stochRescaleCoefficient()

◆ stochRescaleVelocities()

◆ suspend()

◆ tcoupleVelocities()

◆ terminate()

◆ traceBarrier()

◆ write_accelMDG_rest_file()

◆ writeExtendedSystemData()

◆ writeExtendedSystemLabels()

◆ writeFepEnergyData()

◆ writeTiEnergyData()

Friends And Related Function Documentation

◆ CheckpointMsg

◆ Node

◆ ScriptTcl

Member Data Documentation

◆ accelMDdV

◆ accelMDdVAverage

◆ adaptTempAutoDt

◆ adaptTempBetaMax

◆ adaptTempBetaMin