Simbody
3.7

This is a feedbackcontrolled force that uses Nose'Hoover chains to maintain a particular temperature Tb, as though the system were immersed in an infinite heat bath at that temperature. More...
Public Member Functions  
Thermostat (GeneralForceSubsystem &forces, const SimbodyMatterSubsystem &matter, Real boltzmannsConstant, Real bathTemperature, Real relaxationTime, int numExcludedDofs=6)  
Define a global thermostat (one that affects all degrees of freedom) at a given default temperature and relaxation time. More...  
Thermostat ()  
Default constructor creates an empty handle. More...  
Thermostat &  excludeMobilizedBody (MobilizedBodyIndex) 
TODO: not implemented yet. More...  
Thermostat &  setDefaultNumChains (int numChains) 
Set the default (state independent) number of Nose'Hoover chains. More...  
Thermostat &  setDefaultBathTemperature (Real bathTemperature) 
Set the default (state independent) bath temperature. More...  
Thermostat &  setDefaultRelaxationTime (Real relaxationTime) 
Set the default (state independent) relaxation time. More...  
Thermostat &  setDefaultNumExcludedDofs (int numExcludedDofs) 
Set the default number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N; if you don't call this it is assumed that 6 dofs should be excluded. More...  
int  getDefaultNumChains () const 
Get the initial value for the number of chains that will be used for the "number of chains" State variable. More...  
Real  getDefaultBathTemperature () const 
Get the initial value for the bath temperature that will be use for the "bath temperature" State variable. More...  
Real  getDefaultRelaxationTime () const 
Get the initial value for the bath temperature that will be use for the "bath temperature" State variable. More...  
int  getDefaultNumExcludedDofs () const 
Get the initial value for the number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N. More...  
Real  getBoltzmannsConstant () const 
Can't change the value of Boltzmann's constant after construction; this is the value being used. More...  
const Thermostat &  setNumChains (State &, int numChains) const 
Set the actual number of Nose'Hoover chains to be used. More...  
const Thermostat &  setBathTemperature (State &, Real Tb) const 
Set the bath temperature which serves as the target temperature for the thermostat. More...  
const Thermostat &  setRelaxationTime (State &, Real t) const 
Set the relaxation time which determines how long the system will take to equilibrate to the bath temperature. More...  
const Thermostat &  setNumExcludedDofs (State &, int numExcludedDofs) const 
Set the actual number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N. More...  
int  getNumChains (const State &) const 
Obtain the current number of Nose'Hoover chains in use. More...  
Real  getBathTemperature (const State &) const 
Obtain the current bath temperature, in units which are determined by the value of Boltzmann's constant that was supplied on construction of this Thermostat force element. More...  
Real  getRelaxationTime (const State &) const 
Obtain the current relaxation time. More...  
int  getNumExcludedDofs (const State &) const 
Get the current value for the number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N. More...  
int  getNumThermalDofs (const State &) const 
Return the number of thermal degrees of freedom being used in the definition of temperature for this thermostat. More...  
Real  getCurrentTemperature (const State &) const 
Return the temperature of the controlled degrees of freedom via the definition T = 2*ke / (N*Kb) where N is the number of thermal degrees of freedom. More...  
void  initializeChainState (State &) const 
This is a solver that initializes thermostat state variables to zero. More...  
void  setChainState (State &, const Vector &) const 
Set the thermostat state variables to particular values. More...  
Vector  getChainState (const State &) const 
Return the current values of the thermostat chain variables. More...  
Real  calcBathEnergy (const State &state) const 
Calculate the total "bath energy" which, when added to the system energy, should yield a conserved quantity (assuming all other forces are conservative). More...  
Real  getExternalPower (const State &state) const 
Get the amount of power being applied by the thermostat to the system; sign is positive when energy is coming from the bath. More...  
Real  getExternalWork (const State &state) const 
Get the amount of work that has been done by the bath on the system since an arbitrary start time. More...  
void  setExternalWork (State &state, Real work) const 
Set the current value of the work done by the bath to an arbitrary value; normally zero for initialization. More...  
void  initializeSystemToBathTemperature (State &) const 
Set the controlled system to a set of randomized velocities which yields the bath temperature. More...  
void  setSystemToTemperature (State &, Real T) const 
Set the controlled system to a set of randomized velocities which yields a particular temperature. More...  
Public Member Functions inherited from SimTK::Force  
void  disable (State &) const 
Disable this force element, effectively removing it from the System for computational purposes (it is still using its ForceIndex, however). More...  
void  enable (State &) const 
Enable this force element if it was previously disabled. More...  
bool  isDisabled (const State &) const 
Test whether this force element is currently disabled in the supplied State. More...  
void  setDisabledByDefault (bool shouldBeDisabled) 
Normally force elements are enabled when defined and can be disabled later. More...  
bool  isDisabledByDefault () const 
Test whether this force element is disabled by default in which case it must be explicitly enabled before it will take effect. More...  
void  calcForceContribution (const State &state, Vector_< SpatialVec > &bodyForces, Vector_< Vec3 > &particleForces, Vector &mobilityForces) const 
Calculate the force that would be applied by this force element if the given state were realized to Dynamics stage. More...  
Real  calcPotentialEnergyContribution (const State &state) const 
Calculate the potential energy contribution that is made by this force element at the given state. More...  
Force ()  
Default constructor for Force handle base class does nothing. More...  
operator ForceIndex () const  
Implicit conversion to ForceIndex when needed. More...  
const GeneralForceSubsystem &  getForceSubsystem () const 
Get the GeneralForceSubsystem of which this Force is an element. More...  
ForceIndex  getForceIndex () const 
Get the index of this force element within its parent force subsystem. More...  
Public Member Functions inherited from SimTK::PIMPLHandle< Force, ForceImpl, true >  
bool  isEmptyHandle () const 
Returns true if this handle is empty, that is, does not refer to any implementation object. More...  
bool  isOwnerHandle () const 
Returns true if this handle is the owner of the implementation object to which it refers. More...  
bool  isSameHandle (const Force &other) const 
Determine whether the supplied handle is the same object as "this" PIMPLHandle. More...  
void  disown (Force &newOwner) 
Give up ownership of the implementation to an empty handle. More...  
PIMPLHandle &  referenceAssign (const Force &source) 
"Copy" assignment but with shallow (pointer) semantics. More...  
PIMPLHandle &  copyAssign (const Force &source) 
This is real copy assignment, with ordinary C++ object ("value") semantics. More...  
void  clearHandle () 
Make this an empty handle, deleting the implementation object if this handle is the owner of it. More...  
const ForceImpl &  getImpl () const 
Get a const reference to the implementation associated with this Handle. More...  
ForceImpl &  updImpl () 
Get a writable reference to the implementation associated with this Handle. More...  
int  getImplHandleCount () const 
Return the number of handles the implementation believes are referencing it. More...  
Additional Inherited Members  
Public Types inherited from SimTK::PIMPLHandle< Force, ForceImpl, true >  
typedef PIMPLHandle< Force, ForceImpl, PTR >  HandleBase 
typedef HandleBase  ParentHandle 
Protected Member Functions inherited from SimTK::Force  
Force (ForceImpl *r)  
Use this in a derived Force handle class constructor to supply the concrete implementation object to be stored in the handle base. More...  
Protected Member Functions inherited from SimTK::PIMPLHandle< Force, ForceImpl, true >  
PIMPLHandle ()  
The default constructor makes this an empty handle. More...  
PIMPLHandle (ForceImpl *p)  
This provides consruction of a handle referencing an existing implementation object. More...  
PIMPLHandle (const PIMPLHandle &source)  
The copy constructor makes either a deep (value) or shallow (reference) copy of the supplied source PIMPL object, based on whether this is a "pointer
semantics" (PTR=true) or "object (value) semantics" (PTR=false, default) class. More...  
~PIMPLHandle ()  
Note that the destructor is nonvirtual. More...  
PIMPLHandle &  operator= (const PIMPLHandle &source) 
Copy assignment makes the current handle either a deep (value) or shallow (reference) copy of the supplied source PIMPL object, based on whether this is a "pointer sematics" (PTR=true) or "object (value) semantics" (PTR=false, default) class. More...  
void  setImpl (ForceImpl *p) 
Set the implementation for this empty handle. More...  
bool  hasSameImplementation (const Force &other) const 
Determine whether the supplied handle is a reference to the same implementation object as is referenced by "this" PIMPLHandle. More...  
This is a feedbackcontrolled force that uses Nose'Hoover chains to maintain a particular temperature Tb, as though the system were immersed in an infinite heat bath at that temperature.
There are two parameters, the temperature Tb and a "relaxation time" t which controls how tightly the temperature is maintained. This thermostat is particularly useful in molecular simulations but can be applied to any mechanical system also.
Temperature is defined here as T = (2*KE) / (N*kB) where KE is the system kinetic energy, N is the number of coupled degrees of freedom (mobilities minus active, nonredundant constraints, minus up to 6 rigid body dofs for the system as a whole), and kB is Boltzmann's constant in appropriate units.
We use a Nose'Hoover chain to achieve excellent statistical mechanics properties with a continuous force. At equilibrium the temperature will have a Boltzmann distribution; the relaxation time controls how long it takes the system to reach equilibrium with the bath. Smaller values of relaxation time produce faster response but can make the system stiff and will normally require smaller step sizes; larger values will take longer to equilibrate but will run faster.
This Force does not produce any potential energy. However, there is a "bath energy" available through a separate call which can be used in combination with the system energy to construct a conserved quantity; this is described further below.
The current system temperature is defined
T = (2*KE) / (N*kB)
where KE is the kinetic energy of the moving bodies whose N degrees of freedom are being controlled (not necessarily all the bodies in the system), and kB is Boltzmann's constant. Our goal here is to control T so that it follows a Boltzmann distribution around the specified bath temperature Tb.
For an mchain Nose'Hoover chain, we will define m auxiliary "thermostat" state variables c[i], 0<=i<m, with units of 1/time. The 0'th thermostat variable c[0] is used to generate a generalized force f applied to the system mobilities u:
f = c[0] * M * u
where M is the system mass matrix and u is the vector of generalized speeds. (Note that in Simbody the M*u product is formed in O(n) time; M itself is never formed.) The c variables should be initialized to zero at the start of a simulation. Ideally, you should initialize the u's so that they are already at the right temperature, but if not you should still make them nonzero – you can see above that if you have no velocities you will get no Nose'Hoover forces.
If m==1, we have the standard Nose'Hoover method except with a relaxation time specified instead of the thermal mass parameter, as in reference [2]:
cdot[0] = (T/Tb  1)/t^2
Otherwise, for m > 1 we have:
cdot[0] = (T/Tb  1)/t^2  c[0]*c[1] cdot[1] = N*c[0]^2  1/t^2  c[1]*c[2] cdot[i] = c[i1]^2  1/t^2  c[i]*c[i+1] (2<=i<m1) cdot[m1] = c[m2]^2  1/t^2
For comparison with the literature where thermal mass parameters Qi are used, we use Q0 = N kB Tb t^2 and Qi = kB Tb t^2, i > 0. That is, the first thermostat that controls the N thermal degrees of freedom is N times "heavier" than the subsequent ones, each of which controls only the one dof of its nextlower thermostat. See refs [1] and [2].
In addition there is a set of state variables si given by sdot[i]=c[i]. Together these permit us to define a "bath energy" which can be combined with system energy to produce a conserved quantity. Bath energy is KEb+PEb where
KEb = 1/2 kB Tb t^2 (N c[0]^2 + sum(c[i]^2)) PEb = kB Tb (N s[0] + sum(s[i]))
where kB is Boltzmann's constant, Tb the bath temperature, N the number of thermal degrees of freedom in the temperature definition, and the sums run from 1 to m1. Note that you must request the bath energy separately; we do not return any potential energy for this force otherwise.
[1] Martyna, GJ; Klien, ML; Tuckerman, M. Nose'Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 97(4):26352643 (1992).
[2] Vaidehi, N; Jain, A; Goddard, WA. Constant Temperature Constrained Molecular Dynamics: The NewtonEuler Inverse Mass Operator Method. J. Phys. Chem. 100:1050810517 (1996).
SimTK::Force::Thermostat::Thermostat  (  GeneralForceSubsystem &  forces, 
const SimbodyMatterSubsystem &  matter,  
Real  boltzmannsConstant,  
Real  bathTemperature,  
Real  relaxationTime,  
int  numExcludedDofs = 6 

) 
Define a global thermostat (one that affects all degrees of freedom) at a given default temperature and relaxation time.
The number of Nose'Hoover chains is given a default value.

