(360e) Monte Carlo Simulation Using Reversible Mapping between Local Energy Minima

In this paper we demonstrate that a bijective
minimum-to-minimum mapping method [1] dubbed ?MinMap? can be used to construct
new Monte Carlo (MC) moves for molecular simulations.  The method promises to
significantly improve the efficiency of complex moves that allow large jumps in
configuration space.  It is also sufficiently general that it appears to be
comparatively easy to implement for a wide variety of complex molecular systems
and processes.

A MinMap MC move
consists of the following steps:

1. Random selection of a
   small active subset of the molecular system, typically a region of
   space that contains n=5-10 atoms
2. Quenching of the n active atoms by energy minimization while the
   coordinates of the remaining non-active atoms are kept fixed; the atoms
   are additionally constrained to the original active region of space by a suitably
   defined potential
3. Mutation of the active atoms by an appropriately defined transformation
   of their coordinates, composition and/or connectivity, as illustrated by
   the examples below
4. Quenching of the
   mutated atoms
5. The
   mutation/quenching cycle may be broken up into an arbitrary number of
   incremental steps, until a new, energy minimized configuration is obtained
   that completes the desired transformation
6. Reversal of the
   mutation/quenching steps from the new minimum; if the original minimum is
   not recovered then the MC move is not reversible and must be rejected
7. Generation of a set
   of k random excitations of the active atoms about the new
   minimum, where the degree of excitation is governed to be the same as that
   for the original configuration from the original minimum, e.g. as defined
   by the norm of the 3n-dimensional hypersphere of displacement
   vectors [1]
8. Generation of a
   matching set of k excitations from the original minimum, where the
   first of these excitations is the original configuration itself
9. Selection of one of
   the new excitations as the trial configuration
10. Acceptance or
    rejection of the trial configuration using Metropolis sampling [2].

If the selection of these excitations is governed
by their Boltzmann weights, then the sampling corresponds to a generalised form
of the familiar configuration bias method [3].  The new configuration can then
be accepted with probability given by

Acc = Min[1,
r.w_n/w_o]

where w_o and w_n
are the respective Rosenbluth weights [3] for the old and new sets of excitations,
and the ratio r denotes the bias associated with the forward versus
reverse transformations.  If this transformation is symmetric then there is no
bias and r = 1.

The scope of the MinMap approach was
evaluated by implementing a series of new MC moves and testing their simulation
performance for the atomic Lennard-Jones fluid.  Here we consider the following
moves that transform the coordinates of a randomly selected ?cluster? of atoms located
within a sphere of radius a centred at x:

• Inversion
  of the cluster of atoms about x
• Reflection
  of the cluster about a plane passing through x perpendicular to a
  randomly selected coordinate axis
• Rotation
  of the cluster about an axis passing through x and the nearest fixed
  atom

In addition we consider a Tube move,
where the active region is a periodic cylinder of radius a and
effectively infinite length, that is randomly located and oriented parallel to
one of the coordinate axes; the active atoms are displaced an arbitrary
distance (e.g. half the simulation box length) along the cylinder.  Note that
all of the transformations described above are symmetrically reversible.

MC simulations were performed on a pre-equilibrated
Lennard-Jones fluid system with N=151 atoms of diameter σ
and well depth ε.  The canonical or NVT ensemble was
employed for fixed volume V at temperature T=1.5 and number
density ρ=N/V=0.9 in reduced units (σ=
ε =1).  The performance of the new moves was analysed by measuring the
acceptance ratio acc, the average number of displaced atoms n_acc
and the cpu time per move for an SGI Altix 3000 IA64 server.  The atomic mean
squared displacement per cpu second can be used to calculate an effective
diffusivity D which here provides a suitable measure of configurational
sampling rate.  The runs reported here used a single-step mutation and k=20
excitations, which were found to be near optimal in each case.  The notation R(x)
is used to denote the ratio of each quantity x relative to the equivalent
simulation without move relaxation. 

Table 1.  Efficiency of Monte Carlo moves with and without MinMap
relaxation.

The results are summarised in Table 1 for
different move sizes a.  Conservative bounds are shown for cases where
no successful unrelaxed moves could be observed.  It can be seen that the cpu
cost associated with the many minimizations is considerable and increases with
the number of displaced atoms.  However, the enhanced diffusivities show that it
is more than outweighed by the increase in acceptance rates, especially for
larger moves.  Note also that with MinMap the number distributions of active
atoms in accepted moves are similar to those for all attempted moves, whereas
for unrelaxed transformations the distributions differ such that smaller moves
are markedly more likely to be successful, giving R(n_acc)>1.

MinMap-relaxed particle insertions and
deletions were investigated by simulating the Lennard-Jones fluid in the grand
canonical (GC) ensemble.  Here the acceptance ratios for insertions and
removals from an N-particle system are given by r=zV/(N+1)
and r=N/(zV) respectively, where the activity z is calculated
from the desired excess chemical potential [3].  Each simulation consists of a 1:1:1
mix of insertions, deletions and simple single atom displacement moves.  The
inclusion of single atom displacements actually has little effect on the MinMap
simulations, because successful insertions or deletions also induce small local
rearrangements.  Nonetheless they are included here so as to provide optimal circumstances
for the conventional unrelaxed GCMC.  The mutation sequence used for each MinMap
insertion move is to add a repulsive soft sphere of radius 0.6, enlarge it in steps
of 0.1 and then convert it into a full Lennard-Jones atom, where each step in
this sequence is followed by re-minimization [1].  Deletion moves follow the
reverse sequence.  Note that here the stepwise minimization reduces the cpu
time, as well as improving the acceptance ratio by increasing the reversibility
of the mapping. The sampling efficiency is measured by the number of accepted
insertions and deletions per unit cpu, which reflects the rate at which
composition-dependent properties are equilibrated. 

Results for GCMC simulations with a=1.028
and k=20 are given in Table 2, where the quoted average number of displaced
atoms n_acc does not include the
inserted or deleted atom.  We see that the method again leads to significant
increases in the acceptance rate for atom insertions and deletions.  However
this is offset substantially by the greater cpu cost.  The performance of the MinMap
simulations is comparatively moderate because the optimum size for the relaxed
region surrounding single atom insertions and deletions is smaller than for the
moves described in Table 1 (increasing a does improve the acceptance
rate but not the overall sampling efficiency).  As a result, enhanced configuration
sampling performance is only observed at low temperatures and high densities. 
This of course reflects the conditions for which such enhancement is most
desirable.

Table 2.  Efficiency of atom insertions/deletions
in Grand Canonical simulations with and without MinMap relaxation.

In summary, the reversible mapping between
local energy minima can substantially accelerate many Monte
Carlo simulation moves, even for simple atomic systems
such as the Lennard-Jones fluid.  Moreover the method provides the biggest
advantages in the most difficult conditions, such as complex moves involving
many atoms, high densities and/or low temperatures.  It provides a convenient
framework for implementing MC transformations that would otherwise be impossibly
inefficient, which should be of value for simulating more complex molecular
systems.

[1]  D.N. Theodorou.
J.Chem.Phys. 124, 034109, 2006.

[2]  N. Metropolis,
A.W. Rosenbluth, M.N. Rosenbluth, A.N. Teller, and E. Teller. J.Chem.Phys.
21, 125-130, 1987.

[3]  D. Frenkel and
B. Smit, Understanding Molecular Simulation: From Algorithms to Applications,
Academic Press: San
Diego, 1996.