Monte Carlo Simulations of Lattice Models for Single Polymer Systems

Monte Carlo Simulations of Lattice Models for Single Polymer Systems

Abstract

Single linear polymer chains in dilute solutions under good solvent conditions are studied by Monte Carlo simulations with the pruned-enriched Rosenbluth method up to the chain length . Based on the standard simple cubic lattice model (SCLM) with fixed bond length and the bond fluctuation model (BFM) with bond lengths in a range between and , we investigate the conformations of polymer chains described by self-avoiding walks (SAWs) on the simple cubic lattice, and by random walks (RWs) and non-reversible random walks (NRRWs) in the absence of excluded volume (EV) interactions. In addition to flexible chains, we also extend our study to semiflexible chains for different stiffness controlled by a bending potential. The persistence lengths of chains extracted from the orientational correlations are estimated for all cases. We show that chains based on the BFM are more flexible than those based on the SCLM for a fixed bending energy. The microscopic differences between these two lattice models are discussed and the theoretical predictions of scaling laws given in the literature are checked and verified. Our simulations clarify that a different mapping ratio between the coarse-grained models and the atomistically realistic description of polymers is required in a coarse-graining approach due to the different crossovers to the asymptotic behavior.

hsu@mpip-mainz.mpg.de

1Introduction

In the theoretical study of polymer physics [1], computer simulations provide a powerful method to mimic the behavior of polymers covering the range from atomic to coarse-grained scales depending on the problems one is interested in [3]. The generic scaling properties of single linear and branched polymers in the bulk or confinement under various solvent conditions have been described quite well by simple coarse-grained lattice models (i.e., random walks (RWs), non-reversible random walks (NRRWs), self-avoiding random walks (SAWs), or interacting self-avoiding random walks (ISAWs) on a regular lattice, regarding the interactions between non-bonded monomers). As an alternative one can use coarse-grained models in the continuum, such as a bead-spring model (BSM) (where all beads interact with a truncated and shifted Lennard-Jones (LJ) potential while the bonded interactions are captured by a finitely extensible nonlinear elastic (FENE) potential) using Monte Carlo and molecular dynamics simulations [3]. On the one hand, however, as the size and complexity of a system increases, detailed information at the atomic scale may be lost when employing low resolution coarse-graining representations. On the other hand, the cost of computing time may be too high if the system is described at high resolution. Therefore, more scientific effort has been devoted to developing an appropriate coarse-grained model which can reproduce the global thermodynamic properties and the local mechanical and chemical properties such as the intermolecular forces between polymer chains [5]. While these models are already known since a long time, the present work is the first study presenting precise data on conformational properties of these model, when a bond angle potential is included.

In this paper we deal with linear polymer chains in dilute solutions under good solvent conditions, and describe them by lattice models on the simple cubic lattice. Although coarse-grained lattice models neglect the chemical detail of a specific polymer chain and only keep chain connectivity (topology) and excluded volume, the universal behavior of polymers still remains the same in the thermodynamic limit (as the chain length ) [2], Two coarse-grained lattice models, the standard simple cubic lattice model (SCLM) and the bond fluctuation model (BFM) [3], are considered for our simulations. The SCLM is often used for the test of new simulation algorithms, and the verification of theoretically predicted scaling laws due to its simplicity and computational efficiency. The BFM has the advantages that the computational efficiency of lattice models is kept and the behavior of polymers in a continuum space can be described approximately. The model thus introduces some local conformational flexibility while retaining the computational efficiency of lattice models for implementing excluded volume interactions by enforcing a single occupation of each lattice vertex.

redThe excluded volume effect plays an essential role in any real polymer chain, while in a dilute solution under a theta solvent condition, or in concentrated polymer solutions such as melts, and glasses, the real polymer chain behaves like an ideal chain. The excluded volume constraint can easily be incorporated in the lattice models by simply forbidding any two effective monomers occupying the same lattice site (cell). Tries et al. have successfully mapped linear polyethylene in the melt onto the BFM and their results are in good agreement with experimental viscosimetric results quantitatively without adjusting any extra parameters [15]. Varying the backbone length and side chain length of the bottle-brush polymer based on the BFM, a direct comparison of the structure factors between the experimental data for the synthetic bottle-brush polymer consisting of hydroxyethyl methacrylate (PMMA) as the backbone polymer and flexible poly(n-butyl acrylate) (PnBA) as side chains in a good solvent (toluene) and the Monte Carlo results is given [17]. Furthermore, the lattice models have also been used widely to investigate the conformational properties of protein-folding [18] and DNA in chromosomes [19] in biopolymers. For alkane-like chains the angles between subsequent effective bonds are not continuously distributed, but only discrete angles are allowed. Therefore, the lattice model such as the SCLM where only the discrete angles and are allowed is an ideal model for studying alkane-like chains of different stiffnesses.

A direct comparison of simulation results between the lattice models and the off-lattice BSM is also possible, e.g. linear polymers [22] and ring polymers [23] in a melt, adsorption of multi-block and random copolymers [24], semiflexible chains under a good solvent condition [25], the crossover from semiflexible polymer brushes towards semiflexible mushrooms as the grafting density decreases [26]. Results from these two coarse-grained models are qualitatively the same. Namely, they both show the same scaling behavior, but the amplitudes and the critical or the crossover points can vary depending on the underlying models. However, our main motivation was to understand the microscopic differences between the two lattice models, SCLM and BFM, for describing the conformations of polymers and to provide this information for the further development of a multi-scale approach for studying polymers of complex topology and polymer solutions at high concentration based on the lattice models. Therefore, we only focus on the two coarse-grained lattice models, SCLM and BFM, here. In the mapping between atomistic and coarse-grained models, it turns out that a bond angle potential also needs to be included in the coarse-grained models, as is done here.

The outline of the paper is as follows: Section 2 describes the models and the simulation technique. Section 3 and Section 4 review the properties of flexible chains and semiflexible chains, respectively. Polymer chains described by SAWs, and by RWs and NRRWs in the absence of excluded volume effects are studied and compared with theoretical predictions. Finally our conclusions are summarized in Section 5.

2Models and simulation methods

