^{1,a)}and Phillip L. Geissler

^{1}

### Abstract

We introduce a “virtual-move” Monte Carlo algorithm for systems of pairwise-interacting particles. This algorithm facilitates the simulation of particles possessing attractions of short range and arbitrary strength and geometry, an important realization being self-assembling particles endowed with strong, short-ranged, and angularly specific (“patchy”) attractions. Standard Monte Carlo techniques employ sequential updates of particles and can suffer from low acceptance rates when attractions are strong. In this event, collective motion can be strongly suppressed. Our algorithm avoids this problem by proposing simultaneous moves of collections (clusters) of particles according to gradients of interaction energies. One particle first executes a “virtual” trial move. We determine which of its neighbors move in a similar fashion by calculating individual bond energies before and after the proposed move. We iterate this procedure and update simultaneously the positions of all affected particles. Particles move according to an approximation of realistic dynamics without requiring the explicit computation of forces and without the step size restrictions required when integrating equations of motion. We employ a size- and shape-dependent damping of cluster movements, motivated by collective hydrodynamic effects neglected in simple implementations of Brownian dynamics. We discuss the virtual-move algorithm in the context of other Monte Carlo cluster-move schemes and demonstrate its utility by applying it to a model of biological self-assembly.

We thank Edward H. Feng, Matthew B. Francis, Michael Hagan, Robert L. Jack, Chad D. Paavola, Sander Pronk, and Jonathan D. Trent for important discussions. The similarity in name between our algorithm and the virtual-move parallel tempering scheme of Ref. 33 is due to our lack of imagination.

I. INTRODUCTION

II. A “VIRTUAL-MOVE” MONTE CARLO CLUSTER ALGORITHM

A. Making collective moves

B. Making collective moves according to potential energy gradients

III. AVOIDING UNPHYSICAL KINETIC TRAPS IN THE FACE OF STRONG INTERACTIONS

IV. APPLICATION OF VMMC TO SELF-ASSEMBLY

V. CONCLUSIONS

### Key Topics

- Diffusion
- 24.0
- Brownian dynamics
- 21.0
- Polymers
- 19.0
- Self assembly
- 19.0
- Cluster dynamics
- 12.0

## Figures

(Color online) An illustration of the difficulty encountered by standard Monte Carlo methods in the presence of very strong, short-ranged interactions. Here only particles and interact and do so strongly. Single-particle Monte Carlo schemes generate a diffusion of the dimer via relative moves of and . Because such moves [e.g., ] are suppressed by the exponential of the resulting energy change, collective modes of motion are under-represented for any acceptance rate less than unity. Moving particles and simultaneously can restore collective motion. We shall discuss how this can be done in order to approximate a realistic dynamics.

(Color online) An illustration of the difficulty encountered by standard Monte Carlo methods in the presence of very strong, short-ranged interactions. Here only particles and interact and do so strongly. Single-particle Monte Carlo schemes generate a diffusion of the dimer via relative moves of and . Because such moves [e.g., ] are suppressed by the exponential of the resulting energy change, collective modes of motion are under-represented for any acceptance rate less than unity. Moving particles and simultaneously can restore collective motion. We shall discuss how this can be done in order to approximate a realistic dynamics.

(Color online) An illustrative collective move. We define a pseudocluster (shaded particles) using an iterative linking scheme (see text). A particular realization of links is denoted by bold black lines. We subject the pseudocluster to a rotation or a translation (or both). Here a proposed translation moves from contact with cluster in state to contact with cluster in state . The interface between the pseudocluster and its environment in state is labeled and is defined by the set of all pairwise interactions between white and shaded particles.

(Color online) An illustrative collective move. We define a pseudocluster (shaded particles) using an iterative linking scheme (see text). A particular realization of links is denoted by bold black lines. We subject the pseudocluster to a rotation or a translation (or both). Here a proposed translation moves from contact with cluster in state to contact with cluster in state . The interface between the pseudocluster and its environment in state is labeled and is defined by the set of all pairwise interactions between white and shaded particles.

