Monte Carlo simulations of a polymer chain conformation. The effectiveness of local moves algorithms and estimation of entropy

A linear chain on a simple cubic lattice was simulated by the Metropolis Monte Carlo method using a combination of local and non-local chain modifications. Kink-jump, crankshaft, reptation and end-segment moves were used for local changes of the chain conformation, while for non-local chain rearrangements the "cut-and-paste" algorithm was employed. The statistics of local micromodifications was examined. An approximate method for estimating the conformational entropy of a polymer chain, based on the efficiency of the kink-jump motion respecting chain continuity and excluded volume constraints, was proposed. The method was tested by calculating the conformational entropy of the undisturbed chain, the chain under tension and in different solvent conditions (athermal, theta and poor) and also of the chain confined in a slit. The results of these test calculations are qualitatively consistent with expectations. Moreover, the obtained values of the conformational entropy of self avoiding chain with ends fixed over different separations, agree very well with the available literature data. Figure Visualization of the neighborhood of two local chain microconformations containing a bead (indicated by the arrow) which a) can and b) cannot be moved by the kink-jump. In red there are marked the possible trajectories of chain fragment which can block the adjacent site


Introduction
The knowledge of the Helmholtz free energy A (or Gibbs free energy G) of a polymer is essential for understanding many complex processes involving macromolecules, as for instance, the adsorption or grafting of polymer chains on interfaces [1][2][3][4][5][6], the formation of polymer bridges [7,8] or protective layers in processes of flocculation or steric stabilization of colloidal suspensions [8,9], the formation of secondary and tertiary structure of proteins [10][11][12], the release of the genetic material from the viral capsids [13][14][15][16] and many others. However, in such complex processes as those mentioned above, it is not possible to analytically compute the changes in the free energy and, particularly, in its component-the conformational entropy. This problem can be overcome by the use of simulation methods which provide the ability to study complex systems in great detail [17][18][19]. The application of these methods requires a construction of a simulation model which accurately represents the processes being modeled. The polymer model appropriate for investigation of the above-mentioned phenomena should include the excluded volume and energetic interactions. It should also account for deformation of a polymer molecule from its equilibrium conformation in the presence of external constraints (e.g., a macromolecule confined in a nanopore or attached to an impenetrable surface) or/and forces (e.g., a macromolecule bridging two particles undergoing Brownian motion). These requirements impose certain restrictions on the choice of the method and algorithm for the estimation of free energy and conformational entropy of a polymer molecule. Aiming at the maximum possible simplicity of the model and taking into account the features it should include, a linear homopolymer chain simulated on a lattice by the Metropolis Monte Carlo (MMC) method [20,21] was chosen for the study. Constant length of the chain and fixed positions of terminals were assumed. The latter assumption arises from the adopted strategy for the examination of relationship between the Helmholtz free energy of a polymer molecule and the deformation extent, based on simulation of conformation of a chain having arbitrarily imposed fixed distance between its ends.
In this study an attempt was made to develop a method permitting estimation of the Helmholtz free energy of a polymer chain on the basis of probabilities of micromodifications of chain conformation in the MCC simulation. The adopted elementary moves of the chain followed the generalized Verdier-Stockmayer algorithm [21][22][23]. In this algorithm a number of local moves can be performed depending on a local chain conformation. The MC simulations using Verdier-Stockmayer algorithm were shown to reproduce some real kinetic properties of macromolecules [24,25]. It was proven, however, that all local chain length-conserving algorithms are nonergodic [26,27]. Although, it is generally assumed that the systematic errors following from nonergodicity are insignificant when compared with statistical ones [28][29][30][31], particularly for relatively short (up to 100 monomers) self-avoiding walks (SAWs) [32]. However, because the model is intended for use in systems in which the nonergodicity error may be more pronounced due to the stretching of polymer chain and/or the presence of space constraints, a hybrid algorithm combining the local and non-local moves was employed in simulations. The use of non-local algorithm allows elimination of the nonergodicity problem [27,32,33]. In this study, the "cutand-paste" method [27] was chosen for non-local modifications of the chain. The use of the bond fluctuation algorithms [28] was abandoned because their essential drawback revealed against the model demands. This drawback is ambiguous definition of the monomer neighborhood, which implies problems when the effect of the energy of monomer-monomer or monomer-solvent molecule interactions on the chain conformation has to be considered.
The aim of the study was to examine the statistics of local micromodifications in the simulation of chain conformation by means of the MMC method. The simulations were performed both for free undisturbed chains and for those whose two terminals were fixed over a selected distance. Three types of solvent conditions (i.e., athermal, theta and poor) were considered. Also, the chain in the presence of geometric constraints (the chain in a slit formed by two parallel impenetrable surfaces) was studied. As a probe of a local chain conformation and a local environment the kink-jump algorithm was used. The calculated values of the Helmholtz free energy of a SAW chain were confronted with the results available in literature [34], obtained by the expanded ensemble Monte Carlo method.

