A particle-based Ising model
Quentin Novinger, Antonio Suma, Daniel Sigg, Giuseppe Gonnella, Vincenzo Carnevale
AA particle-based Ising model.
Quentin Novinger, Antonio Suma,
2, 1
Daniel Sigg, Giuseppe Gonnella, and Vincenzo Carnevale Institute for Computational Molecular Science, Temple University, Philadelphia, PA 19122, USA Dipartimento di Fisica, Universit`a degli Studi di Bari and INFN,Sezione di Bari, via Amendola 173, 70126 Bari, Italy dPET, Spokane, WA, USA (Dated: September 25, 2020)We characterize equilibrium properties and relaxation dynamics of a two-dimensional lattice con-taining, at each site, two particles connected by a double-well potential (dumbbell). Dumbbells areoriented in the orthogonal direction with respect to the lattice plane and interact with each otherthrough a Lennard-Jones potential truncated at the nearest neighbor distance. We show that thesystem’s equilibrium properties are accurately described by a two-dimensional Ising model with anappropriate coupling constant. Moreover, we characterize the coarsening kinetics by calculating thecluster size as a function of time and compare the results with Monte Carlo simulations based onGlauber or reactive dynamics rate constants. I. INTRODUCTION
Many complex system can be described as a collec-tion of interacting diffusing agents with internal degreesof freedom associated with a discrete set of states. Ex-amples range from bacterial communities, to layers ofcells in tissue development, to biomolecular systems [1–3], where the internal states can represent either emissionof electric or chemical signals, differences in motility, orconformational changes of the agents structure. Impor-tantly, the internal state of the agent/particles may af-fect the type and/or strength of pairwise interactions.This results in non-trivial and often interesting statis-tical properties under equilibrium or out-of-equilibriumconditions.Notable examples of biomolecular systems showingthese features are biological membranes. Cell mem-branes are complex mixtures of amphiphilic molecules,namely lipids, which show a complex phase behavior [4–7]. Of particular relevance are two macroscopic phases:the liquid disordered, L d , and liquid ordered, L o , thatco-exists in ternary mixtures containing approximatelyequal amounts of cholesterol, unsaturated and saturatedlipids. At the de-mixing transition, while unsaturatedlipids partition into a L d phase, saturated lipids mostlysegregate in cholosterol-rich, hexatically-ordered ( L o ) do-mains in which the hydrophobic tails are in the extendedconformation and give rise to a thicker layer compared tothe surrounding L d phase [8]. Interesting and still partlyunanswered questions then concern the role of integralmembrane proteins, which may have different affinitiesfor the two macroscopic phases: does the presence ofproteins change the phase behavior of the membrane? Istheir lateral distribution or their internal conformationaltransitions affected by fluctuations in lipid composition?To address these questions, it is customary to define astatistical field that describes locally the two macroscopicphases L o and L d , mapping their two components ontothe celebrated Ising model (see, e.g. [9]) with Hamilto- nian: H ( σ ) = − (cid:88)
We consider a two-dimensional triangular lattice com-posed of N d pairs of particles (dumbbells), for a totalof N = 2 N d particles, oriented along the direction or-thogonal to the lattice plane (z). The primitive vec-tors of the triangular lattice are n = [ σ min , ,
0] and n = [ σ min , √ σ min , σ min = 2 σ is the min-imum of the Lennard-Jones potential described below.This lattice implies 6 neighbors per dumbbell (Fig.1A).The first layer of particles composing a dumbbell (type1 particles), with Cartesian coordinates r i = [ x, y, z ]( i = 1 , . . . N d ), is constrained on the plane z = 0, whileparticles in the second layer (type 2 particles), with coor-dinates r i ( i = 1 , . . . N d ), can move along z , see Fig.1B-C. The x and y coordinates of all particles are fixed. Type2 particles interact with particles of the same type witha Lennard-Jones potential: V ( r ) = (cid:40) (cid:15) (( σr ) − ( σr ) ) − δ ≤ σ cut x > σ cut , (2)where σ cut = 1 . σ is the cutoff distance, and δ =4 (cid:15) (( σσ cut ) − ( σσ cut ) ) ensures that the potential is contin-uous at σ cut . The potential is thus attractive between theminimum σ min and the cutoff, and repulsive for r < σ min (Fig.1D). Since the coordinates of type 1 particles arefixed over time, no interaction potential is defined be-tween them. When dumbbells are allowed to diffuse inthe x-y plane (a case that we do not discuss in this pa-per) an interaction potential analogous in form to V isintroduced also between particles of type 1.Particles composing a dumbbell, r i and r i , are con-nected by a quartic bond potential. This functional formresults in two energy minima, which can be approxi-mately mapped onto the spin states of an Ising model(Fig.1E). In particular, the bond potential is defined as: V ( r ) = (cid:15) (cid:104) a ( r − c ) − b ( r − c ) (cid:105) + δ (3) where a, b are positive control parameters. The two po-tential minima are positioned at r min = c ± (cid:113) b a , thebarrier height is h = b a and the constant δ = (cid:15) b a conveniently sets the minimum of the potential to zero(Fig.1F). We choose a = 1699 . b = 247 .
5, such thatthe distance between the two minima is 0 . σ and thebarrier is h = 9 (cid:15) ; under these conditions, the two po-tential wells are sufficiently narrow and well separatedthat a meaningful mapping can be established betweenthe z coordinate of the particle and the two discrete spinstates. The constant c = 1 . σ , ensures that the twominima are at r = 1 σ and r = 1 . σ (Fig. 1F). Thisspecific choice of distances ensures that, when two neigh-boring particles of type 2 are in different minima, theirinteraction energy corresponds to the Lennard-Jones po-tential at its inflection point (Fig. 1D-E), located at σ inf = ( ) = 1 . σ , i.e. where the attractive forceis at its maximum.For practical reasons, in numerical simulations we con-sidered a modified version of V ( r ) that greatly increasesthe rate of transitions between the two spin states. Wetruncated the potential around the local maximum, be-tween the two values of r corresponding to an energy of (cid:15) ( r = 1 . σ and r = 1 . σ ). We then defined aquadratic potential in this region so that an energy bar-rier of (cid:15) would be present (the total barrier height fromthe minima is thus 2 (cid:15) ). The overall piecewise potentialreads: V mod ( r ) = (cid:40) − (cid:15)d ( r − c ) + 2 (cid:15) r ≤ r ≤ r (cid:15) [ a ( r − c ) − b ( r − c ) ] + δ otherwise (4)where d = 1 / ( r − c ) . Note that the choice of V mod ,which is continuous at r and r (Fig. 1F), is a com-promise between the contrasting requirements of havinga surmountable energy barrier at thermal equilibriumand enforcing a negligible occupancy of the intermedi-ate states (to obtain an effective two-state behavior).In conclusion, the z position of type 2 particles canbe mapped onto a spin variable, with z = 1 σ represent-ing spin s = − z = 1 . σ spin s = +1. The energymaximum is located at z = 1 . σ with a barrier of height2 (cid:15) . The time evolution of the system is analyzed by in-tegrating an equation of motion of the Langevin type: m ¨ r i = − γ ˙ r i − ∂V mod ∂ r i − (cid:88) i (cid:54) = j ∂V ∂ r ij + (cid:112) k B T γ η i ( t ) , (5)where r ij = r i − r j , r i = r i − r i . The sum is referredto the dumbbell index. T and γ are the temperature andfriction, respectively, of the thermal bath in contact withthe system, m the particle’s mass, η i is an uncorrelatedGaussian noise with zero mean and unit variance. γ isset to 0.5, which corresponds to an intermediate regimeof friction (see Sec. IV), in order for spin-flip to occurefficiently over time. σ min A t y p e t y p e xy z B σ cut D V σ min V mod σ cut E V (r)V (r) mod σ s = -1s = +1 C z fixed σ inf r σ1.54σ r σ inf
0 2 4 6 8 10 0.8 1 1.2 1.4 1.6 1.8V mod [ ε ] r [ σ ] F -1-0.5 0 0.5 1 0.8 1 1.2 1.4 1.6 1.8 V [ ε ] r [ σ ] FIG. 1. A) Top view of the system (on the x-y plane), showing the triangular lattice, along with σ cut = 1 . σ ; σ is the particlesdiameter and σ min = 2 σ . B) Side view of the system (x-z plane). Each dumbbell consists of a bottom particle (type 1,grey) and a top particle (type 2, red) whose z coordinate varies. The former are fixed in position, while the latter interactwith potential V . Particles within each dumbbell interact through the potential V mod . C) Instantaneous configurationof the dumbbells obtained by using the VMD software [13]. D) Potential V ( r ), with σ cut , σ min and the inflection point σ inf = ( ) σ highlighted. E) Side view (on the x-z plane) of two interacting dumbbells in two different minima of V mod .The minima are located at z = σ (spin s = −
1) and z = 1 . σ ( s = +1). The type 2 particles are at a distance of σ inf . F)Dumbbell potential V (purple), and its modified version V mod (blue) with the quadratic function used between r = 1 . σ and r = 1 . σ . The two minima correspond to different z values of type 2 particle. B. Simulations setup and numerical methods
The potentials and the equation of motion describedin Sec. II A were implemented in a modified version ofthe LAMMPS software [14]. The reference units are thestandard Lennard-Jones reduced units m , σ and (cid:15) , allset to unity, as well as k B = 1. Thus, the time units are τ LJ = (cid:113) mσ (cid:15) . We considered temperatures in the range T ∈ [0 . , . N d = 900 dumbbells. The x, y, z positions of type 1 particles were initialized at the latticesites as prescribed in the previous section, while all type2 particles were initialized to z = 1 (spin -1), and thenallowed to equilibrate.The dynamics of single spin-flip, described and char-acterized in Sec. IV A, was studied in a 6 by 6 lattice(36 dumbbells), and the small cluster spin-flip dynamics(Sec. IV B) was studied in a system at least twice as largeas the cluster-side.For the purpose of studying the cluster growth(Sec. IV C), we simulated a 1024 by 1024 system for a total of N d = 1024 dumbbells. Type 1 particles wereinitialized as before, while type 2 particles had the z position randomly assigned to the z = 1 or z = 1 . v is used to coarse-grain the config-uration space and a short-ranged (Gaussian) repulsivepotential, V ( v ( t )), centered on the instantaneous valueof the collective variable is added at fixed time-intervalsto the unperturbed Hamiltonian of the system. It canbe shown that the time integral of the biasing potentialconverges asymptotically to the underlying free energysurface [19]. The collective variable of choice here is theaverage magnetization per dumbbell, defined as M = 1 N d (cid:112) b/ a (cid:88) i ( z i − c ) , (6)where the sum over i is referred to the dumbbell in-dex, and the magnetization distribution is concentratedin the interval M ∈ [ − ,
1] (altough, strictly, M is un-bound). This rescaling allowed for easier comparison tothe properties of the Ising model. Gaussians potentialswere added every 1,000 timesteps for all temperatures,while the Gaussian height was set to 0 . (cid:15) and theGaussian width to 0.001. The latter was chosen to besmaller than the collective variable standard deviationat all temperatures. The potentials of mean force (pmfs)were computed by averaging the final third of the trajec-tory, sampled every 10 steps.Metadynamics simulations were run for a maximumof 5 · timesteps. Simulations to compute the aver-age magnetization and the heat-capacity were run with-out metadynamics for a maximum of 2 · timesteps.Simulations for sampling single spin events lasted 10 timesteps. The largest system used to characterize clus-ter growth was run for about 2 · timesteps. C. Kramers formula for rate constants
In Sec. IV A, we compared the molecular dynamics(MD) spin flip kinetics with the Kramers’ intermediate-friction formula [20], which approximates the rate of es-cape of a Brownian particle from a free energy well. Theapproximate rate reads: k = ω R πω b − γ m + (cid:115)(cid:18) γ m (cid:19) + ω b exp (cid:18) − ∆ FT (cid:19) , (7)where γ is the friction coefficient and F ( M ) is the pmf.Parameters derived from the pmf are the reactant wellfrequency ω R = (cid:112) F (cid:48)(cid:48) ( M R ) /m , the barrier frequency ω b = (cid:112) − F (cid:48)(cid:48) ( M b ) /m , and the barrier height ∆ F = F ( M b ) − F ( M R ). Eq. 7 approaches the well-known largefriction form for large γ > mω b . D. Dynamical Ising model
To compare the kinetics of MD with that of the Isingmodel, we considered a two-dimensional triangular latticeof spins governed by a master equation using reactive orGlauber dynamics. Reactive dynamics [21] are character-ized by an Arrhenius-type rate constant k = νe − ∆ E/ T ,whereas in Glauber dynamics [22] flip rates are saturable,described by k = νe − ∆ E/ T / ( e − ∆ E/ T + e ∆ E/ T ). Therate factor ν determines the time scale of kinetics, while∆ E = 4 J (3 − u ) is the change in system energy for a spinto flip from − u = 0 , .., −
1) are obtained by changing the sign of ∆ E .The extended triangular lattice in the shape of a paral-lelogram can be mapped to the traditional square Isinggrid with 2 diagonal interactions in addition to the usual4 orthogonal interactions for a total of 6 neighboring spininteractions. Ising model simulations were performed inreal time using the Gillespie algorithm for kinetic Monte Carlo (MC) [23], in which waiting times between spinflips are given by − ln r/ (cid:80) i k i , where r is a uniform ran-dom number between 0 and 1, and (cid:80) i k i is the sum oftransition rates taken over all spins i . A second ran-dom number determines which spin is flipped weightedby k i / (cid:80) i k i . The Gillespie algorithm is a true kineticMonte Carlo simulation with exponentially-distributedwaiting times between transition events. E. Computing cluster size
In Sec. IV C we computed the clusters size through thestructure factor. Here we describe the numerical imple-mentation.The spin were first discretized to -1 or 1 (using z = c asmid point). The bounding square rectangle was changedinto a parallelogram conforming to the skew-geometry ofthe primitive cell. For each frame, the triangular latticewas digitized onto a √ N d × √ N d square matrix, usingthe same mapping of the Monte Carlo simulations. Thedigitized matrix was used to compute fast Fourier trans-form (FFT), with discrete vectors in real and reciprocalspace ranging from 1 to √ N d and −√ N d / √ N d / k = (cid:113) k + (1 / k − k ) , (8)with k i the coordinate in reciprocal space, in order toconform to the original Euclidean distance in the trian-gular lattice. The average cluster size L ( t ) was obtainedfrom: L ( t ) √ N σ min = (cid:80) k h ( k, t ) /n ( k ) (cid:80) k h ( k, t ) k/n ( k ) ≈ π (cid:80) k h ( k, t ) /n ( k ) (cid:80) k h ( k, t ) , (9)where we approximate n ( k ) ≈ πk , h ( k ) are the summedbin values of the structure factor in radial intervals ( k − . , k + 0 .
5] and n ( k ) are the bin counts. III. EQUILIBRIUM PROPERTIES
Here we show that the properties of this collection ofdumbbells, as defined in Sec. II, closely match those ofthe Ising model. In particular, we will show how themagnetization and the heat capacity vary as a functionof the temperature T . A. Phase diagram
We performed metadynamics simulations of the sys-tem, for temperatures T ∈ [0 . , . zz FIG. 2. Equilibrium configurations sampled via molecular dynamics at six different temperatures, T =0 . , . , . , . , . T = 0 . z coordinate. each T the free energy as a function of the average mag-netization M , as described in Sec. II. We found that,upon increasing the temperature, the system transitionsfrom a fully ordered configuration in which all spins arealigned, to a configuration in which the average magneti-zation is zero. This is apparent from Fig. 2 showing themost probable conformations for six temperatures. At T = 0 . , . , . , .
33 the conformations are magne-tized and at T = 0 .
45 the conformation is clearly dis-ordered. Instantaneous configurations are rendered byshowing type 2 particles from a top view colored accord-ing their z coordinate.Fig. 3 shows the free energy as a function of M for T =0 . , . , . , . , . , . , .
36. Below T = 0 . k B T . At T = 0 . k B T allowing the system toexplore the two oppositely aligned ferromagnetic states,a sign that the temperature is close to the critical one.By further increasing the temperature, the free energyprofile flattens out completely and starts to develop asingle minimum at M = 0, meaning that the system isin the disordered state.Fig. 4 shows the magnetization values measured fromthe positions of the free energy minima for tempera-tures below T = 0 .
33 (red curve), as well as the aver-age magnetization values obtained from long simulationsperformed without metadynamics. Consistent with the -3-2-1 0 1 2 3 4 5 -1 -0.5 0 0.5 1 F r ee E ne r g y [ ε ] M 0.310.320.330.3350.340.350.36
FIG. 3. Free energy as a function of the average magnetization M for 6 temperatures (see key), computed using metadynam-ics. Two distinct minima are visible with an energy barrierbetween them larger than 2 k B T up until T = 0 . insight obtained from the free energy profiles, we observethat around T ∼ .
33 the magnetization quickly dropsoff to M = 0. Deviations from this value at larger tem-peratures are due to limited sampling (see error bars). C v [ k B ] T [ ε /k B ] M FIG. 4. Magnetization and heat capacity as a function ofthe temperature, obtained from unbiased simulations. Thered bar identifies the region where the critical temperature islocated.
B. Heat capacity
The heat capacity was calculated from unbiased tra-jectories as C v = < E > − < E > k B T , (10)where E is the system’s total potential energy. The vari-ance of E was obtained using the blocking analysis inorder to have uncorrelated data (the maximum decorre-lation time observed is 10 τ LJ ).Results are shown in Fig. 4. Similar to the Ising model,a peak in the specific heat capacity around the criticaltemperature (between T c ∈ [0 . , . k B T thc J = 3 . J is the couplingconstant (Eq. 1). In our case, 2 J corresponds to thedifference in interaction energy between the “parallel”and “anti-parallel” case: 2 J = [ V ( σ inf ) − V ( σ min )] =0 . T thc = 0 . IV. DYNAMICAL PROPERTIES
From a dynamical point of view, the trajectory ob-tained from integration of a Langevin-type equation ofmotion is conceptually distinct from the time series ob-tained from a Monte Carlo (MC) simulation. Thus it isinteresting to investigate the dynamical differences be- tween MC and molecular dynamics for these Ising-likesystems.
A. Single spin-flip dynamics
First, we characterize the spin-flip dynamics for a sin-gle spin. In particular, we focus on cases where the sixneighbor spins have their magnetization fixed, Fig 5A.This is obtained by constraining the positions to z > c for s = +1 and z < c for s = −
1, using for each spin arepulsive quadratic potential at z = c with force constant k = 10 (cid:15)/σ . Seven cases are possible, with a total of 0to 6 neighbor spins with opposite magnetization to thecentral one that has to flip (and a complementary num-ber of 6 to 0 neighbor spins with same magnetization asthe central one). The temperature will be from now onfixed at T = 0 .
185 for MD, which is about 0 . T c and wechoose similarly T /T c = 0 .
56 for MC (we use in this casethe theoretical value for T c ).Fig 5B-C show a single spin-flip trajectory over timewith three neighbors with same magnetization and threewith opposite one. The value of z , and consequently of M , varies continuously in time. Fig 5D shows the freeenergy profile, as a function of z , of a spin with 0 to 3neighbor spins with negative magnetization. The othercases have a free energy which is the mirror image around z = c .Fig. 5E shows the rate (inverse time) for the centralspin to invert its magnetization, for the seven possiblecases. The flipping-rate increases exponentially, as ex-pected due to the linear decrease in the energy barrierobserved in Fig 5D. This results are in agreement with theprediction obtained from Kramers’ intermediate-frictionrate constant formula (Eq. 7) using MD-derived pmfs tocompute the formula parameters.Fig. 5F shows how the flipping rate with 3 neighborspins of opposite magnetization is affected with changing γ , along with Eq. 7. Indeed, the chosen value of γ = 0 . ν , J ,and T with arbitrary units, the Ising rate constants forthe reactive dynamics can be aligned with MD-derivedflip rates (Fig. 5E). The Glauber rates, obtained usingthe same ν as the reactive dynamics, underestimate reac-tion rates for all conditions (Fig. 5E). The two rates cancoincide for the specific case of equal number (3) of oppo-sitely magnetized neighbor particles by scaling Glauberrates by a factor of 2. B. Cluster-flip dynamics
Next, we characterized the amount of time needed foran hexagonal cluster to completely invert its magnetiza-tion when surrounded by a bulk of spins with oppositemagnetization (Fig. 6A). Note that only spins inside the z [ σ ] t[ τ LJ ] ABC D
Single spin-flip dynamics E F [ ε ] z [ σ ] -6 -5 -4 -3 -2 -1
0 1 2 3 4 5 6 r a t e /t [ τ L J - ] Number of neighbors with opposite magnetizationMDGlauber MCReactive MCKramers pred. -4 -3 -2
1 10 r a t e /t [ / τ L J ] γ [m/ τ LJ ] MDKramers pred. z [ σ ] t[ τ LJ ] F FIG. 5. A) Illustration of the setup used to study single spin-flip dynamics. The central spin is allowed to flip, while thesix neighbor spins have fixed spins, see text for detail. B) Time evolution of z for the central spin, with setup as in A). C)Zoom-in of B) over a smaller time interval. D) Free energy profile, as a function of z , with 0 to 3 neighbor spins with negativemagnetization. E) Rate of central spin-flip with 0 to 6 neighbor spins with opposite magnetization, along with Glauber andreactive dynamics MC rates adjusted to fit the MD, see text. In green, Kramers prediction based on Eq. 7. F) Rate of centralspin-flip with same setup as A), as function of γ . The temperature is fixed to T = 0 . cluster were considered when monitoring the magnetiza-tion.Fig.6B reports on the cluster-flip time. The time scaleslinearly with the cluster size for sufficiently large clusters,while a non-trivial dependency is observed for low clus-ter sizes, which is dependent on the chosen dynamics.Although the rates between MD and reactive dynamicswere perfectly matched (the curves start from the samepoint) the cluster-flip times do not match for large clus-ters. The ratio between the two slopes is 2 .
75. Instead,for Glauber MC dynamics these times coincide for largeclusters.
C. Cluster growth kinetics
Last, we characterize the time evolution of cluster size L ( t ) by performing a similar analysis to the one reportedin ref. [24]. Fig. 7A shows the evolution of the 1024 spin system starting from a random configuration for thethree different dynamics considered.We computed the structure factor in the recipro-cal space defined by the contravariant vectors b = [ σ min , − √ σ min ,
0] and b = [0 , √ σ min , n · b = n · b = 0. The structure factor is the spectraldensity of the spin matrix, where the spins have beenconsidered in their discretized version s : S F ( k ) = | (cid:88) r ( s ( r ) − (cid:104) s (cid:105) ) e − πi k · r | , (11)where r are the positions in real space ranging from[1 , . . . , √ N d ] n i , k the ones in the reciprocal space, rang-ing between [ −√ N d / , . . . , √ N d / b i / √ N d ) and (cid:104) s (cid:105) isan average over all configuration’s spins. For each frame,we computed S F ( k , t ) and the average cluster size L ( t )was obtained from: L ( t ) = (cid:82) dk S F ( k , t ) /k (cid:82) dk S F ( k , t ) . (12)For the actual numerical implementation details, seeSec. II E.Fig. 7B shows L ( t ), along with a power law (dottedline) with exponent 0.5, which is expected when the or-der parameter is non-conserved as in model A dynamics[25], see also Supplementary Movie S1 (system with 100 spins, T = 0 .
185 starting from a random configuration).
1 10 100 t[ τ L J ] Cluster size (number of spins)MDGlauber MCReactive MC
Cluster-flip dynamics AB FIG. 6. A) Illustration of the setup used to study cluster-flipdynamics. A single cluster is placed in the center of a systemwith opposite magnetization. B) Time required to invert thecluster’s magnetization, as a function of the cluster size, interm of spins number, for MD and reactive and Glauber MC.
Consistent with the observations of Sec. IV B, L ( t ) of MDand Glauber MC overlaps well at long times, while reac-tive dynamics MC does not. The latter overlaps whenthe times are rescaled by the factor 2 .
75 (cyan curve), asmeasured in Sec. IV B.Fig. 7C shows the growth exponent α , L ( t ) ∼ t /α asa function of time, computed as [26]1 α ( t ) = d ln L ( t ) d ln t . (13)The exponent slowly saturates to 0 . T >
0) temperatures [26].Finally, Fig. 7E-F show the system’s energy E ( t ) mea-sured using the Ising model Hamiltonian 1, after spindiscretization for MD configurations, and the total sys- tem’s magnetization | M ( t ) | . The former saturates to apower law (dotted line) with exponent -0.5. V. DISCUSSION
In this work, we investigated a two-dimensional systemof pairs of particles (dumbbells) connected by a double-well potential. This allows the dumbbells to fluctuatebetween two states, which can be interpreted as two spinstates s i = ±
1. Importantly, dumbbells interact witheach other through a Lennard-Jones potential with aninteraction energy that discriminates between “parallel”and “anti-parallel” configurations. In this way, each spincan influence its closest neighbors (six for the triangularlattice considered here). All spins evolve simultaneouslyunder the action of the interaction potential as prescribedby a Langevin-type equation of motion.First, we characterized equilibrium properties, focus-ing on magnetization and heat capacity: the observed be-havior follows closely what expected for an Ising model.Second, we characterized several dynamical properties,such as single spin- and cluster-flip, and the clustergrowth kinetics starting from a random initial configu-ration. We highlighted a rich dynamical behavior, whichdiffers to some extent from two of the most widely con-sidered dynamics, namely reactive and Glauber dynamicsMonte Carlo.Future work will investigate the most general case inwhich dumbbells are allowed to diffuse in the x, y plane,as it occurs in biological systems such as lipid mem-branes. We expect the phase diagram of this off-latticecase to show a non-trivial combination of vapor-liquidand liquid-solid transitions typical of both Lennard-Jonessystems and ferromagnetic order/disorder transitions.
Acknowledgements
This work was funded by the National Insti-tutes of Health (R01GM093290, S10OD020095 andR01GM131048; V.C.), and the National Science Foun-dation through grant IOS-1934848 (V.C.). This researchincludes calculations carried out on Temple University’sHPC resources and thus was supported in part by the Na-tional Science Foundation through major research instru-mentation grant number 1625061 and by the US ArmyResearch Laboratory under contract number W911NF-16-2-0189. [1] Arthur Prindle, Jintao Liu, Munehiro Asally, San Ly,Jordi Garcia-Ojalvo, and Grol M. Sel. Ion channelsenable electrical communication in bacterial communi-ties. Nature, 527(7576):59, 2015. ISSN 1476-4687. doi:10.1038/nature15709.[2] Brian A Camley and Wouter-Jan Rappel. Physical mod-els of collective cell motility: from cell to tissue. Journalof physics D: Applied physics, 50(11):113002, 2017. [3] Ingela Parmryd, David G. Ackerman, and Gerald W.Feigenson. Lipid bilayers: clusters, domains and phases.Essays in Biochemistry, 57:33–42, 02 2015. ISSN 0071-1365. doi:10.1042/bse0570033. URL https://doi.org/10.1042/bse0570033 .[4] Friederike Schmid. Physical mechanisms of micro-and nanodomain formation in multicomponent lipidmembranes. Biochimica et Biophysica Acta (BBA) -
10 100 L ( t ) [ σ ] MDGlauber MCReactive MCShifted Reactive MC M D G l a u b e r M C R e a c t i v e M C t=10 τ LJ t=10 τ LJ t=10 τ LJ A B α ( t ) -1 | M | t[ τ LJ ] CDE E [ ε ] FIG. 7. Cluster growth, starting from an initial random configuration at T = 0 . < T c , as a function of time. A) Threesnapshots are shown at t = 10 , , τ LJ for MD simulations, and the corresponding snapshots for Glauber and reactivedynamics MC. B) Average cluster size L ( t ), C) growth exponent α ( t ), D) Ising model energy and E) absolute value of mag-netization | M | as a function of time for the three different dynamics considered. In cyan, reactive dynamics MC shifted by afactor of 2 .
75, see main text for details. In B) the dashed line shows the power law fit with exponent 0 .
5. In D) it shows apower law with exponent − . Biomembranes, 1859(4):509528, 2017. ISSN 0005-2736.doi:10.1016/j.bbamem.2016.10.021.[5] George A Pantelopulos and John E Straub. Regimesof complex lipid bilayer phases induced by cholesterolconcentration in md simulation. Biophysical journal, 115(11):2167–2178, 2018.[6] Rebecca D Usery, Thais A Enoki, Sanjula P Wickra-masinghe, Michael D Weiner, Wen-Chyan Tsai, Mary BKim, Shu Wang, Thomas L Torng, David G Ackerman,Frederick A Heberle, et al. Line tension controls liquid-disordered+ liquid-ordered domain size transition in lipidbilayers. Biophysical journal, 112(7):1431–1443, 2017.[7] Aurelia R Honerkamp-Smith, Benjamin B Machta, andSarah L Keller. Experimental observations of dynamiccritical phenomena in a lipid membrane. Physical reviewletters, 108(26):265702, 2012.[8] Shachi Katira, Kranthi K Mandadapu, SuriyanarayananVaikuntanathan, Berend Smit, and David Chandler. Pre-transition effects mediate forces of assembly betweentransmembrane proteins. Elife, 5:e13150, 2016.[9] Ofer Kimchi, Sarah L. Veatch, and Benjamin B. Machta.Ion channels can be allosterically regulated by mem-brane domains near a de-mixing critical point. Journalof General Physiology, 150:1769–1777, 2018. doi:https://doi.org/10.1085/jgp.201711900.[10] A Pikovsky, A Zaikin, and M A de la Casa. System SizeResonance in Coupled Noisy Systems and in the IsingModel. Physical Review Letters, 88(5):050601, 2002.ISSN 0031-9007. doi:10.1103/physrevlett.88.050601.[11] Renfrey B Potts. Spontaneous magnetization of a trian-gular ising lattice. Physical Review, 88(2):352, 1952.[12] John Stephenson. Isingmodel spin correlations on thetriangular lattice. Journal of Mathematical Physics, 5(8):1009–1024, 1964. doi:10.1063/1.1704202. URL https://doi.org/10.1063/1.1704202 .[13] William Humphrey, Andrew Dalke, and Klaus Schulten.VMD – Visual Molecular Dynamics. Journal of MolecularGraphics, 14:33–38, 1996.[14] S. Plimpton. Fast parallel algorithms for short-rangemolecular dynamics. J. Comp. Phys., 117:1–19, 1995.[15] T Schneider and E Stoll. Molecular-dynamics study ofa three-dimensional one-component model for distortivephase transitions. Physical Review B, 17(3):1302, 1978. [16] Alessandro Laio and Michele Parrinello. Escaping free-energy minima. Proceedings of the National Academyof Sciences, 99(20):12562–12566, 2002. ISSN 0027-8424.doi:10.1073/pnas.202427399. URL .[17] In-Rok Oh, Eun-Sang Lee, and YounJoon Jung. A min-imalist model of single molecule spectroscopy in a dy-namic environment studied by metadynamics. Bulletinof the Korean Chemical Society, 33, 03 2012. doi:10.5012/bkcs.2012.33.3.980.[18] G. Fiorin, M. L. Klein, and J. H´enin. Using collec-tive variables to drive molecular dynamics simulations.Mol. Phys., 2013. doi: 10.1080/00268976.2013.813594.[19] James F Dama, Michele Parrinello, and Gregory A Voth.Well-tempered metadynamics converges asymptotically.Physical review letters, 112(24):240602, 2014.[20] H.A. Kramers. Brownian motion in a field of force andthe diffusion model of chemical reactions. Physica, 7(4):284 – 304, 1940.[21] Daniel Sigg, Vincent A Voelz, and Vincenzo Carnevale.Microcanonical coarse-graining of the kinetic ising model.The Journal of Chemical Physics, 152(8):084104, 2020.[22] Roy J. Glauber. Timedependent statistics of the isingmodel. Journal of Mathematical Physics, 4(2):294–307,1963. doi:10.1063/1.1703954. URL https://doi.org/10.1063/1.1703954https://doi.org/10.1063/1.1703954