The basic characteristics of linear polymer chains depend on the solvent conditions. Under good solvent conditions the repulsive interactions (the excluded volume effect) and entropic effects dominate the conformation, and the polymer chain tends to swell to a random coil. In the thermodynamic limit, namely the chain length , the partition function scales as

where is the critical fugacity per monomer, is the effective coordination number, and is the entropic exponent related to the topology. In two dimensions [2] , while the best estimate [27] for is . For the standard self-avoiding walks on the simple cubic lattice [28] in one has and the corresponding effective coordination number . The conformations of polymer chains characterized by the mean square end-to-end distance, , and the mean square gyration radius, , scale as [29]:

and

where is the Flory exponent, is the leading correction to the scaling exponent, and are non-universal constants, and is the mean square bond length. The quantities , , and the ratio are universal [31], while the quantities, , , , and , depend on the microscopic realization. In one has , while in the most accurate estimate of the Flory exponent [32] . We use for our data analysis in this paper.

Figure 1: (a) Effective exponents \gamma_{\rm eff}^{(1)} and \gamma_{\rm eff}^{(2)} [computed from Eqs. () and () ] plotted versus N on a semi-log scale. \gamma=\lim_{N \rightarrow \infty} \gamma_{\rm eff}^{(1)}(N) = 
\lim_{N \rightarrow \infty} \gamma_{\rm eff}^{(2)}(N)=1.1578(6). (b) \ln Z_N + N \ln a -(\gamma-1) \ln N with \gamma=1.1578 determined from (a) plotted versus N on a semi-log scale. The best estimate of fugacity \mu=a=0.0117241395(75) is determined by the horizontal curve.
(a) Effective exponents \gamma_{\rm eff}^{(1)} and \gamma_{\rm eff}^{(2)} [computed from Eqs. () and () ] plotted versus N on a semi-log scale. \gamma=\lim_{N \rightarrow \infty} \gamma_{\rm eff}^{(1)}(N) = 
\lim_{N \rightarrow \infty} \gamma_{\rm eff}^{(2)}(N)=1.1578(6). (b) \ln Z_N + N \ln a -(\gamma-1) \ln N with \gamma=1.1578 determined from (a) plotted versus N on a semi-log scale. The best estimate of fugacity \mu=a=0.0117241395(75) is determined by the horizontal curve.
Figure 1: (a) Effective exponents and [computed from Eqs. () and () ] plotted versus on a semi-log scale. . (b) with determined from (a) plotted versus on a semi-log scale. The best estimate of fugacity is determined by the horizontal curve.

Two models are used for simulating linear polymers in the bulk under good solvent conditions. One is the standard SAW on the simple cubic lattice, effective monomers being described by occupied lattice sites, connected by bonds of fixed length . Each site can be visited only once, and thus the excluded volume interaction is realized. The other is the standard bond fluctuation model (BFM). On the simple cubic lattice each effective monomer of a SAW chain blocks all eight corners of an elementary cube of the lattice from further occupation. Two successive monomers along a chain are connected by a bond vector which is taken from the set {,, , , , } including also all permutations. The bond length is therefore in a range between and . There are in total bond vectors and different bond angles between two sequential bonds along a chain serving as candidates for building the conformational structure of polymers. The partition sum of a SAW of steps is

which is simply the total number of possible configurations consisting of monomers.

In the literature there are still no estimates of the fugacity and the entropic exponent for SAWs on the BFM. According to the scaling law of the partition sum , Eq. (Equation 1), the effective entropic exponent obtained from triple ratios [33]

is shown in Figure 1a. It gives . The fugacity is therefore determined by adjusting the value such that the curve of with becomes horizontal for very large (see Figure 1b). We obtain the fugacity and the corresponding effective coordination number listed in Table 1. In Figure 1a we also show the asymptotic behavior of the effective entropic exponent defined by

with our estimate of for comparison.

If the excluded volume effect is ignored completely, a polymer chain behaves as an ideal chain. It is well described by a random walk (RW), a walk that can cross itself or may trace back the same path, or by a non-reversal random walk (NRRW) where immediate back tracing is not allowed. The partition sums of RW and NRRW are given by

where is the coordination number. for the standard RW on the simple cubic lattice, and for the BFM [34]. The Flory exponent is for an ideal chain and its mean square gyration radius .

For the simulations of single RW, NRRW, and SAW chains we use the pruned-enriched Rosenbluth method (PERM) [35]. It is a biased chain growth algorithm with resampling and population control. In this algorithm a polymer chain is built like a random walk by adding one monomer at each step with a bias depending on the problem at hand, and each configuration carries its own weight. The population control at each step is made such that the “bad” configurations are pruned with a certain probability, and the “good” configurations are enriched by properly reweighting, until a chain has either grown to the maximum length of steps, , or has been killed due to attrition. A detailed description of the algorithm PERM and its applications is given in a review paper [36]. The algorithm has the advantage that the partition sum can be estimated very precisely and directly. It is also very efficient for simulating linear polymer chains up to very long chain lengths in dilute solution at and above the -point. Therefore, we apply the algorithm on the two lattice models, SCLM and BFM, in order to check for major differences between these two microscopic models. The longest chain length is in our simulations here.