(Color online) An illustration of virtual-move Monte Carlo. Starting from state the seed particle is assigned a virtual move map, denoted by a thin arrow (a). We propose a link between and by calculating the energy difference of the bond before and after the virtual move of [(b) and (c)]. In this example a link forms, and adopts the virtual move map of ; the latter is returned to its original position (d). We iterate this procedure until no more links form (e) and displace all linked particles simultaneously (bold arrows denote real, not virtual, moves): (f) . The new configuration is proposed as the final state, and the Monte Carlo acceptance probability is evaluated.

(Color online) An illustration of virtual-move Monte Carlo. Starting from state the seed particle is assigned a virtual move map, denoted by a thin arrow (a). We propose a link between and by calculating the energy difference of the bond before and after the virtual move of [(b) and (c)]. In this example a link forms, and adopts the virtual move map of ; the latter is returned to its original position (d). We iterate this procedure until no more links form (e) and displace all linked particles simultaneously (bold arrows denote real, not virtual, moves): (f) . The new configuration is proposed as the final state, and the Monte Carlo acceptance probability is evaluated.

(Color online) Ensuring reversibility for cluster moves. We require that the likelihood of a given particle, say , “pushing” (f) and “pulling” its host cluster is identical. To do so, every time we link two particles we record the likelihood that the link would have formed had the reverse move been proposed. In the reverse move, the linking particle translates in the opposite direction (and/or rotates with the opposite sense). We account for this probability difference via the overall acceptance rate, ensuring superdetailed balance. The acceptance rate for such collective moves is high if the binding energy of particles is high (the regime of interest) and vanishes in the limit of vanishing interaction strength (Ref. 25). We require that the reverse move be initiated by the seed particle of the forward move, because moves having distinct seeds (compare and ) cannot in general balance each other.

(Color online) Ensuring reversibility for cluster moves. We require that the likelihood of a given particle, say , “pushing” (f) and “pulling” its host cluster is identical. To do so, every time we link two particles we record the likelihood that the link would have formed had the reverse move been proposed. In the reverse move, the linking particle translates in the opposite direction (and/or rotates with the opposite sense). We account for this probability difference via the overall acceptance rate, ensuring superdetailed balance. The acceptance rate for such collective moves is high if the binding energy of particles is high (the regime of interest) and vanishes in the limit of vanishing interaction strength (Ref. 25). We require that the reverse move be initiated by the seed particle of the forward move, because moves having distinct seeds (compare and ) cannot in general balance each other.

(Color online) Generalizing VMMC to allow distinct real and virtual moves. From starting state the seed particle is assigned a virtual move defining a leftward translation and a rotation about the axis through its center (here particles are confined to the plane). In this example the link forms but and do not; we therefore record the link-forming probability and the link-failure probabilities [frames (1) and ]. Frame (see text) shows the intermediate (virtual) state for the forward move. We then return the chosen pseudocluster to its original position (not shown) and execute our chosen real move. Our real move consists of the translational and rotational components of the virtual move with the rotation diminished by the square root of the cluster’s moment of inertia about the relevant axis. Here this move yields proposed final state . From we record the probability of making link by executing the reverse *virtual* move, , and the probabilities of failing to form links and [frames (2) and ]. Frame shows the intermediate (virtual) state for the reverse move. The acceptance probability for this procedure is given by Eq. (17).

(Color online) Generalizing VMMC to allow distinct real and virtual moves. From starting state the seed particle is assigned a virtual move defining a leftward translation and a rotation about the axis through its center (here particles are confined to the plane). In this example the link forms but and do not; we therefore record the link-forming probability and the link-failure probabilities [frames (1) and ]. Frame (see text) shows the intermediate (virtual) state for the forward move. We then return the chosen pseudocluster to its original position (not shown) and execute our chosen real move. Our real move consists of the translational and rotational components of the virtual move with the rotation diminished by the square root of the cluster’s moment of inertia about the relevant axis. Here this move yields proposed final state . From we record the probability of making link by executing the reverse *virtual* move, , and the probabilities of failing to form links and [frames (2) and ]. Frame shows the intermediate (virtual) state for the reverse move. The acceptance probability for this procedure is given by Eq. (17).

