Entropy Production beyond the Thermodynamic Limit from Single-Molecule Stretching Simulations

Single-molecular systems are a test bed to analyze to what extent thermodynamics applies when the size of the system is drastically reduced. Isometric and isotensional single-molecule stretching experiments and their theoretical interpretations have shown the lack of a thermodynamic limit at those scales and the nonequivalence between their corresponding statistical ensembles. This disparity between thermodynamic results obtained in both experimental protocols can also be observed in entropy production, as previous theoretical results have shown. In this work, we present results from molecular dynamics simulations of stretching of a typical polymer, polyethylene-oxide, where this framework is applied to obtain friction coefficients associated with stretching at the two different statistical ensembles for two different system sizes, from which the entropy production follows. In the smallest system, they are different up to a factor of 2, and for the bigger system, the difference is smaller, as predicted. In this way, we provide numerical evidence that a thermodynamic description is still meaningful for the case of single-molecule stretching.


■ INTRODUCTION
Small systems, unlike those that are in the thermodynamic limit, do not have an extensive internal energy. 1 Because of the small number of particles, they are subjected to large fluctuations. Consequently, it becomes more challenging to obtain relations for average quantities, which are standard in thermodynamics and statistical mechanics of large systems. Gibbs thermodynamics, as we know it from standard texts, 2 ceases to apply for such systems. In view of the numerous and important applications in nanotechnology, for instance, in nanofluidics 3,4 and biology, 5 this situation poses a problem: there is a need to describe energy conversion on the small scale, but a lack of sufficient theoretical understanding. At the most extreme end of the small scale, we are not able to properly describe statistical averages for single molecules. Doubt has thus been raised on the applicability of standard thermodynamic equations to the stretching of single molecules under all conditions. 6 In general, the energy involved in the stretching of a sufficiently small polymer depends on whether one controls the stretching length or the stretching force. The average force for isometric stretching differs from that for isotensional stretching. In the long polymer limit, they are the same, however, which has been verified experimentally, computationally, and theoretically. A very good discussion of this is given by Suzen et al. 7 In an earlier paper, 8 some of us extended Hill's theory for thermodynamics of small systems 1 to time-dependent stretching processes, by deriving expressions for the entropy production for isometric and isotensional stretching. This leads to rate laws with friction coefficients that depended on the control variables. The aim of the present work is to calculate such friction coefficients and the corresponding entropy production using computer simulations and to verify that they depend on the control variables. This is the first example of a dynamic coefficient in molecular stretching.
We investigate the molecular stretching numerically using molecular dynamics simulations. 9 As a model, we have chosen to use a united-atom model of poly-ethylene oxide (PEO), cf. Figure 1, well-documented in the literature. 10 This molecular model has all standard modes of movement under tension, translation, rotation, torsion, and, eventually, the breaking of bonds, and lends itself to a testing of the theoretical description.
In our simulations, the stretching process can be controlled by the environment in two different ways. The endpoints of the hydrocarbon chain can be controlled by either an external force, i.e., f ext is a constant, or by fixing the terminal positions of the molecule, i.e., l is a constant. These isometric and isotensional ways to operate are illustrated in Figure 1a,b. The figures show molecules that are not fully stretched.
Typically, torsional degrees of freedom are associated with lower energies and forces than bending, which in turn is associated with lower energies and forces than bond stretching. We thus expect the response to the environment to change as each of these different modes of elongating the molecule becomes accessible. From the simulation results, we shall find the appropriate dynamic description and relate the molecular properties to the dissipation.
In the thermodynamic limit, the rate laws of the two modes of operation are the same. Here, we present for the first time detailed numerical evidence that there is a difference in the dynamics in the two cases, as predicted from the method of Hill. 8