Figure 2: Mean square end-to-end distance \langle R^2_e \rangle and gyration radius \langle R^2_g \rangle scaled by (N^\nu \ell_b^2) with the Flory exponent \nu=1/2 for RWs and NRRWs (a), and \nu=0.5876 for SAWs , plotted against N. Here the bond length \ell_b = \langle \vec{b}^2 \rangle^{1/2}: \ell_b=1 (SCLM) and \ell_b=2.72 (BFM).
Mean square end-to-end distance \langle R^2_e \rangle and gyration radius \langle R^2_g \rangle scaled by (N^\nu \ell_b^2) with the Flory exponent \nu=1/2 for RWs and NRRWs (a), and \nu=0.5876 for SAWs , plotted against N. Here the bond length \ell_b = \langle \vec{b}^2 \rangle^{1/2}: \ell_b=1 (SCLM) and \ell_b=2.72 (BFM).
Figure 2: Mean square end-to-end distance and gyration radius scaled by with the Flory exponent for RWs and NRRWs (a), and for SAWs , plotted against . Here the bond length : (SCLM) and (BFM).
Figure 3: Ratio between the mean square end-to-end distance and the mean square gyration radius, \langle R^2_e \rangle/\langle R^2_g \rangle, plotted against the chain segments N, for RWs and NRRWs (a), and for SAWs (b). As N \gg 0, \langle R_e^2 \rangle/\langle R_g^2 \rangle \sim 6.00 in (a) and \langle R_e^2 \rangle/\langle R_g^2 \rangle \sim 6.25 in (b).
Ratio between the mean square end-to-end distance and the mean square gyration radius, \langle R^2_e \rangle/\langle R^2_g \rangle, plotted against the chain segments N, for RWs and NRRWs (a), and for SAWs (b). As N \gg 0, \langle R_e^2 \rangle/\langle R_g^2 \rangle \sim 6.00 in (a) and \langle R_e^2 \rangle/\langle R_g^2 \rangle \sim 6.25 in (b).
Figure 3: Ratio between the mean square end-to-end distance and the mean square gyration radius, , plotted against the chain segments , for RWs and NRRWs (a), and for SAWs (b). As , in (a) and in (b).
Table 1: The estimates of fugacity , the effective coordination number in Eq. (), the amplitudes and in Eqs. () and () determined from the simulation data of and for RWs, NRRWs, and SAWs based on the two lattice models, SCLM and BFM.
model RW NRRW SAW RW NRRW SAW
1/6 1/5 0.21349098(5) 1/108 1/107 0.01172414395(75)
6 5 4.6840386(11) 108 107 85.294106(55)
1.0000(2) 1.4988(4) 1.220(3) 0.9986(2) 1.0714(2) 1.247(5)
0.16666(0) 0.24985(7) 0.1952(4) 0.16645(4) 0.16959(3) 0.1993(6)
Figure 4: (a),(c) Normalized probability distributions of end-to-end distance, h_N(R_e/\ell_b)(={\cal P}(R_e/\ell_b)), plotted versus R_e/\ell_b. (b),(d) similar as (a),(c), but for gyration radius R_g. Data are for RWs (a),(b) and SAWs (c),(d). Several values of chain lengths N are chosen, as indicated.
(a),(c) Normalized probability distributions of end-to-end distance, h_N(R_e/\ell_b)(={\cal P}(R_e/\ell_b)), plotted versus R_e/\ell_b. (b),(d) similar as (a),(c), but for gyration radius R_g. Data are for RWs (a),(b) and SAWs (c),(d). Several values of chain lengths N are chosen, as indicated.
(a),(c) Normalized probability distributions of end-to-end distance, h_N(R_e/\ell_b)(={\cal P}(R_e/\ell_b)), plotted versus R_e/\ell_b. (b),(d) similar as (a),(c), but for gyration radius R_g. Data are for RWs (a),(b) and SAWs (c),(d). Several values of chain lengths N are chosen, as indicated.
(a),(c) Normalized probability distributions of end-to-end distance, h_N(R_e/\ell_b)(={\cal P}(R_e/\ell_b)), plotted versus R_e/\ell_b. (b),(d) similar as (a),(c), but for gyration radius R_g. Data are for RWs (a),(b) and SAWs (c),(d). Several values of chain lengths N are chosen, as indicated.
Figure 4: (a),(c) Normalized probability distributions of end-to-end distance, , plotted versus . (b),(d) similar as (a),(c), but for gyration radius . Data are for RWs (a),(b) and SAWs (c),(d). Several values of chain lengths are chosen, as indicated.

3Conformations of single linear polymer chains: RWs, NRRWs, and SAWs

We employ the pruned-enriched Rosenbluth method (PERM) for the simulations of long single linear polymer chains of chain lengths (segments) up to , modeled by RWs, NRRWs, and SAWs depending on the interactions between monomers. Figure 2a with and Figure 2b with show that the scaling laws, Eqs. (Equation 2) and (Equation 3), are verified as one should expect. The mean square end-to-end distance simply is

The mean square gyration radius is given by

where is the center of mass position of the polymer. The amplitudes and for RWs, NRRWs, SAWs based on the two lattice models, SCLM and BFM, are listed in Table 1. Results of and for RWs from both models follow the same curves although the bond vectors in the BFM are not all along the lattice directions and do not have the same bond length. Here is the root-mean square bond length, for the SCLM and for the BFM. In Figure 2a, values of and for NRRWs, obtained from SCLM for all lengths are significant larger than that obtained from the BFM, since at each step the walker can only go straight or make a L-turn in the SCLM. In Figure 2b, the two curves showing the results of [] with as functions of obtained from the two models intersect at , and finally the amplitude for BFM is larger in the asymptotic regime. The slight deviation from the plateau value is due to the finite size effects. The correction exponent [Eqs. (Equation 2) and (Equation 3)] for these two models is determined by plotting and versus (not shown). One should expect straight lines near if and only if . We obtain for both models, which is in agreement with the previous simulation in Ref. [29] within error bars. The ratio between the mean square end-to-end distance and the mean square gyration radius, is indeed for RWs and NRRWs (Fig. Figure 3a). For SAWs our results give (Fig. Figure 3b). For SAWs on the simple cubic lattice the most accurate estimates of , , and are given in Ref. [32]. Our results are also in perfect agreement with them. However, much longer chain lengths will be needed for a more precise estimate of the plateau value of the ratio in the asymptotic scaling regime. Note that for the behavior is clearly model-dependent.

We include here the RW and NRRW versions of both models not just for the sake of an exercise: often the mapping from an atomistic to a coarse-grained model is to be done under melt conditions, where excluded volume interactions are screened.

The shapes of polymer chains can also be described by the probability distributions of end-to-end distance and gyration radius, and , respectively. Numerically, they are obtained by accumulating the histogram of over all configurations of length , given by

here each configuration carries its own weight . The normalized histogram is therefore,

Results of and for RWs and SAWs obtained from the two models for various values of chain lengths are shown in Figure 4. We see that both models give for the same distributions of and although the mean values of and are slightly different between these two lattice models (Fig. Figure 2b) for SAWs. Note that an angular average over all directions has been included in the accumulating process of the histogram due to spherical symmetry. Thus, the normalized histograms of ,