The model and simulation method
The majority of results reported in this paper concern the polymer chain represented by SAW embedded to a simple cubic lattice of the lattice constant b. One monomer is identified with the help of a site on the lattice (bead), so the chain is an ensemble of N consecutive occupied sites. Mutual positions of beads are limited to the vector set being a permutation of (0, 0, ±b), which corresponds to the lattice coordination number ω equal to 6. The simulations were performed for two different bead numbers N, i.e., N=50 and 100 (the chain of N beads is equivalent to M=N-1 SAW steps, each of the length b) in a box with the edge length considerably greater than the contour length of the longer chain (601b). Periodic boundaries of the simulation box were applied.
The simulations were completed using the Metropolis Monte Carlo (MMC) method [27]. In this study we performed two kinds of simulations: 1) for a free unperturbed chain which was employed as a reference system and 2) for chains of arbitrarily chosen distance L H between terminal beads of fixed coordinates. In the first case the initial chain conformation was generated by the static MC method [27] and next it was randomly modified using algorithms of the following local micromodifications: end-bead move (E), kink-jump (K), crankshaft (C) and reptation (R) motions [28] and two non-local cut-and-paste modifications: the inversion of a chain piece (I) and the reflection of a chain piece through a perpendicular bisector of a straight line joining the first and last bead of this subchain (F) followed by motion (I) (the modification F is necessary for ergodicity [27,33]). In the second case the simulation started from a regular spiral-like conformation which was next equilibrated via four types of modifications: K, C, I and F.
The MC step, defined as the number of modifications needed to give each of the beads the possibility to move once, contained N local micromodifications chosen at random, thus the frequencies of K and C motions were on average the same, whereas those of E or R motions (if applied) were N/2 times lower. The probability of performing the non-local cut-and paste modification, involving simultaneous changes in coordinates of N/3 beads on average, was 10 3 times smaller than that of local modifications. Thus, the ratio between frequencies of attempted local "internal" moves (i.e., K or C motions), micro-moves of end beads (i.e., E or R motions) and non-local cut-endpaste moves per a bead was 1:0.02:0.03, respectively. Because the increase in the frequency of cut-and-paste modifications up to 10 times involved no changes in the results obtained (i.e., the differences between the results were within the statistical errors of the simulation), it was assumed that the frequency applied was sufficient to ensure the ergodicity of the procedure. This frequency of the non-local modifications allowed efficient simulations. On the other hand, the relatively high frequency of local modifications ensured the statistical significance of the probability of K motions, which were used as a probe of chain conformation.
In addition to the simulations performed for pure SAW chains where all molecular interaction energies are assumed as equal to 0, the conformations of chains immersed in the theta solvent (ε PS =ε SS =0, ε PP /k B T=−0.2693 [35]) and in poor solvent (ε PS =ε SS =0, ε PP /k B T=−1) were also simulated. Symbols ε PS , ε SS and ε PP denote the energies of beadsolvent molecule, solvent moleculesolvent molecule interactions and the energy assigned to each non-consecutive pair of beads lying on neighboring lattice sites, respectively. In the course of these simulations, new conformations were accepted or rejected by the Metropolis procedure in the standard way [20,28].
Some simulations were performed for the chain with terminal beads attached to two parallel surfaces forming a slit of the width equal to L H +2b. The surfaces were rigid and impenetrable to the polymer. They were purely repulsive and there was no chain adsorption (athermal surface) except for the irreversible attachment of the chain ends. Thus, the influence of surfaces on the polymer conformation was solely of entropic character.
All the results presented in the paper are averaged over 10 4 conformations. Each conformation is a result of 10 6 MC steps. The first half of motions (i.e., 0.5⋅10 6 MC steps) were used for the initial chain thermalization. For the second half the following parameters characterizing the chain conformation were computed: the probability of performing a given move without violation of the continuity of the chain (i.e., the probability that the skeletal constraint is satisfied) and the probability of move execution without violation of the excluded volume constraint.
The probabilities discussed are described by the symbol P X (Y) Z , where X is the movement limitation (skeletal or excluded volume denoted as S and E, respectively), Y is the type of micromodification (E, K, C or R) and Z (if needed) indicates the state of the chain (Z=0 for the unperturbed chain).
For comparative purposes, the simulations for the conventional random walk were also performed.

