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...
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... | |
| 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 (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 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 (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 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 non-virtual. 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... | |
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
◆ Thermostat() [1/2]
| 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.
◆ Thermostat() [2/2]
|
inline |
Default constructor creates an empty handle.
Member Function Documentation
◆ excludeMobilizedBody()
| 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.
◆ setDefaultNumChains()
| Thermostat& SimTK::Force::Thermostat::setDefaultNumChains | ( | int | numChains | ) |
Set the default (state independent) number of Nose'-Hoover chains.
This is a Topology-stage change.
◆ setDefaultBathTemperature()
| 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.
◆ setDefaultRelaxationTime()
| Thermostat& SimTK::Force::Thermostat::setDefaultRelaxationTime | ( | Real | relaxationTime | ) |
Set the default (state independent) relaxation time.
◆ setDefaultNumExcludedDofs()
| 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.
◆ getDefaultNumChains()
| int SimTK::Force::Thermostat::getDefaultNumChains | ( | ) | const |
◆ getDefaultBathTemperature()
| Real SimTK::Force::Thermostat::getDefaultBathTemperature | ( | ) | const |
◆ getDefaultRelaxationTime()
| Real SimTK::Force::Thermostat::getDefaultRelaxationTime | ( | ) | const |
Get the initial value for the bath temperature that will be use for the "bath temperature" State variable.
◆ getDefaultNumExcludedDofs()
| 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.
◆ getBoltzmannsConstant()
| Real SimTK::Force::Thermostat::getBoltzmannsConstant | ( | ) | const |
Can't change the value of Boltzmann's constant after construction; this is the value being used.
◆ setNumChains()
| 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.
◆ setBathTemperature()
| 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.
◆ setRelaxationTime()
| 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.
◆ setNumExcludedDofs()
| 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.
◆ getNumChains()
| 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.
◆ getBathTemperature()
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.
◆ getRelaxationTime()
Obtain the current relaxation time.
This is a state variable so can be obtained any time after realization of the Model stage.
◆ getNumExcludedDofs()
| 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.
◆ getNumThermalDofs()
| 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).
◆ getCurrentTemperature()
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.
◆ initializeChainState()
| void SimTK::Force::Thermostat::initializeChainState | ( | State & | ) | const |
This is a solver that initializes thermostat state variables to zero.
◆ setChainState()
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.
◆ getChainState()
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.
◆ calcBathEnergy()
Calculate the total "bath energy" which, when added to the system energy, should yield a conserved quantity (assuming all other forces are conservative).
◆ getExternalPower()
Get the amount of power being applied by the thermostat to the system; sign is positive when energy is coming from the bath.
◆ getExternalWork()
Get the amount of work that has been done by the bath on the system since an arbitrary start time.
◆ setExternalWork()
Set the current value of the work done by the bath to an arbitrary value; normally zero for initialization.
◆ initializeSystemToBathTemperature()
| 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.
◆ setSystemToTemperature()
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:
