In this appendix, the implementation of the AFED method for the isomerization reaction of Sec. 3.2 is discussed. The procedure presented here is generally applicable to any problem for which the reaction coordinate is a distance. Other types of reaction coordinates and generalized procedures are discussed in Ref. [15]. As seen in sect(3.2) the Lagrangian of the system can be rewritten as:
where r, the bond length, is the reaction coordinate, 
is the
unit vector along the bond, 
is the kinetic energy of
the bath, and 
is a mass associated with the r
degree of freedom only. In addition, 
and 
.
From Eq. (41), the equations of motion for r and follow directly:
The equations of motion, Eqs. (44), are slightly more complicated that the usual equations of motion encountered in ordinary molecular dynamics, and, therefore, the problem of constructing a reversible, symplectic integrator for them based on the Liouville operator formalism is discussed below.
The Liouville operator corresponding to Eqs. (44) is
In order to generate a reversible integrator for Eqs. (44), the
following factorization of the classical propagator, 
for a time step, 
is constructed:
Application of this operator to the phase space vector gives the
time evolution of the variables,
r, 
, , 
. The operators that act on the variables
r, , yield a simple time translation,
while the action of 
gives the following evolution of
:
where has the value it has obtained when the operator, 
is applied.
The potentially singular term,
,
can be expanded in a Maclaurin series when is small,
since its limit for 
is finite.
In order to complete the scheme, it is necessary to add
an additional force that constraint to be a unit vector.
This is accomplished by modifying the equation of motion for
according to
In Eq. (48), a constraint force 
has been
added to 
, where 
is a Lagrange multiplier that ensures the condition 
is maintained at all time t. The multiplier is calculated using
the standard procedure. The vector is updated in the absence of
the constraint force using Eq. (47) to produce
the unconstrained evolution 
. The vector 
is then constructed according to
Defining 
, the condition
is then imposed, which leads to the following
expression for 
:
Once the multiplier has been determined, is obtained from
Eq. (49), and the velocity 
at the half
step obtained from
The forces are then recalculated at the new positions, and
is updated with the appropriate updated unconstrained force according
to Eq. (47). Finally, the constraint force is applied, and
the resulting update for , leading to the fully updated
velocity is
where 
is the velocity
just before the constraint force at is applied.
Additionally, it is necessary to apply thermostats to ensure canonical
sampling. As discussed in Sec. 2, a thermostat is applied to the
reaction coordinate r, maintaining its temperature at an
elevated value 
. Because of the constraint on , it is
necessary to thermostat the three components of
together and separately from the remaining bath degrees
of freedom [18]. Apart from this, the application of the
thermostat operators follows the procedure described in detail
in Refs. [13] and [18]. Finally, it should be noted that
the formulation of the integrator via the Liouville operator allows the
AFED method to be easily combined with multiple time scale
integration (RESPA) [14] techniques for the bath.
