The impact of bonded interactions on the ground-state geometries of a small flexible polymer

Bonded interactions in coarse-grained models of elastic polymers are commonly represented by the finitely extensible nonlinear elastic (FENE) potential. In this study, we perform parallel multicanonical Monte Carlo simulations to examine the impact of an additional Lennard-Jones term in the bonded potential on the geometry of ground-state structures of a short polymer. Employing microcanonical inflection point analysis and conformational analysis, we construct a hyper-phase diagram and identify ground-state structures with two distinct geometries.


Introduction
The physical and chemical functionality of biopolymers and synthetic variants is intimately related to their conformational structure. Therefore, the systematic study of their dynamic and structural properties is of a fundamental importance and a major topic of current interdisciplinary research. Theoretical approaches are generally not adequate due to the prohibitive complexity of even the simplest molecular systems [1]. On the other hand, computational methods such as the replica-exchange [2], multicanonical [3,4], and Wang-Landau sampling [5], have proven to be most effective. Despite the vast increase in the availability of computational resources over the past decade, simulations of fully atomistic models remain a major challenge. Hence a clever choice of suitable coarse-grained models with a reduced number of parameters is often a prerequisite to a successful computational study. In this paper, we consider a generic model of the flexible elastic homopolymer and examine the effects of adding an extra Lennard-Jones term to the bonded potential. The unexpected consequences for the ground-state structures highlight the importance of a careful choice of model parameters.

Model and Methods
We employ a generic model of a flexible elastic homopolymer to study the effects of the shape and width of the bonded potential on the formation of low-energy structures for a polymer chain with N = 15 monomers. Finite-size effects are essential and the 15mer shows a particularly distinct and clear transition behavior, which makes it a perfect candidate for this systematic study. Non-bonded interactions are represented by the truncated and shifted Lennard-Jones (LJ) potential where The energy scale is set to ǫ = 1 and the van-der-Waals radius to σ = r 0 /2 1/6 , where r 0 = 1.0 is the location of the potential minimum. We select a cut-off radius at r c = 2.5 σ and introduce a shift U LJ (r c ) ≈ −0.0163169 ǫ to avoid discontinuities in the potential. In our model, monomers adjacent in the linear chain interact via a modified bonded potential where R = 3/7 and K = 98/5. In addition to the standard anharmonic FENE (finitely extensible nonlinear elastic) potential [6,7,8], the modified potential contains an additional Lennard-Jones term adjusted by a bond flexibility control parameter η. The bonded potential is shifted by −(ǫ + U shift ) in order to match the minimum energy with the non-bonded interactions. The maximum bond extension is limited by the FENE potential, which diverges as r → r 0 ±R. Increasing the value of η introduces asymmetry to the bonded potential and the energy cost associated with non-optimal bond lengths is increased. In particular, compressed bonds result in high energy penalties as η becomes large. The effects of different η values on the shape and width of the bonded potential is exhibited in Fig. 1. The total energy of a configuration X = ( r 1 , ..., r N ) with the monomer-monomer distances r ij = | r i − r j | is given by In order to enhance the sampling of low-energy configurations, a parallel version of multicanonical sampling [3,4,9] is employed in this simulation. The standard multicanonical runs are performed in K replicas independently with the same initial weight functions but different random seeds. Displacement updates are proposed within a cubic box of edge lengths d = 0.3r 0 and accepted according to the probability where W (E(X)) represents the weight function of a given configuration X. After the ith iteration, since the weights are identical in each thread, the energy histograms obtained for each replica can simply be summed up: The total histograms are combined with the current weights to calculate the weights for the subsequent iteration by utilizing the error-weighted recursive scheme [1,3,4].
To construct the hyper-phase diagram, we use the generalized microcanonical inflectionpoint analysis [1,10]. This approach has the advantage of uniquely and systematically locating the transition energies and thus is commonly employed to study pseudophase transitions in finite-size systems. By applying the principle of minimal sensitivity [11] to the derivatives of microcanonical entropy S(E), the (2n+1)th-order transition (n is a positive integer) is identified from the least sensitive inflection point of the 2nth-derivative of entropy and the positive valley in the (2n + 1)th-derivative curve. For a 2nth-order transition, the least sensitive inflection point in the (2n − 1)th-derivative of entropy together with the negative peak in the 2nth-order derivative curve are utilized to locate the transition energy. The specialty of the first-order-like transition is that it can be distinguished from the back-bending region in the inverse temperature β(E) ≡ T −1 (E) = dS(E)/dE and the corresponding positive peak in the γ(E) = dβ(E)/dE curve.

Canonical and microcanonical analysis
First we discuss the results of canonical statistical analysis applied to our generic model. Heatcapacity curves as functions of temperature ( Fig. 2(a)) are constructed using the time series of data collected in the multicanonical production run. At T ≈ 0.34, broad prominent peaks, indicating the freezing transition, are identified for all simulated η values. At this transition, globular structures change to more compact crystalline or amorphous structures. For η = 0, an additional peak emerges at T = 0.11, suggesting the existence of a solid-solid transition.
With increasing values of η the peak becomes more pronounced as it gradually shifts towards zero temperature and finally disappears when η ≥ 0.1. However, the order of the individual transitions remains ambiguous and the broad peaks of the freezing transition may envelope several transitions which cannot be resolved by the methods of canonical analysis. Therefore we further examine the system using a more systematic and robust approach. One such method that has proven to reliably signal transitions is the microcanonical inflection-point analysis [10].
The microcanonical results are shown in Fig. 2(b-d) for six different values of η. Careful inspection of the first and second derivatives of β(E) in the energy region E ∈ [−45.5, −33] reveals that the broad peak in the canonical specific heat encloses two distinct transition signals; clear indication that the freezing transition is a two-step process. The first signal located at E ≈ −44, is a fourth-order transition indicated by the corresponding least sensitive inflection point in δ(E). The second transition, found at E ≈ −38, is of third order for η ≤ 0.2, but it is classified as a second-order transition for higher values of the bond flexibility control parameter. In agreement with the canonical results, we have also identified signals corresponding to a solid- polymer expands in free space and forms random-coil structures. As the energy decreases, the expanded chain undergoes a second-order collapse transition and enters the "liquid" pseudophase consisting mainly of globular structures. Passing the third/second-order transition associated with the nucleation process, the polymer enters the S ic pseudophase in which incomplete icosahedral structures become dominant. Further decrease in energy weakens thermal fluctuations and allows for the formation of a stable surface layer. The transition associated with the surface formation process is of fourth order. Visual inspection of low-energy structures reveals that icosahedral geometries are dominant. However, for η ≤ 0.1, the additional solidsolid transition suggests the existence of low-energy conformations with unexpected geometric properties. In order to examine the low-energy structures systematically, we carry out a careful structural analysis utilizing a suitable set of order parameters.

Structural analysis
Various order parameters, such as the number of monomer-monomer contacts, radius of gyration, or radial and angular distributions, have proven to provide valuable insight into the thermodynamic and structural properties of polymer systems. Based on the microcanonical results in Fig. 3, we expect the existence of two solid phases when the strength of the bond flexibility parameter is sufficiently small (η ≤ 0.1). We aim to identify the dominant structures in the low-energy phases and to gather additional data supporting the existence of the solid-solid transition line. For this purpose, we employ a set of order parameters exploiting the symmetry  properties of real spherical harmonics [12]. We define a polymer core to consist of monomers within a distance r core < 1.25σ of the central monomer, which has been chosen to be nearest to the center of mass. Let C = { r 1 , ....., r M } be the coordinates of a core with M monomers. Various core geometries can be distinguished using the set of rotationally invariant order parameters where is the average of the real spherical harmonics evaluated at the locations of the core monomers. The connection between the real and complex spherical harmonics is given by Using of the order of 10 6 polymer structures per value of η, we computed Q l up to l = 6 and found that Q 6 can be used most effectively to resolve the geometries of the low-energy conformations. We present the results in the form of intensity plots in Fig. 4. The probability of detecting a structure with a specific value of the order parameter at an energy E is represented by shading; red indicating the maximum probability and black corresponding to zero. In agreement with the microcanonical and canonical results, we detect a single solid phase for η > 0.1, corresponding to  Two distinct lowenergy structures of the elastic 15-mer.
(a) Compact structure with a stable icosahedral core and two monomers displaced onto the incomplete second layer. (b) The bihexagon is the preferred groundstate geometry for η ≤ 0.1.
The shape of the bonded potential has undoubtedly a strong effect on the geometry of the ground state. Having identified the two dominant structure types, we may ask why the additional LJ term in the bonded potential eventually precludes the existence of the bihexagonal phase. The answer is readily obtained by comparing the average bond lengths for the icosahedral and bihexagonal structures. The bihexagon accommodates all monomers into a single shell allowing for a larger number of non-bonded interactions and consequently lower energy. However, the two six-monomer rings of the bihexagon contain significantly compressed bonds (r bond ≈ 0.88r 0 ), which become energetically infeasible as η increases. In contrast, we find near-optimal bond lengths in the icosahedron (r bond ≈ r 0 ), hence the "narrowing" of the bonded potential imposes no additional energetic penalty.

Summary
By means of multicanonical simulations of a generic model for elastic, flexible polymers, we have investigated the structural behavior of a 15-mer upon changing a model parameter η that controls the shape of the bonded potential. For small values of this parameter, a freezing transition into an icosahedral phase precedes a solid-solid transition into low-energy states with bihexagonal geometry. The non-optimal bond lengths found in bihexagonal conformations cause a large energy penalty due to the "narrowing" of the bonded potential if η is increased. Hence only a single solid phase remains for η > 0.1, which is icosahedral. The striking consequences of a relatively small modification to the standard model of elastic, flexible homopolymers illustrate the importance of a careful choice of model parameters.