■ THEORY
The thermodynamic basis for our numerical single-molecule stretching experiments was worked out earlier, 8 when we derived the governing equations for isometric and isotensional experiments on single molecules. In the classical thermodynamic limit, the same set of equations applies to both cases. For small systems, however, there are different sets, as each set depends on how the system is controlled by the environment. 1 An introduction to the general idea of Hill and a more extensive explanation on the structure of nano-thermodynamics can be found in a recent book. 11 In the present work, our system is always just one polymer. The length and therefore the number of monomers and the degrees of freedom vary. A bar will be used above a symbol to denote the average property of an ensemble of systems. We recapitulate the results of earlier 8 to provide a basis for the present step, how the equations can be applied to understand simulations andin a possible next stepexperimental results.
Isometric Experiments. In this experiment, we control the temperature T and the length of the molecule, l. The change in the average internal energy of a system is U̅ , given using the Gibbs equation where S is the system entropy and f̅ is the average internal force working on the terminals, see Figure 1a. The average internal energy can also change by adding heat and work to the system, dU̅ = dQ + f̅ ext dl. The length change is a result of a change in the average external force on the terminals, f̅ ext . By introducing these relations in eq 1, we can identify the entropy change in the surroundings by dS = dQ/T, while the average entropy production per unit of time for the system (one molecule) becomes We now denote the velocity by v ≡ dl/dt and the average change in the force by Δf̅ ≡ f̅ ext − f̅ . The rate law for the isometric case becomes Here, ξ l = ξ l (l) is the friction coefficient specific for the lengthcontrolled case. This is now of primary interest, one of the two coefficients we want to find.
Once we know the friction coefficient, we can compute the entropy production from eq 2, that is, dS/dt = v 2 ξ l (l)/T. The entropy production is proportional to the friction coefficient of the length-controlled case. The entropy production is zero when the external force is balanced by the internal force, f̅ ext = f̅ .
Isotensional Experiments. In isotensional experiments, we control the temperature T and the force of the molecule, f ext . The average internal energy changes as The length of a single molecule is now fluctuating, and l̅ indicates its average. The first law takes the form dU̅ = dQ + f ext dl̅ . By the same reasoning as above, we obtain the entropy production per molecule The controlled change in the force is Δf = f ext − f, resulting in the average stretching velocity v ̅ = dl̅ /dt. The rate law in the force-controlled regime becomes (6) where ξ f = ξ f (f) is the friction coefficient under isotensional conditions, the second target of this study. The entropy production then follows as dS/dt = v ̅ 2 ξ f (f)/T. The entropy production is now proportional to the friction coefficient of the force-controlled case.
In the thermodynamic limit, the two friction coefficients are the same. Away from the limit, this is not the case, as the rate laws depend on the set of the environmental control variables in use.
We shall find below that the stretching simulations of PEO with the smallest molecule under investigation gives a friction coefficient for the case of Figure 1a which is around twice the value of the coefficient for Figure 1b, confirming the prediction The Journal of Physical Chemistry B pubs.acs.org/JPCB Article from the theory that we can expect differences between the two coefficients. Force in the Entropic Regime. Figure 1 illustrates the molecule for relatively small forces, when it is in the entropic regime. In this regime, the molecule behaves similarly to the thermodynamic limit because it has numerous degrees of freedom for movements.
We assume that the molecule to a good approximation can be modeled as a freely jointed chain in the entropic regime with an effective bead length b eff and an effective number of beads N eff , with an unfolded length l unf = N eff b eff . 12 In a system with a solvent, this would correspond to an assumption of theta conditions, that is, the solvent is exactly poor enough to increase the intramolecular forces to perfectly balance out the steric effects. The statistics of the configurations of the system then becomes similar to a random walk, and the radius of gyration, The length b eff is expected to be close to the length of each monomer.
At larger extensions, the forces will first become dominated by unfurling of the torsional degrees of freedom, then the bending, and finally the stretching of the bonds. 13 In these regimes, the force and dynamics typically display nonlinearities.
Helmholz' and Gibbs' Energies. Away from the entropic regime, we expect to be in the small-system regime. In this regime, there is a nontrivial size dependency of properties which is normally extensive. This is due to the fact that fluctuations in the different ensembles are different and lead to different size effects.
For the isometric experiments, there is a fluctuating force for each length. If we let The region for the torsional unfolding is marked with an orange background, and the transition region to the monomer-stretching regime is shown more clearly in the insets. In (d), we see that the entropic region for N = 51 is well-described by a freely jointed chain with N eff = 10 and b eff = 4 Å.
we can compute the Helmholtz energy from That is, the integral along the length axis of the force− elongation curves is shown in Figure 2, giving the area below the curves.
For the isotensional experiments, there is a fluctuating length for each force. If we let the Gibbs energy is given by That is, the integral along the force axis of the force− elongation curves is shown in Figure 2, giving the area above the curves.
In the thermodynamic limit, A and G are related by a Legendre transformation. With Δl = l − l 0 and Δf̅ = f̅ (l) − f̅ (l 0 ), we obtain for sufficiently large systems. 8 Small systems in general deviate from this, and the entropy production in the two ensembles is different. However, eq 12 is still valid when the force is linear in the elongation, like it is in the entropic regime. The nonequivalence between the isometric and isotensional statistical ensembles is the result of the difference between the work done to stretch the molecule, f̅ l and f̅ l, respectively. Considering the nonlinear force−elongation relationship f = al + bl 2 + ..., with a and b two constant parameters, we can easily show that up to a linear order, both works coincide. The nonlinear term, however, breaks down the equality, thus indicating the failure of the thermodynamic limit.
For the entropy production, it is useful to evaluate the from eqs 2 and 5, which is greater than or equal to zero in the second order of l for a specific set of lengths and velocities. From this, one would expect the entropy production for the isometric ensemble to be larger than for the intensional ensemble when the force elongation is nonlinear.

