In order to illustrate the method, we consider a simple model system, which
serves as a paradigm for the general problem of determining free energy profiles.
The model consists
of two degrees of freedom, namely, a reaction coordinate x with
mass
and an additional degree of freedom y with mass
.
The two degrees of freedom are coupled through a potential V(x,y), and
the Hamiltonian of the system is taken to be of the form:
The simple analysis based on Eq. (3) neglects the fact that, in generalized coordinates, the kinetic energy will actually involve a (generally) coordinate-dependent mass metric tensor. A more general analysis and procedure for treating N-particle systems in terms of generalized coordinates will be given in Appendix A.
If
and
are characteristic frequencies of the
x and y motion, then an adiabatic separation between x and
y is achieved by requiring that
.
This can be realized dynamically by choosing the masses such that
. In addition, temperatures
and
associated
with the two degrees of freedom are introduced such that
, as
noted above. In the proceeding discussion, a detailed analysis of the
dynamics and the phase space it generates will be carried out. The goal of
this analysis is to show that, under the usual assumptions of ergodicity, the
resulting phase space distribution function of the reaction coordinate, x,
leads to a free energy profile given by Eq. (1) with
replaced by
.
In order to maintain the two temperatures,
and
in the
system, the dynamics generated by Eq. (3) must
be supplemented by coupling the variables x and y to independent
heat baths or thermostats. Thus, the phase space evolution of
the system is governed by a Liouville operator of the form
where
,
and
and
are
dynamical thermostats, e.g. the Nosé-Hoover chain (NHC) [12] or the
recently introduced generalized Gaussian moment thermostat (GGMT) [13], which
maintain x at temperature
and y at temperature
.
The explicit forms of the thermostat operators are not important for the
present discussion, however, the interested reader is referred to
Refs. [12] and [13] for the detailed expressions of these operators.
For the present discussion, it is sufficient to know that they act on the
momenta
and
in such a way as to control the kinetic
energy fluctuations and generate a Maxwell-Boltzmann distribution
in these variables at the required temperatures and that they involve
a set of additional or ``extended'' phase space variables, denoted
generally by
and
.
It is, again, stressed that two separate
thermostats are needed in order to maintain the two separate
temperatures
and
.
The time evolution of the full phase space
, starting from
an initial condition
, is formally given by
where
is the classical propagator.
Consider the evolution of the system over a time interval
characteristic of the x motion. In order to analyze the dynamics over
such an interval, it is useful to define ``reference system'' Liouville
operators:
and express the total Liouville operator as
where
Using the Trotter theorem, a reversible, symplectic factorization of the classical propagator is constructed according to
In Eq. (9), note that the y-propagator,
,
is left intact. The factorization scheme in Eq. (9) allows the
y-propagator to be evaluated formally exactly, again, using the Trotter theorem.
When Eq. (10) is substituted into Eq. (9), the
result is the formally exact analog of a multiple time scale factorization such
as is discussed in Ref. [14].
The formally exact evaluation of the y-propagator is necessary since, under
the adiabatic conditions of the problem, the time interval,
, is
very long compared to the time scale on which y evolves.
Combining Eqs. (9) and (10)
and acting with the resulting operator on the full phase
space vector
yields the following evolution of
the physical variables [11]:
In Eq. (11),
indicates
the evolution of x under the action of the reference system Liouville operator
up to time t' starting from initial
conditions x',
and
with an analogous meaning for
. Similarly,
indicates the
exact adiabatic, i.e., at fixed x,
evolution of y up to time t' under the action of the
operator in brackets in Eq. (11) starting from
initial conditions y',
,
.
Since
, the motion of the y variable will follow
instantaneously the motion of x dynamically and sample its available
configuration space over the time interval
during which x remains
fixed. In Eq. (11), it can
be seen that the force governing the evolution of x and
is a time average of
over the adiabatic y trajectory.
Assuming that the y motion is ergodic over this time interval, then the time average
of
appearing in Eq. (11) can be replaced
by a configurational average over y according to
where
and
is an effective configurational partition function of x. Note that the
integration over
trivially cancels out of Eq. (12).
Therefore, in the adiabatic limit, an
effective Hamiltonian for the x degree of freedom can be constructed according to
From Eq. (14) follow the statistical mechanical properties of the x degree of freedom. The partition function for x is given by
Therefore, the probability distribution function of x becomes
so that
.
From Eq. (16), it follows that the free energy profile, F(x), which,
by definition is
can be calculated directly from the adiabatically generated probability distribution function, P(x) by
The equivalence between Eqs. (17) and (18)
can be seen by direct substitution of P(x) into Eq. (18).
Thus, the true free energy profile can be obtained, up to an irrelevant constant,
from the probability
distribution function, P(x) generated by the adiabatic dynamics
with
and
via Eq. (18). Note
that Eq. (18) is of the same form as the standard
free energy profile defined in Eq. (1) with T, the
temperature of the ensemble, replaced by
the temperature to which
the x degree of freedom has been heated. Despite the fact that
is Eq. (18), the correct free energy
profile at the temperature,
of the bath/environment is obtained.
Equation (18)
constitutes the central result of this paper.