(Color online) Configurations as a function of time for very attractive disks evolved according to Brownian dynamics (BD), virtual-move Monte Carlo (VMMC), and single-particle Monte Carlo (SPMC) protocols. Cluster sizes and shapes as a function of time are similar under BD and VMMC algorithms. VMMC and SPMC simulations were performed using the same distribution of basic displacements. The acceptance rate for sequential moves of single particles in this regime is low enough to suppress, unphysically, collective modes of relaxation. The computational efficiency of VMMC exceeds that of the BD protocol by more than an order of magnitude.

(Color online) Configurations as a function of time for very attractive disks evolved according to Brownian dynamics (BD), virtual-move Monte Carlo (VMMC), and single-particle Monte Carlo (SPMC) protocols. Cluster sizes and shapes as a function of time are similar under BD and VMMC algorithms. VMMC and SPMC simulations were performed using the same distribution of basic displacements. The acceptance rate for sequential moves of single particles in this regime is low enough to suppress, unphysically, collective modes of relaxation. The computational efficiency of VMMC exceeds that of the BD protocol by more than an order of magnitude.

(Color) Geometry for a schematic model of chaperonin self-assembly (Refs. 9–11). We use a short-ranged, anisotropic pair potential to mimic the tendency of chaperonins to bind equator to equator. We show chaperonin structure determined by homology (image courtesy of Matthew B. Francis and Chad D. Paavola) together with the angles relevant for our chosen interaction. The angle between orientation vectors is ; the angles between orientation vectors and the interunit separation vector are and .

(Color) Geometry for a schematic model of chaperonin self-assembly (Refs. 9–11). We use a short-ranged, anisotropic pair potential to mimic the tendency of chaperonins to bind equator to equator. We show chaperonin structure determined by homology (image courtesy of Matthew B. Francis and Chad D. Paavola) together with the angles relevant for our chosen interaction. The angle between orientation vectors is ; the angles between orientation vectors and the interunit separation vector are and .

(Color online) Kinetic phase diagram for our schematic chaperonin system, evolved using the VMMC algorithm with hydrodynamic damping. We indicate regions of “no assembly” (open squares), “good assembly” (circles), and “bad assembly” (closed squares); see text. A good assembly–bad assembly pair indicates that intermediate building blocks are well formed, but subsequent multiparticle collisions induce kinetic frustration. Configurations corresponding to parameter sets (a), (b), and (c) are shown in Fig. 9–11.

(Color online) Kinetic phase diagram for our schematic chaperonin system, evolved using the VMMC algorithm with hydrodynamic damping. We indicate regions of “no assembly” (open squares), “good assembly” (circles), and “bad assembly” (closed squares); see text. A good assembly–bad assembly pair indicates that intermediate building blocks are well formed, but subsequent multiparticle collisions induce kinetic frustration. Configurations corresponding to parameter sets (a), (b), and (c) are shown in Fig. 9–11.

(Color online) Configurations obtained after (top) and (bottom) VMMC sweeps for our schematic chaperonin system of spheres with sticky equators of strength [parameter set (a) of Fig. 8]. At this attraction strength nucleation is sufficiently rare that self-assembly proceeds by the binding of monomers to a single sheet. The resulting structure is relatively well-formed.

(Color online) Configurations obtained after (top) and (bottom) VMMC sweeps for our schematic chaperonin system of spheres with sticky equators of strength [parameter set (a) of Fig. 8]. At this attraction strength nucleation is sufficiently rare that self-assembly proceeds by the binding of monomers to a single sheet. The resulting structure is relatively well-formed.

(Color online) Configurations obtained using the VMMC algorithm applied to the chaperonin system with equatorial interaction strength [parameter set (b) of Fig. 8]. A modest increase in attraction strength promotes nucleation to a considerable degree. Assembly proceeds both by the addition of monomers to single sheets (top panel, sweeps) and via collisions of multiparticle sheets (bottom panel, sweeps). Sheets often collide awkwardly, producing ill-formed structures. Here relaxation of these structures is sufficiently rapid that assembly is still “good.”