■ MODELS AND METHODS
Although the theory presented above is of general applicability, we choose a specific system for our numerical experiments: a chain of poly-ethylene oxide (PEO) of the form CH 3 where each carbon is grouped with its bonded hydrogen atoms. The PEO monomer consists of one oxygen and two carbons along with their hydrogens. As stated above and illustrated in Figure 1, the endpoints of the chain are controlled by either length (Figure 1a, N, l, T is controlled) or by fixing the endpoints in space (Figure 1b, N, f ext , T is controlled).
The potential energy as a function of the coordinates of the coarse-grained particles has contributions from stretching, bending, and torsion. Using a model that includes these different dynamics allows us to examine the effect of the different modes of stretching and the nonlinearities on the results. The force field is compatible with the LAMMPS 14 simulation package that has been used for all of our computations.
The bond stretching is given using a Morse potential ij ij ij (13) which saturates to a finite value at large separations. The parameters used for the dissociation energies D ij were obtained from density functional computations, 15 and the parameters for α ij were found by requiring the Morse potential to have the same curvature as the harmonic bond, that is, The harmonic force field parameterization is taken from van Zon et al., 16 based on a modification of the explicit atom force field of Neyertz et al. 17 The potentials for the bending and torsion of bonds are and where i, j, k, and l are the atoms joined by consecutive covalent bonds and K ij s , K ijk b , and K ijkl t and r ̅ ij and θ̅ ijk are force constants and reference values, respectively, of stretching (s), bending (b), and torsion (t) energy contributions, selected to reproduce molecular properties measured by spectroscopy or computed by ab initio methods. Note that the sum of the torsional coefficients includes every possible dihedral. Nonbonded interactions were not taken into account, which means that our model polymer is surrounded by an implicit theta solvent. We make this choice because an ideal chain of interacting subunits would deviate from a Gaussian chain even in the thermodynamic limit. 12 The force field parameters we used are presented in Table 1. 10,16,17 The temperature was controlled with a Langevin thermostat, which mimics the viscous aspect of a solvent. During sampling, the relaxation time was set to 1 ps and the temperature was set The Journal of Physical Chemistry B pubs.acs.org/JPCB Article to 300 K. The time step used in the simulations was 1 fs. All quantities presented were averaged over 200 samples. We obtain initial conditions with a low potential energy using a simulated annealing approach. After the initialization setup, all samples are heated up to 2000 K during 0.1 ns before the temperature is slowly decreased during 1 ns.
Case Studies. In the present paper, we present investigations of three different molecule sizes, N = 12, N = 24, and N = 51. Some simulations were also performed with N = 102. The forces varied from 0.01 up to 5 nN or up to the failure limit of the molecule. The length-controlled simulations were sampled evenly in the length, while the force-controlled simulations were sampled evenly on a log scale in the force. This was done to distribute the data points more evenly along the force−elongation curve. To ease the comparisons between system sizes, the molecule length will be presented in units of the longitudinal length divided by the number of bonds l b ≡ l/ (N − 1) and l̅ b ≡ l̅ /(N − 1).

