Thermostats for ``slow'' configurational modes
A.~A.~Samoletov $^{1,2,3}$
C.~P.~Dettmann $^{1}$ and M.~A.~J.~Chaplain $^{2}$
Thermostats are dynamical equations used to model thermodynamic variables such as temperature and pressure in molecular simulations. For computationally intensive problems such as the simulation of biomolecules, we propose to average over fast momentum degrees of freedom and construct thermostat equations in configuration space. The equations of motion are deterministic analogues of the Smoluchowski dynamics in the method of stochastic differential equations.
SmT-epl_arc.tex
\documentclass{article}%
\setlength\textwidth{14.0cm}
\oddsidemargin 30pt
\usepackage{amsmath}%
\usepackage{amsfonts}%
\usepackage{amssymb}%
\usepackage{bm}%
\usepackage{graphicx}

\begin{document}

\title{Thermostats for ``slow'' configurational modes}

\author{A.~A.~Samoletov $^{1,2,3}$,
C.~P.~Dettmann $^{1}$ and M.~A.~J.~Chaplain $^{2}$ }

\date{ }

\maketitle

${}^1${\it Department of Mathematics, University of Bristol -
Bristol
BS8 1TW, UK }%

${}^2${\it Division of Mathematics, University of Dundee - Dundee
DD1 4HN, UK}%

${}^3${\it Institute for Physics and Technology, NASU - Donetsk
83114, Ukraine}


\begin{abstract}
Thermostats are dynamical equations used to model thermodynamic
variables such as temperature and pressure in molecular simulations.
For computationally intensive problems such as the simulation of
biomolecules, we propose to average over fast momentum degrees of
freedom and construct thermostat equations in configuration space.
The equations of motion are \emph{deterministic} analogues of the
Smoluchowski dynamics in the method of stochastic differential
equations.
\end{abstract}

\section{Introduction}

One of the most ambitious challenges in mathematical modelling of
biological processes is to describe dynamics of two major biological
events within a cell - DNA molecule replication and transcription. In
both cases, dynamical properties, especially large amplitude
conformational changes of the double-stranded DNA molecule play a
vital part. The principal physical feature of biological
functioning of biomolecules is that they operate at ambient
physiological temperature, pressure and solvent conditions. Thus,
the surrounding physiological solvent plays the role of
a thermostat, among others. To properly thermostat this dynamics,
especially on a biological (``slow'') time scale, we propose a
suitable and effective temperature control for both the practice of
numerical simulation and the general theory of dynamical systems.

For a recent and comprehensive review of the problem outlined above
we refer to \cite{mp04}. Current approaches commonly use the
Nos\'{e}-Hoover thermostat method (for review see
\cite{T1,T2,T3,T4}). This canonical thermostat method involves
integration of both position and momentum phase space variables.
However, for problems that are related to slow conformational
changes of biomolecules, integration of fast momentum variables
appears superfluous from a theoretical point of view (unobservable
variables) as well expensive in the sense of numerical simulation.

In this Letter, a novel approach to the problem of slow
conformational changes in thermostatting dynamics is presented. The
method is based on an analogy with the derivation of the Smoluchowski
dynamics~(Eq.~(\ref{Smoluchowski}) below) from the Langevin
stochastic dynamics~(Eq.~(\ref{Langevin}) below)
\cite{Gard04,kramers40,as99}. The configurational deterministic
thermostat is constructed so as to effect the temperature control
via certain dynamics of the relaxation rate variable and also
involves a dynamically fluctuated collective force to ensure the
ergodicity property. As the result, the temperature control is
connected to the specific configurational temperature recently
introduced in \cite{rugh97,evans00} in a different context. The new
configurational thermostat can be combined with complementary
temperature control via a dynamically fluctuated virial function
\cite{hamilton90} that helps to enhance the efficiency of the
thermostat and, more importantly because biological molecules are
functioning at a constant pressure, implements the pressure control
into the corresponding dynamics. Moreover, the configurational
thermostat admits stimulation by a chain of thermostats similar to
Nos\'{e}-Hoover \cite{chain} and a stochastically driven method (see
below). To test the new configurational thermostat, corresponding
simulations of a one-dimensional harmonic and the Morse oscillator
dynamics is given, providing a stringent test of the ergodicity
property.