(Color online) Configurations obtained using the VMMC algorithm applied to the chaperonin system with equatorial interaction strength [parameter set (b) of Fig. 8]. A modest increase in attraction strength promotes nucleation to a considerable degree. Assembly proceeds both by the addition of monomers to single sheets (top panel, sweeps) and via collisions of multiparticle sheets (bottom panel, sweeps). Sheets often collide awkwardly, producing ill-formed structures. Here relaxation of these structures is sufficiently rapid that assembly is still “good.”

(Color online) Configuration obtained via VMMC for our schematic chaperonin system with equatorial interaction strength [parameter set (c) of Fig. 8]. Nucleation is so rapid that multiparticle binding events occur frequently. The resulting structures fail to relax before encountering other such structures, leading to aggregates trapped far from equilibrium. (Top) Excluded-volume view. (Bottom) Bond view.

(Color online) Configuration obtained via VMMC for our schematic chaperonin system with equatorial interaction strength [parameter set (c) of Fig. 8]. Nucleation is so rapid that multiparticle binding events occur frequently. The resulting structures fail to relax before encountering other such structures, leading to aggregates trapped far from equilibrium. (Top) Excluded-volume view. (Bottom) Bond view.

(Color online) Configurations obtained via a single-particle Monte Carlo algorithm for our schematic chaperonin system with , and (top to bottom) equatorial interactions , 7, and . For the displacement distribution employed here (uniform, with maximum displacement ) single-particle moves strongly suppress collective modes of motion and very few multiparticle collisions take place. For all attraction strengths structures are well formed, even though VMMC results indicate that for the two larger values of a collective dynamics results in the formation of nonoptimal aggregates via multiparticle collisions (see Figs. 9–11).

(Color online) Configurations obtained via a single-particle Monte Carlo algorithm for our schematic chaperonin system with , and (top to bottom) equatorial interactions , 7, and . For the displacement distribution employed here (uniform, with maximum displacement ) single-particle moves strongly suppress collective modes of motion and very few multiparticle collisions take place. For all attraction strengths structures are well formed, even though VMMC results indicate that for the two larger values of a collective dynamics results in the formation of nonoptimal aggregates via multiparticle collisions (see Figs. 9–11).

(Color online) Configurations from representative trajectories of our test system (see text) at times . (Left panels) Single-particle moves plus rejection-free cluster algorithm, leading to a configuration trapped far from equilibrium. (Right panels) Single-particle moves plus cleaving algorithm. The latter circumvents the kinetic traps associated with rarely testing strong intercluster bonds.

(Color online) Configurations from representative trajectories of our test system (see text) at times . (Left panels) Single-particle moves plus rejection-free cluster algorithm, leading to a configuration trapped far from equilibrium. (Right panels) Single-particle moves plus cleaving algorithm. The latter circumvents the kinetic traps associated with rarely testing strong intercluster bonds.

(Color online) Acceptance rates (vertical axis) as a function of maximum displacement (horizontal axis) for two particles interacting with the Lennard-Jones-esque potential discussed in Sec. III. Simulations in Sec. III were performed using . The solid line shows the ratio of actual cluster displacements to the displacement predicted on the grounds of cluster-move proposal rate, . The deviation from unity is due to the superdetailed balance factor. This ratio is close to unity for the displacement scale we use, indicating that little unwanted suppression of collective motion occurs. We show also the fraction of moves for which relative particle displacements are accepted (, dot-dashed line) and the fraction of moves for which the particles are displaced collectively (, curved dotted line).

(Color online) Acceptance rates (vertical axis) as a function of maximum displacement (horizontal axis) for two particles interacting with the Lennard-Jones-esque potential discussed in Sec. III. Simulations in Sec. III were performed using . The solid line shows the ratio of actual cluster displacements to the displacement predicted on the grounds of cluster-move proposal rate, . The deviation from unity is due to the superdetailed balance factor. This ratio is close to unity for the displacement scale we use, indicating that little unwanted suppression of collective motion occurs. We show also the fraction of moves for which relative particle displacements are accepted (, dot-dashed line) and the fraction of moves for which the particles are displaced collectively (, curved dotted line).

Article metrics loading...

Full text loading...

Commenting has been disabled for this content