As a simple paradigm for solution phase chemical processes, we
consider a simple isomerization reaction of a diatomic molecule
solvated in a Lennard-Jones liquid. The intramolecular potential
of the diatomic is characterized by a single coordinate
, the distance between the two atoms comprising
the diatomic. The intramolecular potential is, again, given
by a quartic double well form:
so that the total potential is given by
In this case, the potential parameters are
10 
K/Å 
, a= 0.22 Å, and 
4.26 Å.
With these parameters, the diatomic can exist in
two stable conformations characterized by bond lengths
of 
, which leads to 
= 4.48 Å
and 
4.04 Å. The two conformations are
separated by a 46.5 kcal/mol
potential barrier at 
.
The diatomic is solvated in a bath of 108 Lennard-Jones particles
with 
90 K, 
3.405 Å and mass m=39.94 amu.
The masses of the two atoms in the diatomic are both the same as the
mass of the Lennard-Jones particles. The
density and temperature of the system are 
0.844
and T= 300 K.
In this example, the reaction coordinate is r, the diatomic bond length. The AFED method requires that the equations of motion be formulated in a coordinate system in which r is an explicit coordinate and is given both a large mass and temperature. To see how this is accomplished, consider the Lagrangian of the system:
where 
is the kinetic energy of the Lennard-Jones
bath particles. First, we transform to center of mass
and relative coordinates for the diatomic, 
and 
,
yielding:
where M=2m and 
are the total and reduced masses
of the diatomic, respectively. Finally, we transform
to spherical polar coordinates 
in the
relative coordinate giving:
where the unit vector 
is the bond orientation
vector. The adiabaticity condition is imposed by
choosing a temperature 
associated with the radial
degree of freedom only and by redefining the Lagrangian
in Eq. (25) according to
where 
is a mass associated with the r
degree of freedom only.
It is important to note that, although it is necessary to work with
the coordinate r explicitly, it is not necessary to
work directly with the angular degrees of freedom 
and
. Rather, we work with the Cartesian components 
the vector 
directly so that the kinetic energy term
involving becomes:
When this is done, the equations of motion generated by
Eq. (26) will be in terms of r and
the three Cartesian components of . In order to ensure that
remains a unit vector as the system evolves in time, it is
necessary to add a simple constraint to the equations of motion
that 
. Details of how the equations
of motion are integrated and how the constraint is implemented
are given in Appendix B.
AFED simulations were carried out using temperatures 
K
and 
K.
At 
, masses of 
10 
amu and
10 amu were used
with a time step of 
5 fs.
At 
, masses of 
amu
and 35953 amu were used.
AFED simulation lengths of various lengths were performed in order
to test convergence. In addition,
bluemoon ensemble simulations were carried out using 11
evenly spaced fixed values of r between 4.0 Å and 4.50 Å.
(This was found to be the minimum number of points needed to
generate an accurate free energy profile.) Each
bluemoon simulation consisted of 
steps with a time step
of 
fs for a total run time of 2.2 ns. Although this is
somewhat longer than is actually needed to generate the free energy
profile via the bluemoon ensemble method, such a long run ensures
a highly converged free energy profile against which the AFED method
can be benchmarked. In Fig. 4, the free energy profiles
obtained from the AFED calculations for the various mass and temperature choices
and that obtained from the bluemoon calculation are shown.
(cf. Eq. (26)) with K (a)
and for K (b).
In both panels, the solid line is the free energy profile obtained
employing the ``bluemoon ensemble'' method [5, 6].
In all cases m=39.94 amu
and the temperature and density are T=300 K and 
,
respectively.
It can be seen that for and 
240 
10 amu,
the agreement between the two methods is good, and for and 
amu, the agreement between the two methods is also good.
In the latter case, a 600 ps run yields the best agreement for 
,
however, even a 200 ps run reproduces the free energy barrier to within
about 5%. For 
, a 600 ps run is insufficiently long
to converge the free energy profile, as can be see in Fig. 4(b), however, if
the trajectory is allowed to run for 2 ns, it is found that the correct
free energy profile is obtained.
In order to compare the efficiency of the AFED and bluemoon methods, the
following analysis was carried out. First, note in Fig. 4, that
the AFED method actually generates a larger range of r
than the bluemoon method. Of course, the bluemoon method can be
made to sample the same range of r by carrying out more simulations.
In this case, 22 evenly spaced points would be required, hence, 22
simulation. Next, an AFED simulation using 22 bins for the histogram
of P(r) was carried out, and the error bar on the
force 
on r in one of the bins
computed and compared to the error bar obtained from the bluemoon ensemble
method for the same value of r. It was found that, in order to obtain an
error bar of size 8 
au, an AFED simulation of 880 ps
was needed compared to a bluemoon calculation of 80 ps. However, since
one AFED simulation generates the full free energy profile, while
22 individual bluemoon simulation are needed, the total simulation
time for the bluemoon method would then be 
ns.
Thus, comparing the total simulation times, shows that the AFED method
is nearly twice as efficient as the bluemoon method. In practice,
one might be content with a restricted range for the reaction coordinate
in the bluemoon method. In the present example, the above analysis shows that,
in this case, the two method are, then, of equal efficiency, however,
the AFED method gives a more complete picture of the free energy profile.