Remark: Braga and Travis \cite{BT05} have recently proposed
a thermostat based on the Smoluchowski equation and configurational
temperature ideas (but retaining the momentum variables), citing an
early preprint of this paper \cite{SCD04}.
\section{Preliminaries}

The early and successful attempts to describe the dynamics of a
mechanical system being in contact with an environment playing a
role of the thermostat is based on the concept of the stochastic
differential equation \cite{Gard04}. The case of one-dimensional
motion of a particle mass $m$ in a potential field $V(q)$ provides the
characteristic example,
\begin{equation}
\label{Langevin}m\dot{q}=p,\quad\dot{p}=-\nabla V(q)-\gamma
p+\sqrt{2D}f(t),
\end{equation}
where the friction coefficient $\gamma$ and the intensity $D$ of the
external random force $f(t)$ are connected by the relation, $D=m\gamma
k_{\mathrm{B}}T$\:; here $f(t)$ is the generalized Gaussian
stochastic process, \textquotedblleft white noise\textquotedblright,
with characteristic cumulants $\left\langle f(t)\right\rangle =0,\:
\left\langle f(t)f(t^{\prime})\right\rangle =\delta(t-t^{\prime})$.
The equilibrium solution of the corresponding Fokker-Planck equation
\cite{Gard04} is known, $\rho_{\infty}\propto \exp\left[  -\left(
p^{2}/2m+V(q)\right)  /(k_{\mathrm{B}}T) \right]  $.

Langevin's equation (\ref{Langevin}) is a prototype of the
Nos\'{e}-Hoover deterministic dynamics (after generalisation of the
last with the term $\bm{\xi}$ as in Eq.~(\ref{VSdyn}) it is
especially evident). But while the Langevin dynamics generate all
the sample trajectories and the corresponding measure, the
Nos\'{e}-Hoover dynamics produce a single sample trajectory with the
correct canonical ensemble statistics when ergodicity holds.

Often details of the dynamics of a system on short time scales are not
needed for a dynamical description of the observable
variables. However, then the Smoluchowski limit of (\ref{Langevin}),
\begin{equation}
\label{Smoluchowski}m\dot{q}=-\tau\nabla
V(q)+\sqrt{2mk_{\mathrm{B}}T\tau }f(t),
\end{equation}
where $\tau=\gamma^{-1}$, is an appropriate formulation
\cite{Gard04,kramers40,as99}. The position variable only is involved
in this equation. Formally, it is supposed that the momentum
variable relaxes to the local equilibrium state. The corresponding
Smoluchowski equation has the Boltzmann distribution as the
equilibrium solution \cite{Gard04}. Eq.~(\ref{Smoluchowski}),
without the random perturbation, appears as a dissipative dynamics
with $V(q)$ playing a role of the Lyapunov function,
$\dot{V}=-\tau\left( \nabla V(q)\right) ^{2}\leq0$. Thus, the full
dynamics is a superposition of relaxation to a minimum of the
potential and random perturbations that occasionally expel the
system outside the vicinity of the minimum. This process
equilibrates the system.

\section{Configurational thermostat}

It is reasonable, in the spirit of deterministic thermostat methods,
to conjecture that it is possible to use the relaxation time $\tau$
for thermostatting configurational degrees of freedom when momentum
variables are still relaxed in their local equilibrium state. Of
course, in this case the sign of $\tau$ is not fixed and $V$ loses
its meaning as a Lyapunov function. In a sense, it means that time
can go back as well as forward.