with

and the normalized histograms of ,

with

where and are the normalization factors.

The probability distribution of end-to-end distance for ideal chains is simply a Gaussian distribution,

Our numerical data for RWs obtained from BFM and SCLM shown in Figure 4 are in perfect agreement with the Gaussian distribution (see Figure 5). From Eqs. (Equation 10), (Equation 6) and (Equation 7) we obtain the normalized factor .

Figure 5: Same as in Fig. a, but data are only for BFM. The predicted distribution 4 \pi (R_e/\ell_b)^2P(R_e/\ell_b) for various values of N are shown by solid curves. Here the distribution function P(R_e/\ell_b) is given by Eq. ().
Figure 5: Same as in Fig. a, but data are only for BFM. The predicted distribution for various values of are shown by solid curves. Here the distribution function is given by Eq. ().
Figure 6: Root mean square radius of gyration, \langle R_g^2 \rangle^{1/2}, and the gyration radius R_{g,m} where P(R_g) has its maximum value, plotted against chain length N for RWs and SAWs. Data are for BFM.
Figure 6: Root mean square radius of gyration, , and the gyration radius where has its maximum value, plotted against chain length for RWs and SAWs. Data are for BFM.
Figure 7: Logarithm of the rescaled probability distribution of gyration radius, \ln(P(R_{g,m})/P(R_g)) as a function of R_g/R_{g,m} for SAWs obtained from the models (a) SCLM and (b) BFM. The best fit of Eq. () with A=1.18 to our data is shown by the solid curve. The dashed curve with A=1.34 given in Ref.  is also shown for the comparison. Several values of chain lengths N are chosen, as indicated.
Logarithm of the rescaled probability distribution of gyration radius, \ln(P(R_{g,m})/P(R_g)) as a function of R_g/R_{g,m} for SAWs obtained from the models (a) SCLM and (b) BFM. The best fit of Eq. () with A=1.18 to our data is shown by the solid curve. The dashed curve with A=1.34 given in Ref.  is also shown for the comparison. Several values of chain lengths N are chosen, as indicated.
Figure 7: Logarithm of the rescaled probability distribution of gyration radius, as a function of for SAWs obtained from the models (a) SCLM and (b) BFM. The best fit of Eq. () with to our data is shown by the solid curve. The dashed curve with given in Ref. is also shown for the comparison. Several values of chain lengths are chosen, as indicated.
Figure 8: Similar as in Fig. , but for RWs. The best fitting of our data gives A=0.97 for both models.
Similar as in Fig. , but for RWs. The best fitting of our data gives A=0.97 for both models.
Figure 8: Similar as in Fig. , but for RWs. The best fitting of our data gives for both models.
Figure 9: Normalized probability distributions of R_g/\ell_b, h(R_g/\ell_b)={\cal P}_N(R_g/\ell_b), plotted versus R_g/\ell_b for various values of N, and for BFM. The fitting functions 4\pi C_{g,N}(R_g/\ell_b)^2 P(R_g/\ell_b) [Eqs. (), (), and ()] with parameters a_1, a_2, and C_{g,N} determined by method (1) and method (2) are shown by curves for comparison in (a) and (b), respectively.
Normalized probability distributions of R_g/\ell_b, h(R_g/\ell_b)={\cal P}_N(R_g/\ell_b), plotted versus R_g/\ell_b for various values of N, and for BFM. The fitting functions 4\pi C_{g,N}(R_g/\ell_b)^2 P(R_g/\ell_b) [Eqs. (), (), and ()] with parameters a_1, a_2, and C_{g,N} determined by method (1) and method (2) are shown by curves for comparison in (a) and (b), respectively.
Figure 9: Normalized probability distributions of , , plotted versus for various values of , and for BFM. The fitting functions [Eqs. (), (), and ()] with parameters , , and determined by method (1) and method (2) are shown by curves for comparison in (a) and (b), respectively.

The theoretical prediction of the gyration radius probability distribution of polymer chains under good solvent conditions in -dimensions suggested by Lhuillier [37] is as follows:

where and are (non-universal) constants, and the exponents and are linked to the space dimension and the Flory exponent by

Here is the des Cloizeaux exponent [38] for the osmotic pressure of a semidilute polymer solution, and is the Fisher exponent [39] characterizing the end-to-end distance distribution.

This scaling form has been verified in the previous Monte Carlo simulation studies of the standard self-avoiding walks on square () and cubic () lattices up to steps using the slithering-snake and pivot algorithms [40]. The two fitting parameters and are actually not independent since at the position where the distribution has its maximum value, i.e. , the corresponding gyration radius (see Figure 6). Using Eq. (Equation 11), the logarithm of the rescaled probability is written as

with

From Eq. (Equation 9), we obtain

Our estimate of for SAWs based on the two lattice models, SCLM and BFM, are shown in Figure 7. As chain lengths , we see the nice data collapse, and the logarithm of the scaled probability of is described by Eq. (Equation 12) with very well. Due to the finite-size effect it is clearly seen that the previous estimate is an overestimate [41].

For an ideal chain the distribution of is no longer a simple Gaussian distribution as shown in Eq. (Equation 10), and the exact expression is quite complicated. Vettorel et al. [10] found out the formula given by Lhuillier [37] is a good approximation for describing the distribution for an ideal chain based on the BSM. Therefore, we also use the same formula for the investigation of the distribution obtained from the two coarse-grained lattice models. Two methods are discussed here. Method (1): We use the formula [Eq. (Equation 12)] which contains only one fitting parameter since for RWs as shown in Figure 6. From our simulations of RWs, we still see the nice data collapse for in the plot of the logarithm of the rescaled distribution of (Fig. Figure 8), but the distribution can only be described by Eq. (Equation 12) well for . Using the least square fit, it gives . Method(2): We assume that the two parameters and in Eq. (Equation 11) are independent. Using Eqs. (Equation 8), (Equation 9), and (Equation 11), values of , , and the normalization factor are determined by the best fit of the normalized histograms obtained from our Monte Carlo simulations. Note that it is not possible to determine and using the second method for since the normalization condition, Eq. (Equation 9), is not satisfied. In Figure 9 we compare our results of for BFM to the fitting function [Eqs. (Equation 8), (Equation 9), and (Equation 11)] with parameters determined by these two different methods. Values of and plotted versus are shown in Figure 10 and listed in Table 2. Our results show that and are almost constants for large , which are comparable with the results obtained for the BSM [10].

