|
Simbody
|
|
Simbody
|
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. | |
| Thermostat & | excludeMobilizedBody (MobilizedBodyIndex) |
| TODO: not implemented yet. | |
| Thermostat & | setDefaultNumChains (int numChains) |
| Set the default (state independent) number of Nose'-Hoover chains. | |
| Thermostat & | setDefaultBathTemperature (Real bathTemperature) |
| Set the default (state independent) bath temperature. | |
| Thermostat & | setDefaultRelaxationTime (Real relaxationTime) |
| Set the default (state independent) relaxation time. | |
| 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 | getDefaultNumChains () const |
| Get the initial value for the number of chains that will be used for the "number of chains" State variable. | |
| Real | getDefaultBathTemperature () const |
| Get the initial value for the bath temperature that will be use for the "bath temperature" State variable. | |
| Real | getDefaultRelaxationTime () const |
| Get the initial value for the bath temperature that will be use for the "bath temperature" State variable. | |
| 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. | |
| Real | getBoltzmannsConstant () const |
| Can't change the value of Boltzmann's constant after construction; this is the value being used. | |
| const Thermostat & | setNumChains (State &, int numChains) const |
| Set the actual number of Nose'-Hoover chains to be used. | |
| const Thermostat & | setBathTemperature (State &, Real Tb) const |
| Set the bath temperature which serves as the target temperature for the thermostat. | |
| 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. | |
| 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. | |
| int | getNumChains (const State &) const |
| Obtain the current number of Nose'-Hoover chains in use. | |
| 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. | |
| Real | getRelaxationTime (const State &) const |
| Obtain the current relaxation time. | |
| 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. | |
| int | getNumThermalDofs (const State &) const |
| Return the number of thermal degrees of freedom being used in the definition of temperature for this thermostat. | |
| 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. | |
| void | initializeChainState (State &) const |
| This is a solver that initializes thermostat state variables to zero. | |
| void | setChainState (State &, const Vector &) const |
| Set the thermostat state variables to particular values. | |
| Vector | getChainState (const State &) const |
| Return the current values of the thermostat chain variables. | |
| 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). | |
| 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. | |
| 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. | |
| 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. | |
| void | initializeSystemToBathTemperature (State &) const |
| Set the controlled system to a set of randomized velocities which yields the bath temperature. | |
| void | setSystemToTemperature (State &, Real T) const |
| Set the controlled system to a set of randomized velocities which yields a particular temperature. | |
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.
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.
[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).
| 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& 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 |
| 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 (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.
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.
| 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.
1.7.3