First, consider the simple dynamical equations for $\mathcal{N}$
degrees of freedom,
\begin{equation}
\label{iso-S}m\mathbf{\dot{q}}=-\tau\bm{\nabla} V(q),
\end{equation}
where $\tau$ is a constant, but by analogy with Nos\'e-Hoover will
be endowed with its own equation of motion below. By the change of
variables, $\mathbf{x}=\sqrt{m}\mathbf{q}$, it is possible to
exclude all masses from the formulae in what follows, but we prefer
to save the physical notation. Requiring that the corresponding
Liouville equation has the Boltzmann distribution,
$\rho_{\infty}\propto\exp\left[  -V(q)/(k_{\mathrm{B}}T)\right] $,
as a steady state solution, we arrive at the condition that involves
the temperature in the dynamics,
\begin{equation}
\label{T}\sum\frac{1}{m}\left[  \Delta
V(q)-\frac{1}{k_{\mathrm{B}}T}\left( \bm{\nabla} V(q)\right)
^{2}\right]  =0,
\end{equation}
where $\sum$ denotes summation over all particles of the system.
After time averaging, Eq.~(\ref{T}) appears as the definition of the
recently introduced so-called configurational temperature
\cite{rugh97,evans00}. Currently it is used in molecular dynamics
simulations \cite{delh04}. In a more general context, in the case of
a presupposed anisotropy in the system, let us assume that $\tau$ in
(\ref{iso-S}) is a matrix, $\tau\to\bm{\Gamma}$. Then the dynamics
take the form, $m\mathbf{\dot{q}}=-\bm{\Gamma}\bm{\nabla} V(q)$, and
the condition that involves the temperature in the dynamics is
\[
\sum\frac{1}{m}\left[  (\bm{\nabla},\bm{\Gamma} \bm{\nabla})
V(q)-\frac {1}{k_{\mathrm{B}}T}\left(  \bm{\nabla} V(q),\bm{\Gamma}
\bm{\nabla} V(q) \right)  \right]  =0.
\]
Note that this condition involves the presupposed time scales in the
dynamical temperature control. Conventional Nos\'e-Hoover methods do
not allow such a generalisation. In what follows we also consider
$\tau$ as a scalar. On the other hand, it is useful to keep in mind
the possibility of the generalisation.

We now attempt to generate statistics as in the Nos\'e-Hoover scheme
by making $\tau$ an independent variable in~(\ref{iso-S}). It is
easily seen that this is too simple. At an equilibrium point
$\bm{\nabla} V=0$, that is, all forces are zero, the evolution
comes to a halt and no longer fluctuates, irrespective of the time
dependence of $\tau$. For initial conditions with nonzero
forces~(\ref{iso-S}) after a (positive or negative) change of time
variable, it is a gradient flow as defined in~\cite{KH}, and it is
easy to show that all trajectories move along paths in $\mathbf{q}$
with equilibrium points at either end. In short, the system is not
ergodic. Note also that~(\ref{T}) is singular when $\bm{\nabla}
V=0$.

The way to overcome the difficulty is suggested by the Smoluchowski
stochastic equation. In this equation the ergodic motion is ensured
by the random forcing. Hence, we need to add a deterministic analogue
of the random force term in~(\ref{Smoluchowski}). Let us consider,
instead of (\ref{iso-S}), the dynamical equations,
\begin{equation}\label{xi}
m\mathbf{\dot{q}}=-\tau\bm{\nabla}V(q)+\bm{\xi},
\end{equation}
where $\bm{\xi}$ are constant vectors, but they will be endowed with
their own equation of motion below. To detail the nature of vectors
$\bm{\xi}$, three principal cases are possible: (a) All
$\bm{\xi}=\{\bm{\xi}_i\}_{i=1}^N$, where $N$ is number of particle
in the system, can be varied independently; (b) All
$\bm{\xi}_i=\bm{\xi}$ are varied identically; (c) There exist a
preferred direction, $\mathbf{e}$, where $\mathbf{e}$ is a constant
unit vector of physical (\textit{i.e.} three dimensional) space, and
only one variable, $\xi$, $\bm{\xi}=\xi \mathbf{e}$, is
varied. Requiring that the Liouville equation corresponding to
(\ref{xi}) has the Boltzmann distribution as a steady states solution
and varying $\bm{\xi}$ according to cases (a), (b) and (c), together
with temperature control condition (\ref{T}) we correspondingly
obtain the following conditions,
\begin{equation}\label{fxi}
(a)\quad \bm{\nabla}V(q)=0, \qquad (b)\quad \sum
\frac{1}{m}\bm{\nabla}V(q)=0, \qquad (c)\quad \sum
\frac{1}{m}\mathbf{e}\cdot\bm{\nabla}V(q)=0.
\end{equation}
These conditions do not involve temperature but the thermalized
forces acting in the system. All of them, as well as atheir
combination, are candidates for simulating the deterministic analogue
of the random force, chosen according to the problem under
consideration. For example, in respect of the Peyrard-Bishop
dynamical model of the DNA molecule \cite{mp04}, case (c) appears to
be appropriate. The physical sense of the conditions above are the
following. Case (a): the force acting on a particle in the system
equals zero (static equilibrium of forces). Case (b): the static
equilibrium of forces is not required but the total force acting on
the system equals zero (stability of the system). Case (c): stability
of the system in the direction $\mathbf{e}$.