Probabilities of moves
In this paragraph we present the results of calculations of probabilities of local micromodifications considered for chains with both free and fixed ends, obtained on the basis of simulation of chain conformations in athermal conditions. Table 1 summarizes the values of probabilities predicted theoretically and obtained from simulations in the case of free chain. Details of considerations on which the theoretical predictions were based are described below. The probabilities of moves satisfying the chain continuity constraint and of those satisfying the excluded volume constraint are considered separately. The results obtained for the chain with fixed distance between ends are only briefly presented here; their detailed discussion is given in further parts of the paper.

The skeleton effect
The obtained values of P S (R) 0 and P S (E) 0 equal to 1 are easy to predict since by no means the move of end bead of the chain can violate the skeleton constraint.
The prediction of P S (K) 0 and P S (C) 0 is somewhat more complex. The probability P S (K) 0 can be evaluated as follows. First, note that only corner located beads can perform kink-jumps, as schematically elucidated in Fig. 1. The probability that three succeeding beads (indexed as i-1, i, i+1) form the corner is given by the ratio (ω-2)/(ω-1), as can be easily deduced from Fig. 2a. The number of possible corners with vertices on bead i can be additionally reduced by the occupancy of the vertex adjacent sites by other beads. The probability of occupancy of a single such site with a given bead decreases rapidly with its distance from bead i (measured along the chain) and thus, for instance (see Fig. 2b), for beads with indices i-3, i-5…, these probabilities are: 1/(ω-1) 2 , 5/(ω-1) 4 and the like. Hence, the value of P S (K) 0 can be found from the following equation: as equal to 0.763 which is close to the value determined from the simulation results (see Table 1). The calculation of the probability P S (C) 0 , which is equivalent to the probability of occurrence of a local U-shaped conformation, is based on similar considerations as those presented above. At first, let us notice that the formation of such a conformation can be thought Table 1 The probabilities of different types of motions estimated theoretically (theor) and calculated on basis of the simulation (sim) performed for unperturbed chains of N=100 and N=50 in the athermal conditions Quantity Micromodification as performed in two steps (see Fig. 3a): 1) the formation of a corner by beads i-1, i and i+1 (the same as in the case of kink-jump movement), 2) the addition of a third edge by suitable location of (i+2)th bead. Then, the probability of formation of this microconformation is equal to P S (K) 0 ⋅1/(ω-1). However, the second step can be hindered by the occupancy of the "suitable" lattice site by other beads. As it can be deduced from Fig. 3b, the probabilities that this site is blocked by beads with indices i-2, i-4… are equal to 1/(ω-1), 3/(ω-1) 2 … Thus finally, the probability that the randomly chosen part of the chain can perform the crankshaft move satisfying the skeleton constraint reads: The P S (C) 0 calculated from Eq. 2 is equal to 0.103 and roughly agrees with that obtained from the simulation (see Table 1).
The probabilities of motions not violating the chain continuity constraint, calculated directly from the simulation results for the chain with fixed positions of terminal beads, are presented in Fig. 4. As expected, both probabilities tend to zero when L H reaches the contour length of the chain.

