Simbody  3.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SimTK::Force::Thermostat Class Reference

This is a feedback-controlled 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...

#include <Force_Thermostat.h>

+ Inheritance diagram for SimTK::Force::Thermostat:

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...
 
ThermostatexcludeMobilizedBody (MobilizedBodyIndex)
 TODO: not implemented yet. More...
 
ThermostatsetDefaultNumChains (int numChains)
 Set the default (state independent) number of Nose'-Hoover chains. More...
 
ThermostatsetDefaultBathTemperature (Real bathTemperature)
 Set the default (state independent) bath temperature. More...
 
ThermostatsetDefaultRelaxationTime (Real relaxationTime)
 Set the default (state independent) relaxation time. More...
 
ThermostatsetDefaultNumExcludedDofs (int numExcludedDofs)
 Set the default number of system rigid body degrees of freedom (0-6) 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 (0-6) 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 ThermostatsetNumChains (State &, int numChains) const
 Set the actual number of Nose'-Hoover chains to be used. More...
 
const ThermostatsetBathTemperature (State &, Real Tb) const
 Set the bath temperature which serves as the target temperature for the thermostat. More...
 
const ThermostatsetRelaxationTime (State &, Real t) const
 Set the relaxation time which determines how long the system will take to equilibrate to the bath temperature. More...
 
const ThermostatsetNumExcludedDofs (State &, int numExcludedDofs) const
 Set the actual number of system rigid body degrees of freedom (0-6) 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 (0-6) 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 GeneralForceSubsystemgetForceSubsystem () 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...
 
PIMPLHandlereferenceAssign (const Force &source)
 "Copy" assignment but with shallow (pointer) semantics. More...
 
PIMPLHandlecopyAssign (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 non-virtual. More...
 
PIMPLHandleoperator= (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...
 

Detailed Description

This is a feedback-controlled 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.

Theory:

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 m-chain 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 non-zero – 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[i-1]^2 - 1/t^2 - c[i]*c[i+1]   (2<=i<m-1)
     cdot[m-1] = c[m-2]^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 next-lower 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 m-1. Note that you must request the bath energy separately; we do not return any potential energy for this force otherwise.

References:

[1] Martyna, GJ; Klien, ML; Tuckerman, M. Nose'-Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 97(4):2635-2643 (1992).

[2] Vaidehi, N; Jain, A; Goddard, WA. Constant Temperature Constrained Molecular Dynamics: The Newton-Euler Inverse Mass Operator Method. J. Phys. Chem. 100:10508-10517 (1996).

Constructor & Destructor Documentation

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.

Default constructor creates an empty handle.

Member Function Documentation

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 Topology-stage 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 (0-6) 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

Get the initial value for the number of chains that will be used for the "number of chains" State variable.

A new value may be set in a particular State.

Real SimTK::Force::Thermostat::getDefaultBathTemperature ( ) const

Get the initial value for the bath temperature that will be use for the "bath temperature" State variable.

A new value may be set in a particular State.

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 (0-6) 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 Instance-stage 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 Instance-stage 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 (0-6) to be excluded from the calculation of the number of thermal degrees of freedom N.

This sets an Instance-stage 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.

Real SimTK::Force::Thermostat::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.

This is a state variable so can be obtained any time after realization of the Model stage.

Real SimTK::Force::Thermostat::getRelaxationTime ( const State ) const

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 (0-6) 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 (0-6).

Real SimTK::Force::Thermostat::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.

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.

void SimTK::Force::Thermostat::setChainState ( State ,
const Vector  
) const

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.

Vector SimTK::Force::Thermostat::getChainState ( const State ) const

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.

Real SimTK::Force::Thermostat::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).

Real SimTK::Force::Thermostat::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.

Real SimTK::Force::Thermostat::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.

void SimTK::Force::Thermostat::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.

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.

void SimTK::Force::Thermostat::setSystemToTemperature ( State ,
Real  T 
) const

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.


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