It is practical to remark that a thermostatting dynamics more
general than (\ref{xi}) is possible,
\begin{equation}
\label{VSdyn}m\mathbf{\dot{q}}=-\tau\nabla V(q)+\eta
m\mathbf{q}+\bm{\xi},
\end{equation}
where the term $\eta m \mathbf{q}$ is suggested by the virial
thermostatting scheme that first appears in~\cite{hamilton90} in the
context of the harmonic oscillator. Requiring that the corresponding
(\ref{VSdyn}) Liouville equation has the Boltzmann distribution
as a steady state solution we obtain together with (\ref{T}) and
\ref{fxi}) the following condition on the virial function,
\begin{equation}\label{V}
    \mathcal{N}k_{\mathrm{B}}T-\sum\mathbf{q}\cdot\nabla V=0.
\end{equation}
In this case a double temperature control is provided. The $\eta$
term is not a mandatory temperature control for our configurational
thermostat. We could consider only the $\tau$ term by a trivial
modification of~(\ref{VSdyn}) and subsequent equations but consider
a more general equations of the form (\ref{VSdyn}). The reason is that
the virial function involves the pressure in the dynamics and this is
important for biologically oriented models. For the sake of
definiteness and keeping in mind a future application to the
Peyrard-Bishop model of the DNA dynamics we here fix case (c). When
$\xi$ is endowed with its own equation of motion, it provides a
shaking of the system along the direction $\mathbf{e}$.

Now we have three parameters which we group together as a 3-vector
$\bm{\alpha}=(\tau,\eta,\xi)^{\mathrm{T}}$. As in~(\ref{T}),
(\ref{fxi}) and (\ref{V}) we find a stationary solution of the
Liouville equation of the form
\[
\rho_{\infty}\propto\exp\left[  -(V(q)+\bm{\alpha}^{\mathrm{T}%
}\bm{Q}\bm{\alpha}/2)/(k_{\mathrm{B}}T)\right]
\]
where $\bm{Q}$ is a positive definite real symmetric matrix. It is,
of course, not possible to justify consistently that parameters
$\bm{\alpha}$ %are statistically independent from the system variables and
have a Gaussian distribution at equilibrium. It is just a reasonable
assumption that leads to the simplest form of self-consistent
thermostat dynamics. Intuitively, we can justify the Gaussian
character of the parameters at equilibrium by the limiting theorems
of probability theory. Note that when multiple thermostats are
present, they are usually assumed to be uncoupled (diagonal
$\bm{Q}$); mathematically it is not required and we do not feel that
this is physically necessary, thus we consider that the more general
case of coupled thermostats may be useful. For comparison with the
uncoupled case we define $\tilde{\bm{Q}}=diag(Q_{\tau
\tau},Q_{\eta\eta},Q_{\xi\xi})$ to be the matrix with only the
diagonal components of $\bm{Q}$. Eq.~(\ref{T}) becomes
\begin{equation}
\label{VStemp}%
\tilde{\bm{Q}}\bm{g} \equiv\left(
\begin{array}
[c]{c}%
\sum\frac{1}{m}\left[  (\nabla V)^{2}-k_{\mathrm{B}}T\Delta V\right] \\
\mathcal{N}k_{\mathrm{B}}T-\sum\mathbf{q}\cdot\nabla V\\
-\sum\frac{1}{m}\mathbf{e}\cdot\nabla V
\end{array}
\right)  =0,
\end{equation}
which now defines $\bm{g}$.