The excluded volume effect
The probabilities of micromodifications not violating the excluded volume constraint are also not difficult to predict in the case of end-segment and reptation moves. The former move, if permitted, changes the orientation of a selected end bead by 90°or 180°, whereas the latter move cuts-off (with the probability equal to 1) one bead from either end of the chain and reattaches it to the opposite chain end. Thus, one can consider these probabilities to be measures of the number of vacant lattice sites adjacent to the selected end bead. One can also notice, that the product P E (R) 0 (ω-1) is equivalent to the effective coordination number of the lattice which was used for calculation of the conformational entropy by means of the statistical counting method (SCM) based on the static MC procedure [36,37]. Assuming that the value of ω eff may be substituted by the average value of the infinitely long chain (for N→∞ ω eff =4.6839 [38]) the probability P E (R) 0 equal to 0.937 is obtained, which is similar to that found from the simulation for the chain of N=100.
Prediction of the P E (K) 0 probability is more complex. As results from Fig. 5, this probability depends on the bead location in the chain and it assumes the highest values (close to P E (R) 0 and P E (E) 0 ) at the chain ends, rapidly decreasing Fig. 1 The kink-jump motion trials for two different locations of a bead in the chain a the forbidden move, b the permitted move Fig. 2 Illustration of estimation of P S (K) 0 : a possible conformations of a chain fragment composed of three beads with indices i-1, i and i+ 1, b a possible location of (i-3)th bead which reduces the number of corner forming conformations (in redthe chain fragment whose formation is disallowed) Fig. 3 Illustration of the estimation of P S (C) 0 : a a possible conformation of a chain fragment composed of four beads with indices i-1, i, i+ 1 and i+2, b an exemplary possible (i-4)th bead location disallowing the completion of U-shaped conformation by suitable orientation of (i+ 2)th bead (in redthe chain fragment whose formation is disallowed) Fig. 4 The influence of L H value on probabilities P S (K) and P S (C). The horizontal lines indicate corresponding probabilities determined for the free chain (SAW, N=100) to a certain almost constant value when the bead distance from the ends (measured along the chain) increases. Such a character of this dependence is consistent with intuitive expectation, and thus allows a conclusion that P E (K) 0 provides information about the likelihood of vacancy of a site adjacent to ith single bead of the chain.
In order to find a correlation between the probability P E (K) 0 and the number of vacant sites around the chain let us analyze the occupancy of the site adjacent to a single ith bead. The probability that the site is not empty can be considered as the probability of the union of two independent events: the occupation of the site by a bead with index lower than i and the site occupancy by a bead with index larger than i. Assuming that these two events can occur with the same probability X, the site occupation probability is then equal to 2X-X 2 . Now, let us consider the probability Y of occupancy of the goal site of a kink-jump move. Since this site has two neighbors (as shown in Fig. 6a), the probability of its occupancy is again equivalent to the probability of the union of two independent events. Thus, on the basis of the above reasoning, we can write: Now in turn, let us notice that the number n 2 of possible chain trajectories passing through a site neighboring two beads (where the kink-jump move is permitted) is lower than the respective number n 1 for a site neighboring only one bead (where the kink-jump is disallowed). This means that the probability Z of occupancy of a site adjacent to the chain disallowed for the kink-jump move is higher than the probability Y by a factor w: where w=n 1 /n 2 . In order to determine the value of w one should count all possible chain trajectories passing through single lattice sites neighboring one and two beads, respectively. Here, for the sake of simplicity, let us consider two step trajectories only. In such a case we obtain n 2 =20, n 1 =25 and hence w=1.25, as explained in Fig. 6.
(For lattices with other coordination number w=(ω-1)/(ω-2)). Finally, the probability P E (K) 0 that the kink-jump move does not violate the excluded volume restriction, i.e., that the target site is empty, is equal to 1-Z. We should emphasize here, however, that the estimate of P E (K) 0 based on considerations presented above is very rough and valid only for free SAW chains, for which one can assume the long distance excluded volume to be small. Finally, combining Eqs. 4 and 5 and assuming that X=1-P E (R) 0 , we get the following relation: From Eq. 6 we obtain P E (K) 0 equal to 0.734 if the value of P E (R) 0 calculated from simulation results is used or P E (K) 0 =0.714 if the value of P E (R) 0 determined on the basis of Eq. 3 is employed. Taking into account that these results are only rough estimates their agreement with the value of P E (K) 0 =0.695 obtained from simulation (Table 1) seems to be satisfactory.
Having evaluated the probability of vacancy of a lattice site neighboring one bead (assumed to be equal to 0.714), one can find the probability P E (C) 0 that two such adjacent sites are empty (for the crankshaft move to be permitted) as equal to (P E (K) 0 ) 2 , obtaining the value 0.510, which is close to that determined from the simulation.
The results obtained for the chain with fixed ends are presented in Fig. 7. The values of both P E (K) and P E (C) increase with increasing end-to-end distance reaching the value of 1 for the fully stretched chain (i.e., L H →bN ⇒ P E (K)→1).  5 The dependence of the relative probability that a bead adjacent site is empty on the position of the bead along the chain (SAW model, N=100, results obtained from a single simulation) Fig. 6 Visualization of the neighborhood of two local microconformations containing a bead (indicated by the arrow) which a can and b cannot be moved by the kink-jump. Note that for the move to be permitted the target site has to neighbor two other beads, while the other "disallowed" sites neighbor one bead only (see blue lines in the figure). The probability of the site occupancy by other, more distant beads is different in both cases. In red there are marked the possible trajectories of chain fragment composed of three beads passing through the site a allowed and b disallowed for the kink-jump move, whose numbers are equal to 20 and 25, respectively

