Theoretical analysis of the adiabatic dynamics approach next up previous
Next: Model problems and results Up: Free energy profiles via Previous: Introduction

Theoretical analysis of the adiabatic dynamics approach

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 tex2html_wrap_inline1301 and an additional degree of freedom y with mass tex2html_wrap_inline1369 . 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:

  equation814

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 tex2html_wrap_inline1375 and tex2html_wrap_inline1377 are characteristic frequencies of the x and y motion, then an adiabatic separation between x and y is achieved by requiring that tex2html_wrap_inline1387 . This can be realized dynamically by choosing the masses such that tex2html_wrap_inline1389 . In addition, temperatures tex2html_wrap_inline1391 and tex2html_wrap_inline1393 associated with the two degrees of freedom are introduced such that tex2html_wrap_inline1395 , 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 tex2html_wrap_inline1353 replaced by tex2html_wrap_inline1401 .

In order to maintain the two temperatures, tex2html_wrap_inline1391 and tex2html_wrap_inline1393 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

  equation816

where tex2html_wrap_inline1411 , tex2html_wrap_inline1413 and tex2html_wrap_inline1415 and tex2html_wrap_inline1417 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 tex2html_wrap_inline1391 and y at temperature tex2html_wrap_inline1393 . 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 tex2html_wrap_inline1427 and tex2html_wrap_inline1429 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 tex2html_wrap_inline1431 and tex2html_wrap_inline1433 . It is, again, stressed that two separate thermostats are needed in order to maintain the two separate temperatures tex2html_wrap_inline1391 and tex2html_wrap_inline1393 .

The time evolution of the full phase space tex2html_wrap_inline1439 , starting from an initial condition tex2html_wrap_inline1441 , is formally given by

  equation818

where tex2html_wrap_inline1443 is the classical propagator. Consider the evolution of the system over a time interval tex2html_wrap_inline1445 characteristic of the x motion. In order to analyze the dynamics over such an interval, it is useful to define ``reference system'' Liouville operators:

  eqnarray125

and express the total Liouville operator as

equation820

where

equation822

Using the Trotter theorem, a reversible, symplectic factorization of the classical propagator is constructed according to

  equation149

In Eq. (9), note that the y-propagator, tex2html_wrap_inline1451 , is left intact. The factorization scheme in Eq. (9) allows the y-propagator to be evaluated formally exactly, again, using the Trotter theorem.

  eqnarray161

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, tex2html_wrap_inline1445 , 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 tex2html_wrap_inline1441 yields the following evolution of the physical variables [11]:

  eqnarray184

In Eq. (11), tex2html_wrap_inline1463 indicates the evolution of x under the action of the reference system Liouville operator tex2html_wrap_inline1467 up to time t' starting from initial conditions x', tex2html_wrap_inline1473 and tex2html_wrap_inline1475 with an analogous meaning for tex2html_wrap_inline1477 . Similarly, tex2html_wrap_inline1479 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', tex2html_wrap_inline1489 , tex2html_wrap_inline1491 . Since tex2html_wrap_inline1389 , the motion of the y variable will follow instantaneously the motion of x dynamically and sample its available configuration space over the time interval tex2html_wrap_inline1499 during which x remains fixed. In Eq. (11), it can be seen that the force governing the evolution of x and tex2html_wrap_inline1505 is a time average of tex2html_wrap_inline1507 over the adiabatic y trajectory.

Assuming that the y motion is ergodic over this time interval, then the time average of tex2html_wrap_inline1507 appearing in Eq. (11) can be replaced by a configurational average over y according to

displaymath219

  eqnarray225

where tex2html_wrap_inline1517 and

  equation234

is an effective configurational partition function of x. Note that the integration over tex2html_wrap_inline1429 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

  equation239

From Eq. (14) follow the statistical mechanical properties of the x degree of freedom. The partition function for x is given by

eqnarray246

Therefore, the probability distribution function of x becomes

  equation250

so that tex2html_wrap_inline1531 . From Eq. (16), it follows that the free energy profile, F(x), which, by definition is

  equation257

can be calculated directly from the adiabatically generated probability distribution function, P(x) by

  equation261

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 tex2html_wrap_inline1389 and tex2html_wrap_inline1395 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 tex2html_wrap_inline1391 the temperature to which the x degree of freedom has been heated. Despite the fact that tex2html_wrap_inline1391 is Eq. (18), the correct free energy profile at the temperature, tex2html_wrap_inline1393 of the bath/environment is obtained. Equation (18) constitutes the central result of this paper.


next up previous
Next: Model problems and results Up: Free energy profiles via Previous: Introduction

Mark Tuckerman
Mon Mar 26 04:23:46 EST 2001