■ RESULTS AND DISCUSSION
To obtain an intuitive understanding of the behavior of the molecule during stretching, it is useful to study the cylindrical radius R c , defined here as the radius of the smallest longitudinal cylinder that can envelop the molecule, shown for N = 24 in Figure 3. There is a sequence of collapses, to be elaborated on below. Four snapshots illustrate the molecular conformation in these regimes. At small lengths, we have a regime dominated by the entropic elasticity, here, the radius R c is 2.3 Å and relatively constant. When the molecule is stretched above l b = 0.5 Å, the torsional degrees of freedom are the first to be confined, and the molecule is unfolded from a helical to a planar configuration. This transition where the C−O−C−C backbone changes from a trans-gauche (ttg) order to an alltrans configuration (ttt) is elaborated in section Torsional Unfolding. This is followed by the unbending and finally the bond stretching. Especially, in regions where several types of dynamics are at play, there is a nonlinear response to stress.
Various Stretching Regimes. In the force−elongation curves shown in Figure 2 for the systems with N = 12 (a), N = 24 (b), and N = 51 (c), we can again identify the different regimes. The entropic regime is shown more clearly for N = 51, see Figure 2d, where lengths below 0.05 Å are considered to be close to zero. The data in this region are consistent with a linear curve. The range where torsion plays a role is indicated by an orange background. The nonlinear transition zone to the monomer-stretching regime is also displayed in more detail in the insets.
Entropic Regime. A predominantly linear relation between force and length develops when 0.05 Å < l b ,l̅ b < 0.47 Å. This is the entropic regime, for which results for N = 51 are enlarged in Figure 2d. From the slope of this curve, we find the effective length b eff of the neighboring units of the ideal chain that gives the correct force−elongation behavior of the molecule in this regime. Within the accuracy of the data presented in Figure 2d, we see that the elongation behavior in this regime is welldescribed by an ideal freely jointed chain for forces up to about 0.05 nN. With a persistence length b eff /2 of 2 Å, 18 we effectively have N eff = 10 beads. The persistence (Kuhn) length b eff corresponds to approximately twice the length of the individual monomers, explained by the bending and torsion, which effectively stiffen the chain. The force-and lengthcontrolled cases appear identical in this regime, as the force− elongation curve here is well-described by a linear function. These findings are in line with eq 12.
Torsional Unfolding. As the molecule is stretched further, the degrees of freedom are reduced, and the freely jointed chain model is no longer applicable. The torsional degrees of freedom are the first to be confined, and this occurs in the region 0.47 Å < l b ,l̅ b < 1.1 Å, marked with an orange background in Figures 3 and 2. The beginning of the interval was found by looking at the deviation from linearity in Figure  2d, and the end of the interval was found from the inflection point of Figure 3. PEO strands are known to attain a helical shape in the crystalline state, in which the bonds of the C−O− C−C backbone are folded in a trans-gauche (ttg) order. 19 This can be seen in the first two snapshots in Figure 3 and is also the case for PEO dissolved in water. 20 An increase in the force gives rise to a transition from a helical ttg order to an elongated, planar all-trans configuration (ttt), as seen in the last two snapshots in Figure 3.
From Figure 2a−c, we can see a systematic deviation that varies with molecular size. This is emphasized in the insets. For The Journal of Physical Chemistry B pubs.acs.org/JPCB Article N = 12, we observe pronounced oscillations in the force− elongation curve; for N = 24, we observe smaller oscillations; and for N = 51, we observe no oscillation. These oscillations in the length-controlled ensemble are finite size effects that originates from local maxima in the potential of mean force associated with the unfolding of the molecule. Here, the molecule is mechanically unstable, and these modes are not accessible in the force-controlled ensemble. 5 This leads to different fluctuations in the two ensembles. Monomer-Stretching Regime. As the molecule is extended above l̅ b > 1.1 Å, the individual monomers are elongated. The molecule is unbending, and the potentials for the stretching, bending, and torsion give rise to a molecule-specific segment elasticity, 13 increasingly dominated by the stretching of the covalent bonds.
In this region, a small systematic difference appears in the force−elongation curves between the length-controlled and the force-controlled stretching experiments. This can be seen in the inset of Figure 2a These nonlinearities from the stretching of the Morse potentials give rise to different fluctuations in the two ensembles, and we expect to see an effect of the small system size. The differences between the force−elongation curves shown in Figure 2a−c are the largest in the transition regime to the monomer-stretching regime, emphasized in the insets. The differences are small but they are finite and systematic.
Gibbs and Helmholtz Energies. The free-energy differences, and the deviation from the Legendre transform in eq 12, are computed from the force−elongation curves shown in Figure 2a−c, according to section Helmholz' and Gibbs' Energies, and shown in Figure 5. We divide by the work required to stretch the molecule completely, in order to compare the different system sizes. The largest free-energy difference is observed in the transition from the torsionalunfolding regime to the monomer-stretching regime, see the insets of the force−elongation curves in Figure 2a−c. Both in the case of N = 12 and N = 24, there is a clear correspondence between the deviations in the force−elongation curves in this region and the peak in the free-energy difference, as shown in Figure 5. There is a significant deviation from eq 12, with the smallest system showing the largest deviation, as expected.
Friction Laws. Force-Controlled Simulations. We can now use our simulations to estimate the friction coefficient ξ f = ξ f (f) in eq 6. This was done for the systems with N = 24 and N = 51 by perturbing the force and determining the rate of change in the average length. To this end, we first generated 200 independent samples, each equilibrated at 150 different constant forces f 0 for 5 ns. At time t = 0, the force on each of these samples was increased by 140 different force increments in the range 4−28%. The length as a function of time before and after the increase in the force is shown in Figure 6 for three force increments in the system with N = 51, averaged over 200 samples.
From these results, we find that the time scale for the initial linear force response is 0.5 ps for N = 51. As one can see in Figure 6, this does not appear to depend on the magnitude of the force increment. The ratio of the force increment to the increase in the linear response is equal within the accuracy of the data points. A similar investigation of N = 24 results in a time scale of ∼0.2 ps. The time scale for the linear regime is related to the relaxation time of the system, which depends on the length of the molecule. Other time scales in the range 0.1− 1 ps was explored and was found to give similar results, although with increased fluctuations, indicating a reasonably good robustness on this parameter. Continuing with the chosen time scales, the linear response dl̅ /dt was then  what is the case in the thermodynamic limit, the friction coefficient was largely dependent on the value of the force and the length of the polymer. Length-Controlled Simulations. To estimate the friction coefficient ξ l = ξ l (l) in eq 3 for N = 24 and N = 51, we stretch the molecule in a range of velocities and estimate the increase in the force Δf̅ associated with each stretching velocity for each sample. A total of 200 independent samples were first equilibrated at 150 different constant lengths l 0 for 5 ns, and at time t = 0, the samples were stretched at 80 different constant velocities v = dl/dt in the range 20−100 m/s for 1 ps. The force response from the molecule Δf̅ for each stretching velocity was then averaged over the same time scale as used for estimating the linear response in the force-controlled simulations. The resulting force−velocity curves for molecules N = 51 with initial lengths of l b = 0.824 Å and l b = 1.192 Å can be seen in Figure 7. Again, we found the friction coefficient ξ l = ξ l (l) using eq 3 from the slope of these force−velocity curves. The variation in the coefficient with the length of the molecule or the force applied was similar to the results from the isotensional experiments, but the coefficients for forcecontrolled systems were systematically smaller than those for the length-controlled systems. As the fluctuations increased significantly for shorter lengths, only lengths per bonds larger than 0.4 Å are shown. Both curves showed a maximum near the relative length 1.2 Å per bond, where the Morse potential for bond stretching is strongly nonlinear.
The difference in the friction coefficient can be expected from a dynamical investigation of the system, by considering the time scales and following the approach of Just et al. 21 to obtain the general form of the effective slow dynamics. The length of the molecule acts as the slow variable, and the probability distributions of the fast variables of the internal degrees of freedom of the molecule are different for fixed force and fixed length. This also leads to two different damping constants.
Entropy Production. The force-controlled friction coefficient ξ f = ξ f (f) = ξ f (f(l̅ )) found in the section Force-Controlled Simulations and the length-controlled friction coefficient ξ l = ξ l (l) found in the section Length-Controlled Simulations are presented as a function of the length in Figure  8 for molecules N = 24 and N = 51. The difference between ξ f and ξ l is smaller for the largest molecule, as expected from eq 12.
The entropy production is found by multiplying this coefficient with the constant velocity squared over the temperature. The energy dissipation producing heat in the surroundings is the entropy production times the (constant) temperature. Apart from this trivial rescaling factor, the basic properties are considered to be temperature-independent under the assumption of theta conditions. For very short lengths, the entropy production by definition should go to zero. Although the uncertainty in this region is rather high, we emphasize that zero is within the margin of error. In the region of torsional unfolding, the ensemble difference is the largest for the smaller system with N = 24 compared to the bigger system with N = 51. This is as expected from the discussion of the different stretching regimes. The entropy production reaches a maximum around 1.2 Å per bond for both system sizes, well into the monomer-stretching regime. Again, the ensemble difference is significantly larger for the smallest system. This can be explained by the nonlinearity of the Morse potential for the bond stretching, giving rise to   Figure  4, we see that the maxima appear to coincide. Moreover, any coupling to low-frequency tangential phonons can also very quickly dissipate energy in this regime. We have seen above that the magnitude of the friction coefficient differs between the two stretching modes, with the length-controlled process having a higher friction coefficient than the force-controlled process. It follows that the first process dissipates more energy regardless of the length of the molecule, as expected. Note that the force-controlled simulations significantly display larger size dependence than what is seen in the length-controlled simulations.