Conformational entropy of a chain with fixed ends
In this section we present a simplified approach to analytical evaluation of the conformational entropy of an SAW chain in a state of conformation distorted by fixing the ends in a selected distance. Because the ideal chain model based on the RW chain can be solved exactly, our approach will thus be to consider the random walk first and then to discuss the impact of the excluded volume effect.
The number of elementary steps M needed to generate the ideal chain is simply the sum of elemental steps performed in different directions: where MX + , MY + and MZ + denote numbers of steps along x, y and z axes, in positive directions, while MX − , MY − and MZ − are the step numbers in the negative directions. The number of all possible conformations of the chain formed by M steps can be expressed as: There is no preferential direction in equilibrium. Now, let us assume that the chain is extended into+z direction. This elongation will be compensated by contraction in other directions. Assuming the same extent of contraction in all directions the number Ω of all possible conformations reads: Equation 9 can be rewritten in a more general form: where Γ denotes the Euler gamma function. The number of steps should fulfill the following requirements: where L is the average z coordinate of one end bead in the deformed chain if the other end bead is located at z=0. Replacing the number of steps, M, by the number of beads N (M=N-1) and generalizing the above reasoning over different lattice types, one can write the following expression for the conformational entropy of a stretched chain where ω eff (equal to 6 for RW model) can be generalized to the effective coordination number of any other model lattice.
The entropy of an undeformed chain is maximum (S max ) for L=0 when the chain assumes the relaxed coil-like conformation. Applying the Stirling approximation this entropy S max is given by: where C is a constant given by: Table 2 presents the estimates of S max obtained on the basis of Eq. 15 for several models of polymer chains: RW, NRRW and SAW. For the sake of comparison there are also the corresponding values S max * calculated using other approaches (for details see the footnote to Table 2). As seen, in all cases the results show a good agreement (the greater N the better the coincidence between the results). This agreement allows one to assume that Eq. 13 can be also applied for the approximate estimation of conformational entropy of Fig. 7 The dependencies of P E (K) and P E (C) on L H . The horizontal lines indicate the corresponding probabilities determined for the free chain-P E (K) 0 and P E (C) 0 , respectively (SAW model, N=100) uniaxially deformed (elongated or shrunk) chains, at least those modeled by RW, NRRW and SAWs.
Equation 13 relates the chain entropy to the average coordinate L, or, in other words, the magnitude of the average value of the projection of the end-to-end distance L H onto a particular axis (here, onto the z-axis).
In the case of ideal chain the L H distance should fulfill two requirements: i) for an unperturbed chain, for which L=0, the distance L H should be equivalent to the rootmean-squared (rms) end-to-end distance 〈R H 2 〉 1/2 which is equal to (bM 1/2 ) [42] ii) for a fully stretched chain the distance L is equivalent to the contour length of the chain Assuming that there is a linear dependence between L and L H and generalizing the above considerations to any chain model one obtains: where 〈R H 2 〉 ν corresponds to the rms end-to-end distance of an unperturbed chain and v denotes the Flory's exponent [27]. The value 〈R H 2 〉 ν is related to the number of steps M by the relationship: The application of Eq. 22 for determination of the dependence between L and L H requires the knowledge of Flory's exponent for a given model or polymer/solvent system. For the real chains (i.e., with the excluded volume interactions) this exponent is equal to 3/5 for good solvents, 1/2 for theta solvents and 1/3 for poor solvents [27]. Figure 8 presents the dependencies of the conformational entropy of a chain on the distance L H between its ends, calculated on the basis of Eqs. 13 and 22 for some chosen chain models. As seen, the curves go through maxima at L H =〈R H 2 〉 ν , which is in accordance with expectation.
In the calculation of S=f(L H ) for the SAW model, which is presented in Fig. 8, a constant value of ω eff =4.6839 (valid for N→∞) was assumed. However, for the chains with excluded volume interactions the value of ω eff depends on L H , because the probability of nearest-neighbor site ) for the NRRW model: where C'=1.17, γ=7/6 and ω eff,∞ =4.6838 for the 3D linear, infinitely long macromolecule in the athermal solution [38][39][40][41] Fig. 8 Dependencies of conformational entropies vs. fixed end-to-end distance calculated for 2D RW, 3D RW and 3D SAW (N =100). Perpendicular lines indicate rms end-to-end distances of RW and SAW chains equal to 10b and 15.8b, respectively. The rms end-toend distance for RW model is independent on the dimensionality of the lattice occupancy changes with the extent of chain elongation. And thus, one can easily deduce that when the chain is highly stretched the excluded volume effect becomes negligibly small and ω eff tends to 5 if the simple cubic lattice is considered or to ω-1 in a general case.
Unfortunately, except for the limiting case of fully extended chain, the determination of ω eff is not straightforward since there is no analytical way of finding this parameter. In the next section of the paper we propose a method for the estimation of ω eff on the basis of selected observables which are directly accessible from MC simulations. Moreover, to extend the applicability of Eqs. 13 and 22 to the case of chains in constrained environments we propose the employment of such observables in the quantification of the degree of chain deformation instead of L H parameter. The meaning of L H in the case of confined chains is sometimes unclear and, what is more important, not corresponding to the chain deformation degree. For instance, when a linear macromolecule is confined in a L-or U-shaped nanochannel, the straight line connecting ends of the molecule may cross the channel walls and its length will not correspond to the extent of the molecule elongation.

