Potential Energy Landscapes for the 2D XY Model: Minima, Transition States and Pathways
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] N ov DAMTP-2013-67
Potential Energy Landscapes for the 2D XY Model:Minima, Transition States and Pathways
Dhagash Mehta, ∗ Ciaran Hughes, † Mario Schröck, ‡ and David J. Wales § Dept of Mathematics, North Carolina State University, Raleigh, NC 27695, USA. The Department of Applied Mathematics and Theoretical Physics,The University of Cambridge, Clarkson Road, Cambridge CB3 0EH, UK. Institut für Physik, FB Theoretische Physik, Universität Graz, 8010 Graz, Austria. University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, UK.
Abstract
We describe a numerical study of the potential energy landscape for the two-dimensional XY model (withno disorder), considering up to 100 spins and CPU and GPU implementations of local optimization, focusingon minima and saddles of index one (transition states). We examine both periodic and anti-periodic boundaryconditions, and show that the number of stationary points located increases exponentially with increasinglattice size. The corresponding disconnectivity graphs exhibit funneled landscapes; the global minima arereadily located because they exhibit relatively large basins of attraction compared to the higher energy minimaas the lattice size increases. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] . INTRODUCTION Stationary points of a potential energy function, defined as configurations where the gradient vanishes,play a key role in determining many phenomena in physics and chemistry. An extensive framework ofconceptual and computational tools has been developed corresponding to the potential energy landscapeapproach [1–5], with applications to many-body systems as diverse as metallic clusters, biomolecules,structural glass formers, and coarse-grained models of soft and condensed matter. In all these examples,stationary points of a high-dimensional potential energy function are considered. Due to the non-linearnature of the potential energy as a function of coordinates arising in most real world applications,obtaining the stationary points analytically is not feasible. Hence, one has to rely on numerical methodsto obtain the necessary information.In the present contribution we initiate an extensive numerical analysis of the stationary points of awell-known example, the XY model (without any disorder), for a Hamiltonian defined in terms of thepotential energy: H = 1 N d d X j =1 X i [1 − cos( θ i +ˆ µ j − θ i )] , , (1)where d is the dimension of a lattice, ˆ µ j is the d -dimensional unit vector in the j -th direction, i.e. ˆ µ =(1 , , . . . , µ = (0 , , , . . . , i stands for the lattice coordinate ( i , . . . , i d ), and the sum over i represents a sum over all i , . . . , i d each running from 1 to N , and each θ i ∈ ( − π, π ]. Hence d is thedimension of the lattice, and N is the number of sites for each dimension, so the number of θ valuesrequired to specify the configuration is N d . The boundary conditions are given by θ i + N ˆ µ j = ( − k θ i for 1 ≤ j ≤ d , where N is the total number of lattice sites in each dimension, with k = 0 for periodicboundary conditions (PBC) and k = 1 for anti-periodic boundary conditions (APBC). With PBC thereis a global degree of freedom leading to a one-parameter family of solutions, as all the equations areunchanged under θ i → θ i + α, ∀ i , where α is an arbitrary constant angle, reflecting the fact that themodel has global O(2) symmetry. We remove this degree of freedom by fixing one of the variables tozero: θ ( N,N,...,N ) = 0. We have included the factor 1 /N d to facilitate comparisons between systems ofdifferent sizes. In the present contribution we mainly focus on analysis of local minima and the pathwaysbetween them that are mediated by transition states (saddles of index one, with a single negative Hessianeigenvalue [6]). H appears in many different contexts: first, in statistical physics H is known as the XY ModelHamiltonian and is known to exhibit a Kosterlitz-Thouless transition [7]. It describes a system of N classical planar spin variables, where each spin is coupled to its four nearest neighbors on the lattice.This representation is employed in studies of low temperature superconductivity, superfluid helium,hexatic liquid crystals, and Josephson junction arrays. H also corresponds to the lattice Landau gauge2unctional for a compact U (1) lattice gauge theory [8–10], and to the nearest-neighbor Kuramoto modelwith homogeneous frequency, where the stationary points constitute special configurations in phasespace from a non-linear dynamical systems viewpoint [11].The XY model is among the simplest lattice spin models in which an energy landscape approach basedon stationary points of the Hamiltonian in a continuous configuration space is appropriate (unlike, forexample, the Ising model whose configuration space is discrete). Nevertheless, we find that the potentialenergy landscape supports a wide range of interesting features, and proves to be very helpful in analyzingthe characteristic structure, dynamics, and thermodynamics.In Ref. [9, 10] all the stationary points of the one-dimensional XY model were found, includinginteresting classes, such as stationary wave solutions. In Ref. [9, 12, 13], all the stationary pointsfor the one-dimensional model with anti-periodic boundary conditions were characterized. Solving thestationary equations for the XY model in more than one dimension turned out to be a difficult task andhas not been completed to date. In Ref. [9] it was shown how the stationary equations could be viewedas a system of polynomial equations, and the numerical polynomial homotopy continuation (NPHC)method was employed to find all the stationary points for small lattices in two dimensions. This methodwas subsequently used to study the potential energy landscapes of various other models in statisticalmechanics and particle physics [14–23]. In particular, it was employed to study the potential energylandscape of the two-dimensional (2D) XY model [24–26]. In the latter work, along with all the isolatedsolutions for a small 3 × O (2) freedom.In Ref. [25] application of the conjugate gradient method for small N suggested that the number oflocal minima for the 2D XY model would increase exponentially. One of the important results of thecurrent paper is to verify this conjecture and make it more precise by characterizing the landscapes forlarger N values, while improving on the earlier results for the number of minima.All the above-mentioned minimization methods have a common shortcoming because they cannotdeal with even moderately high N (of the order 100 angles). Finding saddles is even harder, and so farthe only available results are for 3 × OPTIM package [27], and a GPU-implementation of the overrelaxation method. These approachescan explore the PELs of the 2D XY model with 100 spins and beyond. In the next two Sections, wedescribe the functionalities of
OPTIM that we have used in our work and the GPU implementation of theoverrelaxation method. 3
I. METHODSA. Geometry Optimization
The
OPTIM program includes a wide variety of geometry optimization tools for characterizing sta-tionary points and the pathways that connect them [27]. The most efficient [28] gradient-only min-imization algorithm implemented in
OPTIM is a modified version of the limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm [29, 30]. Both gradient-only and second derivative-basedeigenvector-following [31, 32] and hybrid eigenvector-following algorithms [33, 34] are available for single-and double-ended [35] transition state searches. Stationary points with any specified Hessian index canalso be located [36].We used
OPTIM to sample minima and transition states for the 2D XY model with both PBC andAPBC, using exclusively single-ended search algorithms. In particular, we refined 500 ,
000 randominitial guesses for all lattice sizes up to N = 10, i.e., a total of 100 spins. In each case there also existdegenerate stationary points related by symmetry operations of the Hamiltonian with θ i → − θ i for all N d angular variables as well as θ i → θ i ± ( π, π, . . . , π ). B. GPU Implementation of Overrelaxation
Another approach that we applied to obtain as many minima as possible employed the (over)relaxationalgorithm exploiting graphics processing units (GPUs), which offer a high level of parallelism and thusenabled us to generate large samples within a practical amount of computer time. The idea of therelaxation algorithm is to sweep over the lattice while optimizing the Hamiltonian locally on each latticesite. Our implementation is based on the code [37].In practice we employed four cards of the NVIDIA Tesla C2070 and launched 1024 thread blocks(i.e. samples) per GPU. We kept the grid size as 1024 blocks per GPU fixed and then cycled over 2 =131072 iterations, resulting in around 0 .
134 billion samples per lattice size. We set the overrelaxationparameter to 1.0, i.e. standard relaxation, to increase the chance of finding minima with small basinsof attraction. For each sample we stored the value of the minimum to which the relaxation algorithmconverged along with the N d corresponding θ -coordinates, and subsequently sorted these values viabitonic sort [38], again accelerated by the GPU. As a stopping criterion we required the largest gradientover all lattice sites to be smaller than 10 − (reduced units). The whole simulation was performed indouble precision to reduce numerical inaccuracies.4 II. RESULTS
We first point out that the global minimum of this model, as it does not have any disorder, is wellknown: for the APBC case, θ i = 0 or π for all N d angular variables are the two global minima of themodel at which H = 0. Similarly, for the PBC case, because we have fixed the O (2) symmetry, theunique global minimum corresponds to θ i = 0 for all angles, again with H = 0. A. Number of Minima and Transition States
In Table I we summarize the number of minima and transitions states located for each N . Here, inaddition to finding minima and transition states for larger lattices, we have also improved our previousresults for the number of minima from Ref. [25] at smaller N with the help of the more powerfulalgorithms. Saddles of index one were only obtained from the OPTIM runs. Since the 2D XY modelpossesses a number of discrete symmetries, as discussed in [25], we also tabulate the number of distinctminima and transition states in the table, i.e. solutions unrelated by symmetries of the Hamiltonian.In contrast to N ≤
8, for larger N our samples will be substantially incomplete, even though we havefound around 5 . N = 16 case. As expected from symmetries of the Hamiltonian[25], for each N with PBC the global minimum is unique, the next minimum is 4-fold degenerate, then N / N is even), then 2 N , 4 N , 2 N ,...., N / N is even), and the highest energyminimum is 4-fold degenerate. For the APBC case the global minimum is 2-fold degenerate, and allother minima are at least N degenerate.The number of minima and transition states located as a function of N are plotted in Figure 1. Theplot clearly shows that the total number of minima (including degeneracies), the number of distinctminima, the total number transition states, and the number of distinct transition states, all increaseroughly exponentially with increasing N , as expected [36, 39]. We also observe abrupt jumps at N = 7and 9 for the PBC case, and at N = 6 for the APBC cases in this plot, though the precise reason ofthis behavior is not clear. It is possible that the jumps are caused by sampling issues, but there couldbe a more subtle explanation; for example, certain lattices for particular values of N may restrict thepossible classes of minima.For both APBC and PBC the energies of the local minima shift towards lower energy as N increases,tending to accumulate near the global minima. This behavior has previously been observed in the 1DXY model with PBC [25]. In this case, the potential energy distribution of the minima has a spike forthe global minimum at H = 0 and a two-fold degeneracy for other minima in H ∈ (0 , certify these solutions using techniques based on Smale’s α -theorem [41]. B. Disconnectivity Graphs
Disconnectivity graphs have provided a particularly useful tool for visualizing potential and free en-ergy landscapes [40, 42–44] in systems ranging from atomic and molecular clusters to soft and condensedmatter and biomolecules [1–3]. In particular, this construction enables the lowest potential or free energybarriers to be faithfully represented, and can help us to understand how observable properties emergefrom features of the landscape [45].To produce a disconnectivity graph we require a kinetic transition network [45–47], which can bedefined by a database of local minima and the transition states that connect them [40]. We then choosea regular energy spacing, ∆ V , and determine how the minima are partitioned into subsets (superbasins[40]) at energies V , V + ∆ V, V + 2∆ V, . . . . These subsets consist of minima that can interconvert viaindex one saddles that lie below the energy threshold. For a high enough threshold all the minima caninterconvert and there is just one superbasin, unless there are infinitely high barriers. As the thresholdenergy decreases the superbasins split apart, and this splitting is represented in the disconnectivity graphby lines connecting subgroups to parent superbasins at the threshold energy above. The superbasinsterminate at the energies of individual local minima, which may be grouped together for degeneratestates related by symmetry operations of the Hamiltonian.The significance of the disconnectivity graph construction stems primarily from the insight it providesinto the global thermodynamics and kinetics of the system in question. For example, if the landscapesupports alternative morphologies separated by a high barrier then we anticipate a separation of relax-ation time scales and associated features in the heat capacity [1–3, 45]. Several limiting cases have beenidentified for the organization of the landscape, distinguishing good ‘structure seeking’ systems, whichexhibit efficient relaxation to the global minimum, from models with glassy characteristics [42]. Thesevisualisations have much in common with the ‘energy lid’ and ‘energy threshold’ approaches of Sibani,Schön, and coworkers [48–52].In the current contribution we have characterized both minima and transition states, which enablesus to construct the first disconnectivity graphs for XY models (Figures 2 and 3). These graphs allcorrespond to the structure expected for efficient relaxation to the global minimum over a wide range of6emperature (or total energy), namely the ‘palm tree’ motif [42]. Locating low-lying minima for these 2DXY models should therefore be relatively straightforward: relaxation following the intrinsic dynamics ofthe system should lead to the global minimum for temperatures of physical interest. This is the patternthat we associate with good structure-seeking systems [2, 4, 42, 45], including ‘magic number’ clusterssuch as buckminsterfullerene, self-assembling mesoscopic structures such as virus capsids, crystallisation,and proteins that fold into functional native states on in vivo time scales.
C. Energy Differences
Experimentally, it is not the absolute value of the energy but rather energy differences that aremeasured. For the 2D XY model with no disorder, we can in principle study dE ik,l = E ik − E il , where E ik is the energy of the k -th index i saddle in order of increasing energy. For the 2D XY model, the globalminimum has energy E ≡ E = 0 yielding dE k, = E k . The sequential energy differences dE ik +1 ,k areparticularly interesting, since all other energy differences can be obtained from them. We plot N vs dE k +1 ,k in Figure 4 and N vs dE k +1 ,k in Figure 5. We find that that the sequential energy differencesdecay towards zero as N increases.This observation is similar to results for the 1D XY model with PBC. There, in the continuum limit,the energies of the local minima are distributed continuously over the range [0 , E = 0, while other energy values are two-fold degenerate. For small N , it appears as if dE k, is decaying towards zero. However, as N becomes large enough, dE k, eventuallystarts to fill in [0 , dE k +1 ,k , that decays towards zero.For every saddle in 1D with PBC, we can build a higher dimensional analogue that has the same energy[25]. Hence, the energy spectrum of the 2D XY model with PBC and no disorder contains at least onecopy of the 1D XY model with PBC. dE k, has been studied in reference [25], where it was found to decay to zero for the smaller latticesizes. The results in Figures 4 and 5, coupled to the previous 1D PBC observations, suggest that that dE k +1 ,k tends to zero, while dE k, fills up a continuous spectrum spanning at least [0 , D. Barrier Heights
The average uphill/downhill barrier between minima and transition states can be defined [53] as h ∆ i = h E ts i − Σ γ n γ ts E γ min / n ts E γ min is the energy of the minimum γ and n γ ts is the number of transition states connected to thatminimum. The naive uphill/downhill barrier is given by h Λ i = h E ts − E min i . As noted in [53], the average over minima in the second term of h ∆ i is usually weighted towards thelower energy minima, since they are connected to more transition states. This organisation makes h ∆ i larger than h Λ i , as we see in the plots of the average barriers in Figures 6. In these plots, only distinctnon-degenerate lowest-energy rearrangements were considered in the averages. IV. CONCLUSIONS
The potential energy landscape has been examined for the two-dimensional XY model (with nodisorder) with both periodic and anti-periodic boundary conditions. Lattices with up to N = 10 latticesites in each direction, i.e. 100 spins, were considered, focusing on the potential energy distributionof minima and the transition states (saddles of index one) that connect them. As expected [36, 39],the number of stationary points increases roughly exponentially with the number of degrees of freedom.Knowledge of the pathways that connect the local minima enables us to construct the first disconnectivitygraphs for the XY model, and hence visualize the potential energy landscape. These graphs reveal thatthe landscape is funnelled in each case, with a well-defined global minimum, and small downhill barriersconnecting this structure to the higher-energy configurations. Hence all of these 2D XY landscapesbelong to the class of systems identified as good ‘structure seekers’, which includes ‘magic number’atomic and molecular clusters, naturally occurring proteins, and self-assembling mesoscopic systems,including crystals [2, 4, 42, 45]. Minimization from random starting points confirms that the globalminimum is readily located in each case; the funnelled organisation of the landscape is reinforced by theexistence of relatively large basins of attraction for the global minima compared to the higher energyminima, and this effect grows with increasing lattice size.Although the samples of stationary points are not exhaustive for the larger lattice sizes, we candraw some further general conclusions. First, for a given lattice size, N , there are more minima forantiperiodic boundaries than for periodic boundary conditions. Second, as N increases the energy rangespanned by the local minima increases, as one might expect from extensivity of the energy. This effectis also visible in the disconnectivity graphs. However, the probability distribution for the energy of thelocal minima tends to shift towards the global minimum for larger lattice sizes.The trends we have identified have far-reaching implications for the thermodynamics and global ki-netics of the 2D XY model, which we will investigate in future work. Given the wide-ranging applications8f this model, which include superconductivity, superfluidity, liquid crystals, Josephson junctions, andthe fundamental importance of this Hamiltonian in lattice gauge theory [8–10], the energy landscapeperspective may provide new insight into a variety of different research fields.DM was supported by the U.S. Department of Energy under contract no. DE-FG02- 85ER40237and DARPA Young Faculty Award. CH acknowledges support from Science and Technology FacilitiesCouncil and the Cambridge Home and European Scholarship Scheme. MS acknowledges support bythe Research Executive Agency (REA) of the European Union under Grant Agreement PITN-GA-2009-238353 (ITN STRONGnet). DJW gratefully acknowledges support from the EPSRC and the ERC. [1] D. J. Wales. Energy Landscapes : Applications to Clusters, Biomolecules and Glasses (Cambridge MolecularScience) . Cambridge University Press, January 2004.[2] D. J. Wales.
Phil. Trans. Roy. Soc. A , J. Phys. Chem. B , Phil. Trans. Roy. Soc. A , Rev. Mod. Phys. , (1) 167–187 (2008).[6] J. N. Murrell and K. J. Laidler. Trans. Faraday. Soc. , J. Phys. C: Solid State Physics , Phys. Rept. ,
203 (2013).[9] D. Mehta.
Ph.D. Thesis, The Uni. of Adelaide, Australasian Digital Theses Program (2009).[10] D. Mehta and M. Kastner.
Annals Phys. , Rev. Mod. Phys. ,
137 (2005).[12] L. von Smekal, D. Mehta, A. Sternbeck and A. G. Williams.
PoS , LAT2007 (2007).[13] L. von Smekal, A. Jorkowski, D. Mehta and A. Sternbeck.
PoS , CONFINEMENT8 (2008).[14] D. Mehta.
Phys. Rev. E , Adv. High Energy Phys. , Eur. Phys. J. Plus ,
91 (2012).[17] M. Kastner and D. Mehta.
Phys. Rev. Lett. , JHEP ,
018 (2012).[19] D. Mehta, J. D. Hauenstein, and M. Kastner.
Phys. Rev. E , Phys. Rev. E , 052143 (2013).[21] B. Greene, D. Kagan, A. Masoumi, D. Mehta, E. J. Weinberg and X. Xiao. Phys. Rev. D
22] D. Martinez-Pedrera, D. Mehta, M. Rummel and A. Westphal. JHEP , 110 (2013).[23] Y.-H. He, D. Mehta, M. Niemerg, M. Rummel and A. Valeanu. JHEP , 050 (2013).[24] D. Mehta, A. Sternbeck, L. von Smekal, and A. G. Williams.
PoS , QCD-TNT09 025 (2009).[25] C. Hughes, D. Mehta, and J. I. Skullerud. Annals Phys. , 188 (2013).[26] R. Nerattini, M. Kastner, D. Mehta and L. Casetti, Phys. Rev. E,
87 032140
J. Phys. Chem. B , Mathematics of Computation ,
773 (1980).[30] D. Liu and J. Nocedal.
Math. Prog. ,
503 (1989).[31] D. J. Wales.
J. Chem. Soc. Faraday Trans. ,
653 (1992).[32] D. J. Wales.
J. Chem. Soc. Faraday Trans. , Phys. Rev. B , Chem. Phys. Lett. ,
185 (2001).[35] S. A. Trygubenko and D. J. Wales.
J. Chem. Phys. , J. Chem. Phys. , Comp. Phys. Commun. , Proc. AFIPS Spring Joint Comput. Conf. ,
307 (1968).[39] F. H. Stillinger and T. A. Weber.
Science ,
983 (1984).[40] O. M. Becker and M. Karplus.
J. Chem. Phys. , J. Chem. Phys. , , 171101 (2013).[42] D. J. Wales, M. A. Miller, and T. R Walsh. Nature ,
758 (1998).[43] S. V. Krivov and M. Karplus.
J. Chem. Phys. , J. Chem. Phys. , Curr. Op. Struct. Biol. , Curr. Op. Struct. Biol. ,
154 (2008).[47] D. Prada-Gracia, J. Gómez-Gardenes, P. Echenique, and F. Fernando.
PLoS Comput. Biol. , Europhys. Lett.
479 (1993).[49] P. Sibani and P. Schriver,
Phys. Rev. B Ber. Bunsenges. Phys. Chem.
J. Phys. Condensed Matter.
143 (1996).[52] J. C. Schön,
J. Phys. Chem. A
J. Chem. Phys. , ABLES N Table I: The number of minima and saddles of index one located for different lattice sizes N × N , withboth PBC and APBC.11 IGURES -1 PSfrag replacements N N u m b e r o f M i n i m a Figure 1: Number of minima as a function of the number of lattice sites, N , for each dimension. Thestraight lines are the corresponding best fits for the data-points, i.e. the number of distinct minima andthe total number of minima including degeneracies increases roughly exponentially with increasing N .12 Sfrag replacements0.00.10.20.30.40.50.60.70.80.91.0 PSfrag replacements0.00.10.20.30.40.50.60.70.80.91.0PSfrag replacements0.00.10.20.30.40.50.60.70.80.91.0 PSfrag replacements0.00.10.20.30.40.50.60.70.80.91.0PSfrag replacements0.00.10.20.30.40.50.60.70.80.91.0 PSfrag replacements0.00.10.20.30.40.50.60.70.80.91.0
Figure 2: Disconnectivity graphs for the 5 ×
5, 6 ×
6, 7 ×
7, 8 ×
8, 9 × ×
10 APBC lattices.Each of the two insets represents an example minimum for the corresponding N × N lattice. Eacharrow in these insets represents the corresponding value of θ i at the lattice-site i . At each lattice site i ,we compute the local energy P dj =1 (1 − cos( θ i +ˆ µ j − θ i )), which is in the range [0 , Sfrag replacements0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3: Disconnectivity graphs for 5 ×
5, 6 ×
6, 7 ×
7, 8 ×
8, 9 × ×
10 PBC lattices. Each ofthe two insets represents an example minimum for the corresponding N × N lattice. Each arrow inthese insets represents the corresponding value of θ i at the lattice-site i . At each lattice site i , wecompute the local energy P dj =1 (1 − cos( θ i +ˆ µ j − θ i )), which is in the range [0 , Sfrag replacements N E n e r g y D i ff e r e n ce dE , dE , dE , dE , dE , PSfrag replacements N E n e r g y D i ff e r e n ce dE , dE , dE , dE , dE , Figure 4: Sequential energy differences between minima when ranked energetically for APBC (left)and PBC (right) as a function of lattice dimension N . dE k +1 ,k is the energy difference between minima k + 1 and k when arranged in increasing order from k = 0. PSfrag replacements N E n e r g y D i ff e r e n ce dE , dE , dE , dE , dE , PSfrag replacements N E n e r g y D i ff e r e n ce dE , dE , dE , dE , dE , Figure 5: Sequential energy differences between transition states k + 1 and k when rankedenergetically for APBC (left) and PBC (right) as a function of lattice dimension N . dE k +1 ,k is theenergy difference between transition states k + 1 and k when arranged in increasing order from k = 0.15 Sfrag replacements N A v e r ag e B a rr i e r Uphill h ∆ i Uphill h Λ i Downhill h ∆ i Downhill h Λ i PSfrag replacements N A v e r ag e B a rr i e r Uphill h ∆ i Uphill h Λ i Downhill h ∆ i Downhill h Λ i Figure 6: Average value of barrier from h ∆ i and h Λ ii