inline 
Default constructor creates an empty handle.
Thermostat& SimTK::Force::Thermostat::excludeMobilizedBody  (  MobilizedBodyIndex  ) 
TODO: not implemented yet.
Remove a body from consideration in the thermostat. Typically this would be the system base body so that overall rigid body translation and orientation is not counted as part of the temperature.
Thermostat& SimTK::Force::Thermostat::setDefaultNumChains  (  int  numChains  ) 
Set the default (state independent) number of Nose'Hoover chains.
This is a Topologystage change.
Thermostat& SimTK::Force::Thermostat::setDefaultBathTemperature  (  Real  bathTemperature  ) 
Set the default (state independent) bath temperature.
This will be interpreted using the value of Boltzmann's constant Kb provided on construction. The units will be Kb/energy, typically Kelvins.
Thermostat& SimTK::Force::Thermostat::setDefaultRelaxationTime  (  Real  relaxationTime  ) 
Set the default (state independent) relaxation time.
Thermostat& SimTK::Force::Thermostat::setDefaultNumExcludedDofs  (  int  numExcludedDofs  ) 
Set the default number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N; if you don't call this it is assumed that 6 dofs should be excluded.
int SimTK::Force::Thermostat::getDefaultNumChains  (  )  const 
Real SimTK::Force::Thermostat::getDefaultBathTemperature  (  )  const 
Real SimTK::Force::Thermostat::getDefaultRelaxationTime  (  )  const 
Get the initial value for the bath temperature that will be use for the "bath temperature" State variable.
int SimTK::Force::Thermostat::getDefaultNumExcludedDofs  (  )  const 
Get the initial value for the number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N.
A new value may be set in a particular State.
Real SimTK::Force::Thermostat::getBoltzmannsConstant  (  )  const 
Can't change the value of Boltzmann's constant after construction; this is the value being used.
const Thermostat& SimTK::Force::Thermostat::setNumChains  (  State &  , 
int  numChains  
)  const 
Set the actual number of Nose'Hoover chains to be used.
This variable controls the number of auxiliary state variables allocated by the Thermostat so invalidates Model stage.
const Thermostat& SimTK::Force::Thermostat::setBathTemperature  (  State &  , 
Real  Tb  
)  const 
Set the bath temperature which serves as the target temperature for the thermostat.
This is given in units defined by the value of Boltzmann's constant (which has units of energy/temperature) that was set on construction. This sets an Instancestage state variable so invalidates Instance and higher stages in the given State.
const Thermostat& SimTK::Force::Thermostat::setRelaxationTime  (  State &  , 
Real  t  
)  const 
Set the relaxation time which determines how long the system will take to equilibrate to the bath temperature.
This sets an Instancestage state variable so invalidates Instance and higher stages in the given State.
const Thermostat& SimTK::Force::Thermostat::setNumExcludedDofs  (  State &  , 
int  numExcludedDofs  
)  const 
Set the actual number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N.
This sets an Instancestage state variable so invalidates Instance and higher stages in the given State.
int SimTK::Force::Thermostat::getNumChains  (  const State &  )  const 
Obtain the current number of Nose'Hoover chains in use.
This is a state variable so can be obtained any time after realization of the Model stage.
Obtain the current bath temperature, in units which are determined by the value of Boltzmann's constant that was supplied on construction of this Thermostat force element.
This is a state variable so can be obtained any time after realization of the Model stage.
Obtain the current relaxation time.
This is a state variable so can be obtained any time after realization of the Model stage.
int SimTK::Force::Thermostat::getNumExcludedDofs  (  const State &  )  const 
Get the current value for the number of system rigid body degrees of freedom (06) to be excluded from the calculation of the number of thermal degrees of freedom N.
This is a state variable so can be obtained any time after realization of the Model stage.
int SimTK::Force::Thermostat::getNumThermalDofs  (  const State &  )  const 
Return the number of thermal degrees of freedom being used in the definition of temperature for this thermostat.
This is the net of the total number of mobilities minus nonredundant constraints minus the number of excluded system rigid body degrees of freedom (06).
Return the temperature of the controlled degrees of freedom via the definition T = 2*ke / (N*Kb) where N is the number of thermal degrees of freedom.
You can call this after Stage::Velocity has been realized.
void SimTK::Force::Thermostat::initializeChainState  (  State &  )  const 
This is a solver that initializes thermostat state variables to zero.
Set the thermostat state variables to particular values.
The Vector's length must be the same as twice the current number of chains called for by the State.
Return the current values of the thermostat chain variables.
The returned vector will have twice the length that getNumChains() would return if called on this same State.
Calculate the total "bath energy" which, when added to the system energy, should yield a conserved quantity (assuming all other forces are conservative).
Get the amount of power being applied by the thermostat to the system; sign is positive when energy is coming from the bath.
Get the amount of work that has been done by the bath on the system since an arbitrary start time.
Set the current value of the work done by the bath to an arbitrary value; normally zero for initialization.
void SimTK::Force::Thermostat::initializeSystemToBathTemperature  (  State &  )  const 
Set the controlled system to a set of randomized velocities which yields the bath temperature.
This ignores the current system velocities. TODO: not implemented yet.
Set the controlled system to a set of randomized velocities which yields a particular temperature.
This ignores the current system velocities. The temperature is interpreted using the value of Boltzmann's constant that was provided on construction of this Thermostat. TODO: not implemented yet.