Calculation of L H and ω eff parameters from simulation results
End-to-end distance As the probability P S (K) is related to a selected separation between chain ends (Fig. 4) one can try to apply the value of this probability for the estimation of L or L H parameters. For all the models studied here (2D RW, 3D RW, 2D SAW and 3D SAW) the dependencies of L versus the ratio P S (K)/P S (K) 0 are of similar character: at L=0 (free chain) the value of P S (K)/P S (K) 0 is approximately 1 and then with increasing L it gradually decreases, approaching zero at L= bM (as seen in Fig. 9a). The shape of the curves suggests an exponential relationship between these two variables. Actually, as can be seen in Fig. 9b, plots of (1-P S (K)/P S (K) 0 ) against L on log-log coordinates give nearly straight lines with slopes dependent on the model considered.
From the detailed analysis of the results collected in Fig. 9 the following relationship was derived: where the exponent 1 assumes values of 0.40, 0.40, 0.27 and 0.33, respectively for 2D RW, 3D RW, 2D SAW and 3D SAW. From the comparison of obtained values of 1 with the corresponding Flory's exponents (and keeping in mind that for RW model the Flory's exponent is independent of the lattice dimension [27,30], one can find the following simple relation between 1 and v: In Fig. 10a the values of L calculated on the basis of probabilities P S (K) obtained from simulations are compared with the corresponding values assumed for different chain models. This comparison indicates a good agreement between the retrieved and assumed values and thus allows a conclusion that Eq. 24 can be used for the calculation of parameter L. Although in the range of small P S (K) values, in which the equation becomes numerically unstable, there is a significant error in the L estimation. The described method of the determination of L is independent of the chain length, as can be concluded from inspection of Fig. 10b.