Figure 10: Parameters a_1 and a_2 in Eq. () plotted versus N. Results are obtained from two different methods mentioned in the text. Here a_{1,\infty} and a_{2,\infty} are taken from Table  for N=10000.
Figure 10: Parameters and in Eq. () plotted versus . Results are obtained from two different methods mentioned in the text. Here and are taken from Table for .
Table 2: Parameters and of the probability distribution , Eq. (), determined by two different methods (1) and (2) mentioned in the text for various values of chain length .
10 20 50 100 200 400 800 1000 2000 4000 8000 10000
(1) 7.12 5.27 4.94 4.52 4.70 4.70 4.82 4.34 4.88 4.74 4.37 4.68
(1) 12.80 14.15 14.46 14.90 14.70 14.73 14.58 15.10 14.42 14.66 15.07 14.73
(2) _ _ 4.14 3.83 3.69 3.60 3.57 3.56 3.58 3.54 3.53 3.53
(2) _ _ 13.25 13.27 13.29 13.29 13.30 13.29 13.35 13.29 13.29 13.28

4Semiflexible chains

We extend our simulations in this section from flexible chains to semiflexible chains. Extensive Monte Carlo simulations of semiflexible polymer chains described by standard SAWs on the simple cubic lattice, with a bending potential , have been recently carried out [42]. Recall that atomistic models of real chains may exhibit considerable chain stiffness due to the combined action of torsional and bond angle potentials. When a mapping to a coarse-grained model is performed, this stiffness is lumped into an effective bond angle potential of the coarse-grained model. In this standard model the angle between two subsequent bond vectors along the chain is either or , and hence in the statistical weight of a SAW configuration on the lattice every bend will contribute a Boltzmann factor ( for ordinary SAWs). is of order unity throughout the whole paper. The partition function of such a standard SAW with bonds ( effective monomers) and bends is therefore,

where is the total number of all configurations of a polymer chain of length containing bends.

We are also interested in understanding the microscopic difference between the standard SAWs and the BFM as the stiffness of chains is taken into account. Since there are bond angles possibly occurring in the chain conformations, the partition function cannot be simplified for the BFM, written as,

where is the bond angle between the bond vector and the bond vector along a chain, and is the number of configurations having the same set but fluctuating bond lengths. In the absence of excluded volume effect, the formulas of the partition function, Eq. (Equation 13) and Eq. (Equation 14), remain the same while semiflexible chains are described by RWs and NRRWs.

4.1Theoretical predictions

There exist several theoretical models describing the behavior of semiflexible chains in the absence of excluded volume effects. We first consider a discrete worm-like chain model [45] that a chain consisting of bonds with fixed bond length , but successive bonds are correlated with respect to their relative orientations,

where is the angle between the successive bond vectors. The mean square end-to-end distance is therefore,

This formula agrees with the prediction for a freely rotating chain (FRC). In the limit the bond-bond orientational correlation function decays exponentially as a function of their chemical distance ,

where is the so-called persistence length which can be extracted from the initial decay of . Equivalently, one can calculate the persistence length from

here instead of we use to distinguish between these two measurements.

For rather stiff () and long chains () we expect that the bond angles between successive bonds along chains are very small (), then Eqs. (Equation 16) and (Equation 18) become

and

Eq. (Equation 19) is equivalent to the mean square end-to-end distance of a freely jointed chain that Kuhn segments of length are jointed together,

being the contour length and in this limit.

In the continuum limit , , but keeping and finite, we obtain from Eq. (Equation 16) the prediction for a continuous worm-like chain,

It gives the same result as that derived directly from the Kratky-Porod model [46] for worm-like chains in ,

where the polymer chain is described by the contour in continuous space. Equation (Equation 21) describes the crossover behavior from a rigid-rod for , where , to a Gaussian coil for , where as shown in Eq. (Equation 20).

For semiflexible Gaussian chains the contour length can also be written as and the mean square end-to-end distance and gyration radius described in terms of and are [46]

One can clearly recognize that Gaussian behavior of the radii is only seen, if the number of the persistence length that fits to a given contour length is large, , while a crossover to rigid-rod behavior occurs for of order unity.

In recent works in Ref. [22], authors have shown that the exponential decay of the bond-bond orientational correlation function, Eq. (Equation 17), and the Gaussian coil behavior, Eq. (Equation 21) for , predicted by the worm-like chain model only hold for and up to some values and , respectively when excluded volume effects are considered. The predictions of a theory based on the Flory-type free energy minimization arguments [2] proposed as an alternative to semiflexible chains with excluded volume interactions have been verified. In this treatment one considers a model where rods of length and diameter are jointed together, such that the contour length . Apart from prefactors of order unity, the second virial coefficient in then can be estimated as

The free energy of a chain now contains two terms, the elastic energy taken as that of a free Gaussian, i.e., , and the repulsive energy due to interactions treated in mean field approximation, i.e. proportional to the square of the density and the volume . Hence,

Minimizing with respect to , we obtain for the standard Flory result

Eq. (Equation 25) holds also for finite and since the contribution of the second term in Eq. (Equation 24) is still important. For the first term in Eq. (Equation 24) dominates, and the chain behaves as a Gaussian coil, , while for even smaller , , the chain behaves as a rigid-rod. Thus, the double crossover behavior of the mean square end-to-end distance is summarized as follows,