Let us realize the main conjecture the configurational thermostat
scheme and allow the components of $\bm{\alpha}$ to fluctuate so
that~(\ref{VStemp}) holds only after time averaging; we write
\begin{equation}
\label{alphadot}\dot{\bm{\alpha}}=\bm{G}
\end{equation}
where $\bm{G}$ is as yet an undetermined vector of functions. Now
requiring the same condition for the solution of the Liouville
equation we find that the only solution is
\begin{equation}
\label{G}\bm{G}=\bm{Q}^{-1}\tilde{\bm{Q}}\bm{g}.
\end{equation}
Thus the only undetermined parameters of our thermostatting scheme
are the components of the positive definite real symmetric matrix
$\bm{Q}$, and in the uncoupled case $\tilde{\bm{Q}}=\bm{Q}$ we have
$\bm{G}=\bm{g}$.

We can now ask whether the addition of new variables $\eta$ and/or
$\xi$ will remove the lack of ergodicity implied by the potential
flow argument applying to~(\ref{iso-S}). A partial answer is
provided by the Frobenius theorem of differential
geometry~\cite{Lang}, which in our case states that an integral
surface exists (hence the dynamics are definitely not ergodic) if a
vector space containing the terms in the equation for $\mathbf{q}$
but smaller than the full phase space is closed under Lie brackets.
For realistic potentials (not the harmonic oscillator) this is very
unlikely since multiple derivatives of $V$ are almost always
linearly independent. If the theorem does not apply we are in the
same situation as for nonthermostatted nonintegrable many particle
systems, which are often assumed to be ergodic, at least for
practical purposes.

Since under the transformation $t\rightarrow-t,
\bm{\alpha}\rightarrow -\bm{\alpha}$ the equations of motion
(\ref{VSdyn}) and (\ref{alphadot}) are still unchanged, they are
time reversible.

To find a mechanically important integral of motion of system
(\ref{VSdyn}), we need to add a redundant variable. Indeed, consider
the balance of the mechanical work along trajectories of
Eqs.~(\ref{VSdyn}) and  (\ref{alphadot}),%
\[
\sum-\nabla V\cdot d \mathbf{q}= d (\bm{\alpha}^{\mathrm{T}%
}\bm{Q}\bm{\alpha}/2)+k_{\mathrm{B}}T\left(  \sum\frac{\Delta V}{m}%
\tau-\mathcal{N}\eta\right)   d t.
\]
To obtain an exact differential equation, it is necessary to set%
\[
\left(  \sum\frac{\Delta V}{m}\tau-\mathcal{N}\eta\right)   d t= d
\theta.
\]
In that case, the following first integral takes place,%
\[
I_{\mathrm{S}}=V\left(  q\right)  +\bm{\alpha}^{\mathrm{T}}%
\bm{Q}\bm{\alpha}/2+k_{\mathrm{B}}T\theta.
\]
Since the origin of the redundant variable $\theta$ is arbitrary, it
is always possible for an arbitrary fixed trajectory to set
$I_{\mathrm{S}} =0$. This integral of motion is useful as a control
parameter in numerical simulations. Besides, it clearly relates to
the equilibrium distribution $\rho_{\infty}$ and thus can be
considered as a first step to reformulation of Eqs.~(\ref{VSdyn}%
)-(\ref{alphadot}) in terms of a free energy functional. Recall that
the corresponding Nos\'{e}-Hoover integral of motion~\cite{T1,T2,T3,T4} is
given by its Hamiltonian - no Hamiltonian is possible here since the
momentum does not appear explicitly.

\section{Stimulated configurational thermostats}\label{stimul}