■ CONCLUSIONS AND PERSPECTIVES
In small-scale systems, away from the thermodynamic limit, standard thermodynamics is no longer valid. In this case, thermodynamic potentials become nonextensive and statistical ensembles are not equivalent. Even if the system is very small, extensivity can be restored, if one considers the set of replicas of the original system as a large-scale system. Such a procedure, proposed by Hill, 1 makes it possible to apply the method of thermodynamics on very small scales. This method, initially proposed when the system is in equilibrium, was extended 8 to nonequilibrium situations for the case of the stretching of a polymer.
In this article, we have shown that dissipation generated at small scales is sensitive to the lack of equivalence between statistical ensembles at small scales. Based on earlier work, 8 we have carried out simulations well beyond the thermodynamic limit. We have simulated the stretching of a single PEO molecule of length N = 12, 24, and 51 under force-controlled and length-controlled ensembles and have extracted friction coefficients for the largest two systems.
We have confirmed systematic finite size effects in the two ensembles of general nature. In the static case, the finite size effects are most pronounced in the region of torsional unfolding and originate in local maxima in the potential of mean force that are accessible only in the length-controlled ensemble. This is visible for N = 24 and even more for N = 12. In the dynamic case, the finite size effect originates in the two ensembles having different fluctuations. This is predicted by theory and confirmed for the first time for the dynamical coefficient. For short polymers with N = 24, the friction coefficient of isometric stretching is roughly twice the value of that of an ensemble with isotensional stretching. The difference between the friction coefficients decreases when the length of the polymer is increased to N = 51.
Our study shows how nonequilibrium properties are affected by the absence of the thermodynamic limit. The method presented could be applied systematically to the study of irreversible processes that take place at small scales. Figure 8. Estimated friction coefficient from the length-controlled and force-controlled simulations for N = 24 and N = 51. We see that the ensemble deviation is most noticeable in the monomer-stretching regime and that the ensemble deviation is more pronounced for the smallest system.
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article