Figure 11: Ratio between mean square end-to-end distance and mean square gyration radius, \langle R_e^2 \rangle/\langle R_g^2 \rangle, plotted against chain lengths (segments) N=L/\ell_b for SCLM (a) and for BFM (b). Data for semiflexible chains described by NRRWs and RWs are shown by symbols and lines, respectively. Here \ell_b=1 for SCLM, and \ell_b=2.72 for BFM.
Ratio between mean square end-to-end distance and mean square gyration radius, \langle R_e^2 \rangle/\langle R_g^2 \rangle, plotted against chain lengths (segments) N=L/\ell_b for SCLM (a) and for BFM (b). Data for semiflexible chains described by NRRWs and RWs are shown by symbols and lines, respectively. Here \ell_b=1 for SCLM, and \ell_b=2.72 for BFM.
Figure 11: Ratio between mean square end-to-end distance and mean square gyration radius, , plotted against chain lengths (segments) for SCLM (a) and for BFM (b). Data for semiflexible chains described by NRRWs and RWs are shown by symbols and lines, respectively. Here for SCLM, and for BFM.
Figure 12: Semi-log plot of the bond-bond orientational correlation function \langle \cos \theta (s) \rangle vs. s\ell_b for SCLM with \ell_b=1 (a) and for BFM with \ell_b=2.72 (b). Data are for semiflexible chains described by SAWs, NRRWs and RWs and for \varepsilon_b=5.30, 4.61, 3.91, 3.51, 3.00, 2.30 from top to bottom in (a), and for \varepsilon_b=10, 7, 5, 3, and 1 from top to bottom in (b). The straight lines indicate fits of the initial decay, \langle \cos \theta(s) \rangle \propto \exp(-s\ell_b/\ell_p) [Eq. ()], for RWs.
Semi-log plot of the bond-bond orientational correlation function \langle \cos \theta (s) \rangle vs. s\ell_b for SCLM with \ell_b=1 (a) and for BFM with \ell_b=2.72 (b). Data are for semiflexible chains described by SAWs, NRRWs and RWs and for \varepsilon_b=5.30, 4.61, 3.91, 3.51, 3.00, 2.30 from top to bottom in (a), and for \varepsilon_b=10, 7, 5, 3, and 1 from top to bottom in (b). The straight lines indicate fits of the initial decay, \langle \cos \theta(s) \rangle \propto \exp(-s\ell_b/\ell_p) [Eq. ()], for RWs.
Figure 12: Semi-log plot of the bond-bond orientational correlation function vs. for SCLM with (a) and for BFM with (b). Data are for semiflexible chains described by SAWs, NRRWs and RWs and for , , , , , from top to bottom in (a), and for , , , , and from top to bottom in (b). The straight lines indicate fits of the initial decay, [Eq. ()], for RWs.

4.2Simulation Results

In order to investigate the scaling behavior of the ratio for semiflexible RWs and NRRWs we plot our data versus for several choices of the stiffness parameter (Fig. Figure 11). As increases, the data increase towards a maximum and then decrease towards a plateau where the prediction for ideal chains holds. At the location of the maximum of , , the corresponding maximum values is . The maximum move monotonously to larger values as chains become stiffer. The deviation between the data for RWs and NRRWs based on the SCLM decreases as the bending energy increases (Fig. Figure 11a), while it is negligible for the simulation data obtained based on the BFM (Fig. Figure 11b) in all cases.

Figure 12 shows the bond-bond orientational correlation function plotted versus the chemical distance covering the range from flexible chains to stiff chains characterized by for the models SCLM and BFM. We compare the data obtained for SAWs, NRRWs, and RWs for various values of . The intrinsic stiffness remains the same for SAWs, NRRWs, and RWs as is fixed. Results obtained from both models verify that the asymptotic exponential decay of is valid only if the excluded volume effect is neglected, i.e., for RWs and NRRWs. For semiflexible SAWs cannot be correct for [22], we rather have

As we have seen in Figure 12, the exponential decay is ill-defined for rather flexible SAWs. Using Eq. (Equation 18) as a definition of the persistence length we can still give the estimate of the persistence length which is approximately the same as the estimate of the decay length for moderately stiff chains and stiff chains. The estimates of and depending on using Eqs. (Equation 17) and (Equation 18) are listed in Table 3 and Table 4. RWs are more flexible than NRRWs, and NRRWs are more flexible than SAWs from the estimates of the persistence lengths and based on the SCLM. Using the BFM the persistence lengths are almost the same in all cases of for RWs and NRRWs, and they are smaller compared with the estimates for SAWs. Note that in Figure 12b data deviate slightly from the fitting straight lines describing the initial exponential decay for RWs and NRRWs as the bending energy increases, i.e., the stiffness of chains increases. For the problem is more severe. Therefore, one should be careful using the BFM for studying rather stiff chains. An alternative way to the determination of the persistence length would be given by the best fit of the mean square end-to-end distance of RWs or NRRWs to Eq. (Equation 21). A simple exponential decay is always found for the probability distribution of connected straight segments for semiflexible chains based on the SCLM [43], while large fluctuations are observed for semiflexible chains based on the BFM due to bond vector fluctuations and lattice artifacts [52]. This is the main reason why the different scenarios of the bond-bond orientational correlation functions between the SCLM and the BFM for stiff chains are seen in Figure 12. Figure 13 shows the locations and the heights of the maximum of (Fig. Figure 11) plotted versus the persistence length for semiflexible RWs based on the SCLM and BFM. Note that , , and all depend on which controls the stiffness of chains. We see that the dependence between and are the same for both models, while for the BFM is slightly larger than that for the SCLM for a fixed value of since chains based on the BFM are more flexible.

Table 3: Two estimates for the persistence length from Eq. () and from Eq. () for semiflexible RWs, NRRWs, and SAWs with various values of based on the SCLM (). Here in our simulations values of are chosen for convenience.
1.0 0.4 0.2 0.1 0.05 0.03 0.02 0.01 0.005
0.0 0.91 1.61 2.30 3.00 3.51 3.91 4.61 5.30
RW 0.84 1.54 2.83 5.37 8.73 12.95 25.67 51.38
NRRW 1.05 1.70 2.97 5.50 8.87 13.09 25.80 51.53
SAW 2.04 3.35 5.96 9.54 13.93 26.87 52.61
RW 0.84 1.54 2.83 5.36 8.73 12.95 25.66 51.37
NRRW 0.62 1.05 1.70 2.98 5.50 8.87 13.08 25.79 51.53
SAW 0.67 1.12 1.81 3.12 5.70 9.10 13.35 26.28 51.52
Table 4: Two estimates for the persistence length from Eq. () and from Eq. () for semiflexible RWs, NRRWs, and SAWs with various values of based on the BFM ().
0.0 1.0 2.0 3.0 5.0 7.0 10.0 15
RW 0.87 1.62 2.54 4.69 7.18 12.09 27.65
NRRW 0.87 1.62 2.54 4.69 7.18 12.09 27.65
SAW 1.91 2.78 4.94 7.39 12.37 27.93
RW 0.87 1.62 2.54 4.63 6.87 10.50 17.73
NRRW 0.21 0.87 1.62 2.54 4.63 6.87 10.50 17.73
SAW 0.61 1.11 1.80 2.65 4.68 6.90 10.52 17.75

