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.

 

Move	a	acc %	nacc	D (s-1)	R(cpu)	R(acc)	R(nacc)	R(D)
Inversion	1.028	10.0	3.6	63	42	88	2.5	16
	1.175	3.6	5.7	29	74	2800	1.5	229
	1.322	0.8	7.5	5.1	133	4200	1.1	671
	1.468	0.1	11.1	0.2	252	>105	-	>104
Reflection	1.028	25.5	3.7	62	48	78	1.9	42
	1.175	8.8	5.6	25	79	895	1.3	988
	1.322	1.9	7.8	5.2	131	5000	1.2	57000
	1.468	0.2	10.7	0.5	223	116000	1.1	437000
Rotation	1.028	15.3	3.7	97	43	98	2.2	11
	1.175	5.3	5.5	33	72	1100	1.2	33
	1.322	1.0	7.8	5.8	126	3700	1.1	1200
	1.468	0.1	10.8	0.4	223	2400	1.3	357
Tube	1.028	12.9	3.0	265	33	1000	1.5	202
	1.175	5.5	4.7	90	66	120000	2.4	6300
	1.322	0.8	7.2	8.9	140	>106	-	>105
	1.468	0.1	10.7	0.4	249	>106	-	>106

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.

 

T	ρ	acc (%)	nacc	R(cpu)	R(acc)	R(acc/cpu)
1.5	0.7	1.9	1.9	29	2	0.06
1.5	0.8	1.0	2.6	26	4	0.17
1.5	0.9	0.47	3.2	27	13	0.48
1.5	1.0	0.16	3.8	31	26	0.82
1.2	0.9	0.29	3.2	28	14	0.49
0.9	0.9	0.085	3.4	28	26	0.94
0.6	0.9	0.0009	3.3	24	120	5.0

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.