Since the configurational dynamics above result in the Gaussian
equilibrium fluctuation of thermostat variables, the latter admits
reinforcing by a chain of equations analogous to the Nose-Hoover
chain thermostat \cite{chain}. The chain method consists in
including a subsidiary sequence of dynamical variables,
$\{\bm{\alpha}_{i}\}$, into a thermostat scheme such that
asymptotically, in the equilibrium distribution, they are
independent Gaussian variables,
\[
\rho_{\infty}\propto\exp[-(V(q)+\frac{1}{2}\bm{\alpha}^{\mathrm{T}%
}\bm{Q}\bm{\alpha}+\sum\nolimits_{(i)}\frac{1}{2}%
\bm{\alpha_i}^{\mathrm{T}}\bm{Q}_{i}\bm{\alpha_i})/(k_{\mathrm{B}}T)].
\]
The corresponding dynamics are not unique. We have obtained a clear
method for generalizing the chain, but since this does not directly
relate to our main topic we do not discuss the details. Instead we
cite the example of the chain rule that has been used in our test
simulation,
\begin{align}
\label{chain}%
m{\mathbf{\dot{q}}}=-\tau\nabla V(q)+\xi\boldsymbol{e},\;
\dot{\tau}=g_{\tau}+\tau_{1}\tau, \; \dot{\tau}_{i}=\frac{1}{Q_{\tau
_{i}}}\left(
k_{\text{\textrm{B}}}T-Q_{\tau_{i-1}}\tau_{i-1}^{2}\right)
+\tau_{i+1}\tau_{i}, \; \dot{\xi} =g_{\xi},
%m{\mathbf{\dot{q}}} =-\tau\nabla
%V(q)+\xi\boldsymbol{e},\quad \dot{\tau}
%& =g_{\tau}+\tau_{1}\tau,\quad \dot{\xi} =g_{\xi}, \quad\\
%\dot{\tau}_{i} &=\frac{1}{Q_{\tau _{i}}}\left(
%k_{\text{\textrm{B}}}T-Q_{\tau_{i-1}}\tau_{i-1}^{2}\right)
%+\tau_{i+1}\tau_{i},\nonumber
\end{align}
where $i=1,\ldots,\mathrm{M},\: \tau_{0}=\tau,\:
\tau_{\mathrm{M}+1}\equiv0$. It is a simple chain of total length
$\mathrm{M}$.

It should also be noted that in spite of its popularity, the
effectiveness of the chain method for computing non-equilibrium
properties has been questioned \cite{branka00}.

It is possible to stimulate the Gaussian fluctuation of the
thermostat variables by the process of Brownian motion. This scheme
has the advantage of ensuring the ergodicity property. The
stimulation, similar to the chain one, may be done in a few ways,
applying to one or more of the variables $\tau$, $\eta$ and $\xi$. In
general we have%
\begin{align}
\label{stst}m\mathbf{\dot{q}} =-\tau\nabla V(q)+\eta
m\mathbf{q}+\xi\mathbf{e}
,\quad\dot{\bm{\alpha}}=\bm{G}-\bm{\Lambda}\bm{\alpha}+\sqrt
{2\bm{D}}\bm{f}(t) ,
\end{align}
where now $\bm{\Lambda}$ and $\sqrt{2\bm{D}}$ are positive definite
real symmetric matrices, and $\bm{f}(t)$ is a vector of independent
white noise components. The Liouville equation corresponding to
(\ref{stst}), averaged over all the realizations of $\bm{f}(t)$, has
the form of the Fokker-Planck equation and the Boltzmann
distribution as a steady state solution only if
\[
{k_{\mathrm{B}}T}\,\bm{\Lambda}=\bm{D}\bm{Q}.
\]
We consider that the most physical case is when the noise is used
only for temperature control, that is, for $\tau$ and $\eta$ only.
%The chain scheme can be done analogously.
We do not establish
extreme generality here because our main aim is the presentation of
the idea of the \emph{deterministic} configurational thermostat. The
effectiveness of Eqs.~(\ref{stst}) for systems far from equilibrium
is not clear.
\begin{figure}[ptb]
\includegraphics[clip,width=120mm]{f1.eps}\caption{Probability distributions
of position variable (shown on background of exact analytical
distribution) of the harmonic (1) and Morse (2) oscillator.
Probability densities are calculated as normalized sojourn
distributions. Correspondingly, thermostat is: (a1)-(a2)
non-stimulated and uncorrelated (Eq.~(\ref{VSdyn}), $\dot
\tau=g_{\tau}, \eta\equiv0, \dot\xi=g_{\xi}$); (b1)-(b2)
non-stimulated and uncorrelated but under double temperature control
(Eq.~(\ref{VSdyn}), $\dot\tau=g_{\tau}, \dot\eta=g_{\eta},
\dot\xi=g_{\xi}$); (c1)-(c2) uncorrelated but stimulated by the
chain rule (Eq.~(\ref{chain}), $\mathrm{M}=1$); (d1)-(d2) stimulated
by stochastic process (Eq.~(\ref{stst}), $\eta\equiv0, D=1$). All
simulations are performed at $Q_{\tau}=Q_{\xi}=1$, $Q_{\eta}=0.1$
for $t=10^{4}$ (squares) and $t=10^{6}$ (black circles). If squares
do not shown then distributions are sufficiently near or coincide
(a1).}%
\label{fig1}%
\end{figure}