Effective coordination number
Equation 3 allows the calculation of the effective coordination number ω eff on the basis of the probability P E (R). However, the applicability of this method of ω eff determination is restricted to a free chain. Searching for an alternative method and having in mind the approximate relation 6 one can try to use P E (K) instead of P E (R) as a local probe of ω eff , and thus omit the requirement of free ends, necessary for P E (R) calculation. Assuming that the validity of Eq. 6 can be extended to deformed chains with fixed ends and by combining it with Eq. 3 one obtains: relates ω eff with P E (K). But, as could have been foreseen, its applicability is limited. For instance, it fails in predicting the results for boundary conditions. Namely, the solution of Eq. 26 does not satisfy the relations which can be written for two limiting states of the chain: fully stretched and fully shrunk. The corresponding relations are as follows: where κ→1 when N→∞ (i.e., when the number of all monomers located at the surface of shrunk coil is negligibly small as compared to N).
In the search for a relationship between ω eff and P E (K), which would be of a more general applicability than Eq. 26, we made use of the results of simulations performed for the models including intra-chain interactions, namely, for polymer chains in theta and poor solvents. The results obtained for such models allow direct and independent determination of both ω eff and P E (K). The attractive intra-chain interactions, which exist in theta (χ=1/2) and poor (χ>1/2) solvents, modify the value of ω eff as compared to that in athermal environment. Moreover, these interactions change the internal energy U of the model to be non-zero. The average value of U is given by: where ε ij denotes the interaction energy of the ith bead with an adjacent jth lattice site. The contact interaction between the beads i and j is taken into account only when they are on the adjacent lattice sites, but not on the adjacent sites along the chain. Depending on whether the site is empty or occupied by another bead, the value of ε ij is equal to ε PS or (ε PP + ε SS )/2, respectively. Because the variables ω eff and P E (K) are correlated (the occupancy of a chain adjacent site by another bead involves a simultaneous decrease in ω eff and increase in U), the function describing this correlation can be used for the calculation of ω eff .
Using the values of U determined from the simulations in which ε PS =ε SS =0 was assumed, the average value of ω eff can be calculated form the following equation: The values of ω eff obtained in this way are summarized in Fig. 11, which shows the relation between ω eff and P E (K) in the form of a function ((ω-1)-ω eff ) 1/2 =f(1-P E (K)). As seen, all the data points lie quite well along a straight line whose slope is equal to 2.0, indicating that the searched relationship between ω eff and P E (K) is well approximated by the equation It is noteworthy that Eq. 31 fulfills Eqs. 27 and 28. It follows from Fig. 12a that, as expected, when L H tends to the contour length of the chain, the value of ω eff calculated on the basis of Eq. 30 tends to ω-1. The most distinct decrease in ω eff accompanying the decrease in L H is seen in the case of a chain in poor solvent conditions, which is also as expected. The results collected in Fig. 12b indicate Fig. 10 The plots of the values of L calculated on the basis of probabilities P S (K) against the assumed L values: a for different chain models, b for two different chain lengths (3D SAW). The slope of the straight line is equal to 1 Fig. 11 The values of the expression ((ω-1)-ω eff ) 1/2 plotted against 1-P E (K), obtained for a chain of N=100 in theta and poor solvents. The regression coefficient of the straight line indicated in the figure is equal to 2 that the chain length has only a minor influence on the obtained ω eff values, at least in the studied N range (N=50 and N=100).