The scaling plots for testing the applicability of the worm-like chain prediction, Eq. (Equation 21) and Eq. (Equation 23) to our data of and are shown in Figure 14. The persistence length in Eq. (Equation 21) for various values of are extracted from the exponential fit of Eq. (Equation 17) for NRRWs (see Tables Table 3 and Table 4). Since the worm-like chain model is formulated in the continuum, care has to be taken to correctly take into account the lattice structure of the present model, particularly in the rod limit. Assuming that a rigid rod consisting of monomers is located at , , , along the x-axis on the simple cubic lattice, the mean square gyration radius is:

Therefore, due to the lattice structure, the mean square gyration radius is rescaled by instead of in order to compare with the theoretical predictions in Figure 14c,d. For semiflexible RWs and NRRWs the data are indeed very well described by the worm-like chain model. As increases, we observe the crossover behavior from a rigid-rod regime to a Gaussian coil regime. The plateau value in the Gaussian regime corresponds to the persistence length in Figure 14a,b and in Figure 14c,d. For SAWs the deviation from the prediction becomes more prominent as chains are more flexible since the excluded volume effects are more important.

Figure 13: Location N_h and height h of the maximum of \langle R_e^2 \rangle/\langle R_g^2 \rangle (see Fig. ) plotted versus the persistence length \ell_p/\ell_b for RWs based on the SCLM and BFM.
Figure 13: Location and height of the maximum of (see Fig. ) plotted versus the persistence length for RWs based on the SCLM and BFM.
Figure 14: Log-log plots of rescaled mean square end-to-end distance \langle R_e^2 \rangle /(2 \ell_b L) (a),(b) and rescaled mean square gyration radius \langle R_g^2 \rangle /(2 \ell_b (L+2\ell_b)) (c),(d) versus {N=L/\ell_b} for semiflexible chains described by SAWs, NRRWs and RWs based on the SCLM with \ell_b=1 (a),(c) and BFM with \ell_b=2.72 (b),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for NRRW are taken from Tables  and .
Log-log plots of rescaled mean square end-to-end distance \langle R_e^2 \rangle /(2 \ell_b L) (a),(b) and rescaled mean square gyration radius \langle R_g^2 \rangle /(2 \ell_b (L+2\ell_b)) (c),(d) versus {N=L/\ell_b} for semiflexible chains described by SAWs, NRRWs and RWs based on the SCLM with \ell_b=1 (a),(c) and BFM with \ell_b=2.72 (b),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for NRRW are taken from Tables  and .
Log-log plots of rescaled mean square end-to-end distance \langle R_e^2 \rangle /(2 \ell_b L) (a),(b) and rescaled mean square gyration radius \langle R_g^2 \rangle /(2 \ell_b (L+2\ell_b)) (c),(d) versus {N=L/\ell_b} for semiflexible chains described by SAWs, NRRWs and RWs based on the SCLM with \ell_b=1 (a),(c) and BFM with \ell_b=2.72 (b),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for NRRW are taken from Tables  and .
Log-log plots of rescaled mean square end-to-end distance \langle R_e^2 \rangle /(2 \ell_b L) (a),(b) and rescaled mean square gyration radius \langle R_g^2 \rangle /(2 \ell_b (L+2\ell_b)) (c),(d) versus {N=L/\ell_b} for semiflexible chains described by SAWs, NRRWs and RWs based on the SCLM with \ell_b=1 (a),(c) and BFM with \ell_b=2.72 (b),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for NRRW are taken from Tables  and .
Figure 14: Log-log plots of rescaled mean square end-to-end distance (a),(b) and rescaled mean square gyration radius (c),(d) versus for semiflexible chains described by SAWs, NRRWs and RWs based on the SCLM with (a),(c) and BFM with (b),(d). Data for various values of are shown, as indicated. Solid curves refer to the theoretical prediction, Eq. (), for WLC. The values of the persistence length for NRRW are taken from Tables and .

Note that one should not consider the correction factor relative to the Kratky-Porod model in Eq. (Equation 26) as a “lattice artefact”: in a real stiff polymer (e.g. an alkane-type chain) one also has a sequence of discrete individual monomers (separated by almost rigid covalent bonds along the backbone of the chain) lined up linearly (like in a rigid rod-like molecule) over about the distance of a persistence length. Furthermore, we compare simulation results of the ratio multiplied by as a function of to the theoretical prediction, the ratio between Eq. (Equation 22) and Eq. (Equation 23), in Figure 15. We see the nice data collapse for RWs and NRRWs in the Gaussian regime () and the increase of deviations from the master curve as the stiffness of chains decreases in Figure 15a,b. The ratio as for a rigid-rod, while as for a Gaussian coil. For SAWs we still see the nice data collapse in Figure 15c,d, but in both rigid-rod and Gaussian coil regimes the deviations from the master curve become more prominent as chains are more flexible. For the deviation is due to the excluded volume effects, and finally as for SAWs. Note that in both models the ratio of the mean square end-to-end and gyration radii exceed its asymptotic value still significantly even if is as large as .