\section{Test numerical simulations}

The harmonic oscillator is both a simple and an important physical
system. At the same time, it reveals the ergodicity problem in the
canonical ensemble simulation. For this reason, it is important to
test the Smoluchowski thermostat method capable of generating
the Boltzmann distribution for a single harmonic oscillator in one
dimension, $V(q)=q^{2}/2$. Then it is reasonable to simulate
another one-dimensional system, `good' from the point of view of the
Frobenius theorem, and to compare results. We choose for this purpose
the Morse oscillator,
$V(q)=V_{0}(1-\exp(-aq))^{2}+\mathrm{k}q^{2}/2$. Simulations are
performed using global parameters $m=1$ and $k_{\mathrm{B}}T=1$, and
the Morse potential parameters $V_{0}=0.25$, $a=2$,
$\mathrm{k}=0.25$. Figure~\ref{fig1} shows the probability
distribution of the position variable calculated with the four
simplest configurational thermostats. Note the effectiveness of the
double temperature control. The ability of the Smoluchowski thermostat
to reproduce the correct distribution function, $\rho(q)$, even with
absolute minimum of this thermostat capacity, demonstrates its great
potential for application.

\section{Conclusion}

An innovative constant temperature thermostat, the configurational
thermostat, exclusively involving dynamics of the configurational
variables has been introduced. It poses the general problem of the
derivation of a thermostatting dynamics for slow dynamical variables
and outlines at least one way of the solution. For practical
purposes, the new thermostatting scheme can easily be combined with
complementary temperature control via dynamical fluctuation of the
virial function. This combination helps to enhance the efficiency of
the thermostat temperature control as well as to implement pressure
into the dynamics. It is also relevant to emphasise the appearance of
the dynamically fluctuated forces and correlated temperature control
in the presented thermostatting scheme. We finally remark that
the proposed method is applicable to thermostatting such
macro/meso-scale models as the Ginzburg-Landau or reaction-diffusion
dynamical equations.

 AS is grateful for support from the University of Bristol and the
University of Dundee, from CCP5, and from the Royal Society
(London).

\begin{thebibliography}{99}                                                                                               %


\bibitem{mp04}{Peyrard M.},  {\it Nonlinearity}, {17} (2004) R1.

\bibitem{T1} {Morriss G.P. \and Dettmann C.P.},
{\it Chaos}, {8}(1998) {321}.

\bibitem{T2} {Tuckerman M.E. \and
Martyna~G.J.}, {\it J. Phys. Chem. B}, {104}(2000) {159}.

\bibitem{T3}  {Kusnezov D., Bulgac A. \and Bauer W.},
{\it Ann. Phys.}, {204} (1990) {155}.

\bibitem{T4} {Hoover Wm.G., Aoki K., Hoover C.G. \and De Groot
S.V.},
{\it Physica D}, {187} (2004) {253}.

\bibitem{Gard04}  {Gardiner C.W.},
{\it Handbook of Stochastic Methods, 3rd~ed.}  (Springer-Verlag,
Berlin) 2004.

\bibitem{kramers40} {Kramers H.A.},  {\it Physica}, {7} (1940) {284}.

\bibitem{as99}  {Samoletov A.A.},  {\it J. Stat. Phys.}, {96} (1999) {1351}.

\bibitem{rugh97}  {Rugh H.H.}, {\it Phys. Rev. Lett.} {78} (1997) {772}.

\bibitem{evans00}  {Jepps O.G., Ayton G. \and Evans D.J.},
{\it Phys. Rev. E}, {62} (2000) {4757}.

\bibitem{hamilton90}  {Hamilton I.P.}, {\it Phys.Rev. A}{42}, {1990}{7467};
 {L'Heureux I. \and Hamilton I.} {\it Phys. Rev.
E}{47} (1993) {1411}.

\bibitem{chain}  {Martyna G.J., Klein M.L. \and Tuckerman M.},
{\it J. Chem. Phys.}{97} (1992) {2635}.

\bibitem{BT05}  {Braga C. \and Travis K.P.},  {\it J.
Chem. Phys.}{123} (2005) {134101}.

\bibitem{SCD04}  {Samoletov A., Chaplain M. \and Dettmann C.},
 {\it
arXiv/physics/0412163}, {2004}.

\bibitem{delh04}  {Delhommelle J.},
{\it Phys. Rev. B}, {69} (2004) {144117};  {Delhommelle J., Petravic
J. \and Evans D.J.},  {\it J. Chem. Phys.}{120} (2004) {6117}.

\bibitem{KH}  {Katok A. \and Hasselblatt B.}
{\it Modern Theory of Dynamical Systems}  (Cambridge University
Press, Cambridge)  (1995).

\bibitem{Lang}  {Lang S.}
{\it Introduction to Differentiable Manifolds}  (Springer-Verlag,
Berlin)  (2002).

\bibitem{branka00}  {Branka A.C.}, {\it Phys. Rev. E} {61}(2000){4769}.

\end{thebibliography}

\end{document}

physics0412163.extracted_bib
\begin{thebibliography}{99}                                                                                               %


\bibitem{mp04}{Peyrard M.},  {\it Nonlinearity}, {17} (2004) R1.

\bibitem{T1} {Morriss G.P. \and Dettmann C.P.},
{\it Chaos}, {8}(1998) {321}.

\bibitem{T2} {Tuckerman M.E. \and
Martyna~G.J.}, {\it J. Phys. Chem. B}, {104}(2000) {159}.

\bibitem{T3}  {Kusnezov D., Bulgac A. \and Bauer W.},
{\it Ann. Phys.}, {204} (1990) {155}.

\bibitem{T4} {Hoover Wm.G., Aoki K., Hoover C.G. \and De Groot
S.V.},
{\it Physica D}, {187} (2004) {253}.

\bibitem{Gard04}  {Gardiner C.W.},
{\it Handbook of Stochastic Methods, 3rd~ed.}  (Springer-Verlag,
Berlin) 2004.

\bibitem{kramers40} {Kramers H.A.},  {\it Physica}, {7} (1940) {284}.

\bibitem{as99}  {Samoletov A.A.},  {\it J. Stat. Phys.}, {96} (1999) {1351}.

\bibitem{rugh97}  {Rugh H.H.}, {\it Phys. Rev. Lett.} {78} (1997) {772}.

\bibitem{evans00}  {Jepps O.G., Ayton G. \and Evans D.J.},
{\it Phys. Rev. E}, {62} (2000) {4757}.

\bibitem{hamilton90}  {Hamilton I.P.}, {\it Phys.Rev. A}{42}, {1990}{7467};
 {L'Heureux I. \and Hamilton I.} {\it Phys. Rev.
E}{47} (1993) {1411}.

\bibitem{chain}  {Martyna G.J., Klein M.L. \and Tuckerman M.},
{\it J. Chem. Phys.}{97} (1992) {2635}.

\bibitem{BT05}  {Braga C. \and Travis K.P.},  {\it J.
Chem. Phys.}{123} (2005) {134101}.

\bibitem{SCD04}  {Samoletov A., Chaplain M. \and Dettmann C.},
 {\it
arXiv/physics/0412163}, {2004}.

\bibitem{delh04}  {Delhommelle J.},
{\it Phys. Rev. B}, {69} (2004) {144117};  {Delhommelle J., Petravic
J. \and Evans D.J.},  {\it J. Chem. Phys.}{120} (2004) {6117}.

\bibitem{KH}  {Katok A. \and Hasselblatt B.}
{\it Modern Theory of Dynamical Systems}  (Cambridge University
Press, Cambridge)  (1995).

\bibitem{Lang}  {Lang S.}
{\it Introduction to Differentiable Manifolds}  (Springer-Verlag,
Berlin)  (2002).

\bibitem{branka00}  {Branka A.C.}, {\it Phys. Rev. E} {61}(2000){4769}.

\end{thebibliography}