Conformational entropy and energy of the chain in exemplary systems
In this section, the results of calculations of the conformational entropy and energy of a polymer, obtained by means of the method based on the use of probabilities P S (K) and P E (K), are presented. The calculations were performed for different selected end-to-end distances and solvent conditions, as well as in the presence of geometrical constraints. Figure 13a presents the entropies S of a SAW chain of N=50 obtained at different L H values. These results are compared with the corresponding results available in literature and obtained with the expanded ensemble Monte Carlo method [34]. Practically, in the whole range of L H values a very good agreement is observed. It is also noteworthy that as expected, the maximum on the curve S(L H ) occurs at the L H value corresponding to the rms distance between ends of a free chain: R H /bM ≈0.20 for SAW of N=50 and v=3/5 (Eq. 23). For the sake of comparison with the results reported in ref. [34], the stretching forces F required to    extend the chain to assumed end-to-end separations were also calculated. The corresponding values of F were found from the following relation: As illustrated in Fig. 13b, the obtained results are in fairly good agreement with those of ref. [34].
Using the method described in this paper, the values of Helmholtz free energies A of the chain of N=100 in different solvent conditions (i.e., of different values of Flory's interaction parameter: χ=0, χ=1/2 and 1/2<χ≈2) were calculated from: The obtained results, together with corresponding results for U and S, are shown in Fig. 14. As seen, for large extensions (large L H ) the conformational entropy of the chain is practically independent of the solvent condition (Fig. 14a). It can be explained by the fact that in the case of strongly stretched chain the probability of bead-bead contacts (which mimic inter-monomer interactions) is very small and thus its influence on the chain entropy is also small. For intermediate and low extensions the entropy of the chain in solvents of χ≥1/2 is smaller than in the athermal condition, because of smaller values of ω eff (compare Fig. 12a). As seen in Fig. 14b, in non-athermal conditions the decrease in L H is accompanied by a monotonic decrease in internal energy of the chain. The superposition of energetic and entropic contributions produces the A(L H ) dependencies exhibiting minima (Fig. 14c). In spite of the low precision of estimation of positions of these minima (which is because of relatively small length of the examined chain, N=100), it can be found out that they occur at L H values roughly equal to rms end-to-end distances of the free chain in corresponding solvent conditions (the values of R H /bM calculated from Eq. 23 are 0.16, 0.10 and 0.05 for athermal, theta and poor solvent conditions respectively).
The described method of estimation of the conformational entropy of macromolecules was also tested for a chain in the space with geometrical constraints. The entropies of SAW chain with its ends attached to two opposite parallel plane surfaces which were also assumed to be impenetrable, infinite and non-interacting with the chain except its terminals were calculated at different separations between the surfaces. The results obtained are collected in Fig. 15 and compared with the entropies found for the chain in the open space when selected distance L H between ends was equal to that in the presence of the surfaces. For the confined chain in the range of small plane separations (L H +2b≤〈R H 2 〉 v ) the calculated entropies are clearly smaller as compared to those of the chain in the open space. This can be elucidated by the excluded volume effect of the surfaces, which reduces the value of ω eff . From Fig. 15, it can also be deduced that for a chain with free ends the average end-to-end distance in the equilibrium conformation is longer in the presence of the surfaces.

Conclusions
The probabilities of different types of micromodifications of chain conformation which were found from the performed MMC simulations agree reasonably well with the results obtained on the basis of approximate theoretical considerations regarding the nearest environment of a local chain conformation. The proposed method for the estimation of conformational entropy of a polymer chain using the probabilities of kink-jump moves, P S (K) and P E (K), provides results consistent with literature [34] when applied to a deformed chain (i.e., whose two ends are fixed over a distance L H different from the most probable end-to-end separation of an undisturbed chain). In the range of small deformations the entropy estimation accuracy was low because of the numerical instability of Eq. 31. The dependencies of the conformational entropy S and free energy A of the chain on the imposed end-to-end distance L H , calculated on the basis of probabilities P S (K) and P E (K) for different solvent conditions as well as in the presence of the geometric constraints in the simulation space, qualitatively capture major features of the expected behavior of a polymer chain under such "test" conditions.