Figure 15: Semi-log plots of [(L+2\ell_b)/L]\langle R_e^2 \rangle /\langle R_g^2 \rangle versus n_p=L/\ell_p for semiflexible chains described by RWs and NRRWs (a),(b) and SAWs (c),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, the ratio between Eq. () and Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for RWs, NRRWs, and SAWs are taken from Tables  and , respectively.
Semi-log plots of [(L+2\ell_b)/L]\langle R_e^2 \rangle /\langle R_g^2 \rangle versus n_p=L/\ell_p for semiflexible chains described by RWs and NRRWs (a),(b) and SAWs (c),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, the ratio between Eq. () and Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for RWs, NRRWs, and SAWs are taken from Tables  and , respectively.
Semi-log plots of [(L+2\ell_b)/L]\langle R_e^2 \rangle /\langle R_g^2 \rangle versus n_p=L/\ell_p for semiflexible chains described by RWs and NRRWs (a),(b) and SAWs (c),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, the ratio between Eq. () and Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for RWs, NRRWs, and SAWs are taken from Tables  and , respectively.
Semi-log plots of [(L+2\ell_b)/L]\langle R_e^2 \rangle /\langle R_g^2 \rangle versus n_p=L/\ell_p for semiflexible chains described by RWs and NRRWs (a),(b) and SAWs (c),(d). Data for various values of \varepsilon_b are shown, as indicated. Solid curves refer to the theoretical prediction, the ratio between Eq. () and Eq. (), for WLC. The values of the persistence length \ell_p/\ell_b for RWs, NRRWs, and SAWs are taken from Tables  and , respectively.
Figure 15: Semi-log plots of versus for semiflexible chains described by RWs and NRRWs (a),(b) and SAWs (c),(d). Data for various values of are shown, as indicated. Solid curves refer to the theoretical prediction, the ratio between Eq. () and Eq. (), for WLC. The values of the persistence length for RWs, NRRWs, and SAWs are taken from Tables and , respectively.

Recently, Huang et al. [53] performed Brownian dynamics simulations on two-dimensional (2D) semiflexible chains described by a BSM including the excluded volume interactions. Varying the chain stiffness and chain length their results confirmed the absence of a Gaussian regime in agreement with the results from semiflexible SAWs based on the SCLM [44], and with observations from experiments of circular single stranded DNA adsorbed on a modified graphite surface [?]. The rescaled mean square end-to-end distance, , in terms of for both models on the lattice and in the continuum turns out to be universal from the rigid-rod regime up to the crossover regime () irrespective of the models chosen for the simulations. In the 2D SAW regime, different amplitude factors result from the different models [25].

In , we indeed see the nice data collapse for semiflexible RWs, NRRWs, and SAWs in the plot of versus (cf. Fig. a,b) from rod-like regime crossover to the Gaussian regime for (not shown), and the data obtained from the two lattice models are well described by the Kratky-Porod scaling function, Eq. (). For the BSM in the continuum we should expect the same universal behavior. Although for semiflexible SAWs the second crossover from the Gaussian regime to the SAW regime for is rather gradual and not sharp, the relationship [56] between the crossover chain length and the persistence length , , holds for these two models here. It would be interesting to check whether such a scaling law would also hold for the BSM.

The structure factor is an experimentally accessible quantity measured by neutron scattering. We therefore also estimate by where denote the positions of the monomers in a chain, and the structure factor is normalized such that . In order to compare the results of obtained for fully flexible RWs, NRRWs, and SAWs based on the SCLM and BFM, we plot versus ( for SCLM, and for BFM) in Fig. a. We see that for , while for the power law ( for SAWs, and for RWs and NRRWs) holds. The lattice artifact sets in at . Due to the local packing the first peak appears at for the SCLM, while at for the BFM as increases. In Fig. b we show the results for semiflexible SAWs of different stiffnesses based on the BFM. The Gaussian regime where for large values of and then crosses gradually over to as expected for rigid rods. [57].

Figure 16: (a)(b) Log-log plot of structure factor S(q) versus q\ell_b. Data are for fully flexible chains based on the SCLM and BFM in (a), and for semiflexible chains based on the BFM including 6 choices of the stiffness in (b). (c) Rescaled structure factor qLS(q) plotted versus ql_p. Data are for semiflexible chains based on the SCLM and BFM including 4 choices of the stiffness each. In (a)(b), the straight lines indicate the rod-like behavior at large q (slope =-1) the SAW behavior for flexible chains (slope =-1/\nu, with \nu=0.588), and the Gaussian behavior (slope =-2). In (c), the formulas proposed by Kholodenko {Eqs. ()-()}, the Debye function {Eq. ()} for Gaussian chains, and qLS(q) \rightarrow \pi for a rigid-rod  are also shown for comparison.
(a)(b) Log-log plot of structure factor S(q) versus q\ell_b. Data are for fully flexible chains based on the SCLM and BFM in (a), and for semiflexible chains based on the BFM including 6 choices of the stiffness in (b). (c) Rescaled structure factor qLS(q) plotted versus ql_p. Data are for semiflexible chains based on the SCLM and BFM including 4 choices of the stiffness each. In (a)(b), the straight lines indicate the rod-like behavior at large q (slope =-1) the SAW behavior for flexible chains (slope =-1/\nu, with \nu=0.588), and the Gaussian behavior (slope =-2). In (c), the formulas proposed by Kholodenko {Eqs. ()-()}, the Debye function {Eq. ()} for Gaussian chains, and qLS(q) \rightarrow \pi for a rigid-rod  are also shown for comparison.
(a)(b) Log-log plot of structure factor S(q) versus q\ell_b. Data are for fully flexible chains based on the SCLM and BFM in (a), and for semiflexible chains based on the BFM including 6 choices of the stiffness in (b). (c) Rescaled structure factor qLS(q) plotted versus ql_p. Data are for semiflexible chains based on the SCLM and BFM including 4 choices of the stiffness each. In (a)(b), the straight lines indicate the rod-like behavior at large q (slope =-1) the SAW behavior for flexible chains (slope =-1/\nu, with \nu=0.588), and the Gaussian behavior (slope =-2). In (c), the formulas proposed by Kholodenko {Eqs. ()-()}, the Debye function {Eq. ()} for Gaussian chains, and qLS(q) \rightarrow \pi for a rigid-rod  are also shown for comparison.
Figure 16: (a)(b) Log-log plot of structure factor versus . Data are for fully flexible chains based on the SCLM and BFM in (a), and for semiflexible chains based on the BFM including choices of the stiffness in (b). (c) Rescaled structure factor plotted versus . Data are for semiflexible chains based on the SCLM and BFM including choices of the stiffness each. In (a)(b), the straight lines indicate the rod-like behavior at large (slope ) the SAW behavior for flexible chains (slope , with