Article pubs.acs.org/JPCB

Quantitative Study of Fluctuation Effects by Fast Lattice Monte Carlo Simulations. V. Incompressible Homopolymer Melts Pengfei Zhang, Delian Yang, and Qiang Wang* Department of Chemical and Biological Engineering, Colorado State University, Fort Collins, Colorado 80523-1370, United States S Supporting Information *

ABSTRACT: Using fast lattice Monte Carlo (FLMC) simulations (Wang, Q. Sof t Matter 2009, 5, 4564) and the corresponding polymer lattice field theories, including the lattice self-consistent field and Gaussian-fluctuation (LGF) theories, we studied a model system of incompressible homopolymer melts on a hexagonal lattice, where each lattice site is occupied by a total of ρ0 ≥ 1 polymer segments. We generalized the cooperative motion algorithm (Pakula, T. Macromolecules 1987, 20, 679), as well as the related vacancy diffusion algorithm (Reiter, J.; Edling, T.; Pakula, T. J. Chem. Phys. 1990, 93, 837), originally proposed for the self- and mutual-avoiding walk (where ρ0 = 1) to the case of ρ0 > 1, where our generalized algorithm is highly efficient (i.e., nearly rejection-free). On the other hand, we extended the method of Wang (Wang, Z.-G. Macromolecules 1995, 28, 570) to calculate various single-chain properties in LGF theory. Direct comparisons between FLMC and LGF results, both of which are based on the same Hamiltonian (thus without any parameter-fitting between them), unambiguously and quantitatively reveal the effects of non-Gaussian fluctuations neglected by the latter. We found that FLMC results approach LGF predictions with increasing ρ0, and that the leading order of non-Gaussian fluctuation effects on the singlechain properties is inversely proportional to ρ20. Our work suggests that theories capturing the first-order non-Gaussian fluctuation effects may give quantitative agreement with FLMC simulations of incompressible homopolymer melts at ρ0 ≥ 2 in two and three dimensions. the related vacancy diffusion algorithm8) to two-dimensional (2D) and three-dimensional (3D) lattice models with multiple occupancy of lattice sites where ρ0 > 15 and focus on incompressible (and monodisperse) homopolymer melts in 2D, where the system fluctuations/correlations are expected to be stronger than in 3D. Polymers in 2D (such as DNA molecules adsorbed to charged lipid bilayers, ultrathin polymer layers at the interface of two fluids or in lubricated contacts between the solids, etc.) are fascinating systems and of great importance in biology, tribology, material processing, and membrane physics. Their applications include surface and interface modification, colloid stabilization, polymer adhesion and lubrication, polymer intercalation, and biopolymers at cell membranes. Homopolymer melts in 2D have been studied by many groups using theories9−18 and simulations.8,19−49 Most of these studies, however, use hard excluded-volume interactions (e.g., the Lennard-Jones potential or SMAW)8,9,11,12,14−16,18−49 and thus compressible melts or incompressible solutions (where the overall polymer volume fraction ϕ < 1). To the best of our knowledge, only chain dimensions (i.e., the chain mean-square end-to-end distance and radius of gyration) were briefly

1. INTRODUCTION This is the fifth paper in our series of study,1−4 where fast lattice Monte Carlo (FLMC) simulations 5 are directly compared with the corresponding polymer lattice field theories based on the same model system (Hamiltonian), thus without any parameter-fitting, to unambiguously and quantitatively reveal the effects of fluctuations and correlations either neglected or treated approximately in the theories. In the first paper in this series (referred to as Paper I1 hereafter), we studied a model system of compressible homopolymer melts (or equivalently, homopolymers in an implicit, good solvent). Here we study incompressible melts, a model system commonly used in polymer field theories.6 This requires development of new lattice simulation techniques not used and leads to new results not found in our previous work.1−4 Due to the small compressibility of real polymer melts, they are usually modeled as incompressible to reduce the number of parameters. The incompressibility constraint requires the total number of polymer segments at any spatial position r be constant [i.e., ρ̂(r) = ρ0]. While this constraint is commonly enforced in polymer field theories,6 in particle-based simulations it cannot be used in continuum but only on a lattice. Conventional lattice models correspond to ρ0 = 1 (i.e., the selfand mutual-avoiding walk, SMAW), and the cooperative motion algorithm (CMA) proposed by Pakula7 is the only Monte Carlo (MC) method to simulate such incompressible and monodisperse melts. In this work, we generalize CMA (and © 2014 American Chemical Society

Received: July 23, 2014 Revised: September 18, 2014 Published: September 18, 2014 12059

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

Article

BFM2, and BFM3 lattices in 3D. 52,53 For an initial configuration satisfying the incompressibility constraint, our generalized CMA proceeds as follows: (1) Randomly choose a lattice site r0 and a segment, the position of which along its chain contour is denoted by s, from the set of “available” segments located at the nearest neighbors of r0; an inner segment is available if the number of broken bonds is less than 2 when it is moved to r0, and an end segment is always available. The trial move is unsuccessful if there is no available segment; otherwise, we proceed to the next step. (2) Move the chosen segment s to r0 (referred to as “hopping”). If this hopping breaks its bond, for example, with segment s + 1, also reptate the segments after s along the chain contour in the direction from N to 1 until a kink or an end segment t is encountered [that is, segment s′ + 1 moves to the old position of segment s′ for all s′ ∈ (s,t)]; similar reptation but in the opposite direction is performed if the hopping of segment s breaks its bond with segment s − 1. Let r1 be the old position of segment s if no reptation is performed or that of segment t otherwise. If r1 = r0, the trial move is completed; otherwise, we have ρ0 + 1 segments at r0 and ρ0 − 1 segments at r1 (referred to as the “vacancy”) and proceed to the next step. (3) Randomly choose a segment from the set of available segments located at the nearest neighbors of r1 (which excludes those leading to the reverse of the previous movement; again, the trial move is unsuccessful if there is no available segment), and move it to r1 in a similar way to step (2) above. This makes the vacancy diffuse in the system (and updates r1). This step is repeated until r1 coincides with r0 so that the incompressibility constraint is satisfied again. Steps (1)−(3) then form a trial move of our generalized CMA; in Section 1 of the Supporting Information, we generalize the related vacancy diffusion algorithm proposed by Reiter et al. for conventional lattice models8 to the ρ0 > 1 case. For homopolymer melts studied here, unsuccessful trial moves are rejected and successful trial moves are always accepted. We find that the fraction f of successful trial moves increases monotonically toward 1 with increasing ρ0; for example, with N = 40 and L = 20, we have f = 0.890204, 0.996207, and 0.999997, respectively, at ρ0 = 1, 2, and 4. This indicates that our generalized CMA is highly efficient (i.e., nearly rejection-free) for ρ0 > 1. We prepare the initial configurations by consecutively arranging all chains along the y direction as straight rods; when a chain reaches the boundaries of the simulation box in the y direction, we fold it back into the adjacent lattice layer. We use (1−4) × 107 Monte Carlo steps (MCS) in each simulation run, where one MCS is defined as on average one segment move for each segment in the system. We find that the average number of segment moves, ns, in each CMA trial move increases with increasing ρ0; for example, with N = 40 and L = 20, we have ns ≈ 384, 800, 1003, 1134, 1212, and 1255, respectively, at ρ0 = 1, 2, 4, 8, 16, and 32. We monitor the variation of mean-square chain end-to-end distance with the number of MCS to ensure the equilibration of our simulations. We also estimate the error bar of each ensemble-averaged quantity (described below) as three times its standard deviation, with the statistical correlation among samples collected after equilibration taken into account through the correlation function method.54 To quantify the fluctuation and correlation effects in our incompressible homopolymer melts, we calculate the bond− bond correlation function Cb, the mean-square end-to-end

reported in previous MC simulations of incompressible melts (i.e., ϕ = 1),8,33,43 where SMAW was used. The mean-field limit, on the other hand, corresponds to ρ0 → ∞, where neither fluctuations nor correlations present, the chain conformations are ideal, and the system is structureless. The lattice Gaussian-fluctuation (LGF, or equivalently the firstorder perturbation50) theory quantitatively describes the system at large ρ0 where only Gaussian fluctuations present. By direct comparisons between FLMC and LGF results at various ρ0, we therefore unambiguously quantify non-Gaussian fluctuation effects on various properties of 2D incompressible homopolymer melts. Our work suggests that theories capturing the firstorder non-Gaussian fluctuation effects may give quantitative agreement with FLMC simulations of incompressible homopolymer melts at ρ0 ≥ 2 in 2D (thus also in 3D).

2. MODEL AND METHODS 2.1. Model System. We consider a system of incompressible homopolymer melts of n chains each having N segments on a hexagonal (HEX) lattice of totally V = Ld sites, where L is the total number of lattice layers in each direction and d = 2 the number of dimensions. Each polymer segment occupies one lattice site, and each lattice site is occupied by a total of ρ0 = nN/V polymer segments, which is chosen to be a positive integer. Periodic boundary conditions are applied in all directions. The spatial coordinates of lattice site (i,j) are given by x = (√3i/2)b and y = (i/2+j)b, where {i,j} = 0,...,L−1, and b is the lattice spacing. The configuration integral of the system is n

A(n , N , V ) =

N

n

∏ ∏ ∑ ·exp(−β ∑ hkC) ∏ δ ρ ̂(r),ρ

0

k = 1 s = 1 R k ,s

k=1

r

(1) th

where Rk,s denotes the spatial position of the s segment on the kth chain, and ρ̂(r) ≡ Σkn = 1ΣsN= 1δr,Rk,s is the microscopic number density of polymer segments at lattice site r. In addition, hCk is the Hamiltonian of the kth chain due to its connectivity; βhCk = 0 if the chain connectivity is maintained and ∞ otherwise, where β ≡ 1/kBT with kB being the Boltzmann constant and T the thermodynamic temperature. Finally, the last term in eq 1 enforces the incompressibility constraint with δ denoting the Kronecker δ function. Clearly, this system has no interesting thermodynamic properties, and the incompressibility constraint makes the intra- and interchain correlations trivially dependent on each other, thus providing a simple but useful model for polymer melts. Also, the system fluctuations/correlations are mainly controlled by the number of chains per volume of Rde , ρ0Rde /N, where Re denotes the root-mean-square chain end-to-end distance. Since Re ∼ N1/2 for large N in both 2D9 and 3D51 melts, we expect the effects of N to be rather small in 2D (in contrast to the 3D case) and therefore set N = 40 in this work. 2.2. Fast Lattice Monte Carlo (FLMC) Simulations with Cooperative Motion Algorithm (CMA). Our FLMC simulations are performed in a canonical ensemble, where n and ρ0 (thus V) are fixed. For trial moves, we generalize CMA originally proposed by Pakula for conventional lattice models (i.e., ρ0 = 1)7 to the ρ0 > 1 case. We define an inner segment s (s ≠ 1 or N) at lattice position Rs as a kink segment when Rs+1 − Rs−1 belongs to the set of allowed bond vectors on the lattice; this allows our generalized CMA to be applicable to HEX and STAR lattices in 2D, as well as the face-centered cubic (FCC), 12060

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

Article

Now consider the fluctuations around ω*(r). Under the random-phase approximation,9 we have

distance R2e (Ns) of a partial chain having Ns segments, and the single-chain structure factor S1 from each collected sample after equilibration. In particular, N n(N − Ns − 1)

C b(u) ≡

n

N − Ns − 1

∑ ∑ k=1

ln Q [ω] ≈ −

bk , s ·bk , s + Ns b

s=1

2

(2)

1 ≡ n(N − Ns + 1)

n k=1

ALGF = .n ∏ [2πρc N 2P(k)]−1/2

2

and

=

1 nN 1 nN

n

N

n

N

N

∑k = 1 ∑s = 1 ∑t = 1 exp[− − 1 k· (R k , s − R k , t)] ∑k = 1 ∑s = 1 exp( − −1 k·R k , s)|2 (4)

⟨A⟩ ≡

where k denotes the wave vector. For an ideal chain, Cb,0(u) = 0, R2e,0(Ns) = (Ns −1)b2, and N

S1,0(k) =

N − B(k)[2 + NB(k) − 2B (k)] N[1 − B(k)]2

1 n

n

∑ Ak

=−

k=1

A(η) = .n(η) ∏

∑ exp(−

r

−1 k·b)/z

where b is an allowed bond vector and z the total number of allowed bond vectors.55 For HEX lattice, we have z = 6 and BHEX(k) = [cos k2 + 2 cos(√3k1/2) cos(k2/2)]/3 with k1 and k2 being the x and y components of the wave vector k, respectively; note that the L2 discrete values of k also form a HEX lattice with the coordinates of lattice site (i,j) in k-space given by k1 = 4π(i − j/2)/√3L and k2 = 2πj/L. Since the anisotropy of HEX lattice is rather small,52 we only examine S1(k) ≡ Σ|k|=k S1(k)/Nk, where the summation is over all k satisfying |k| = k and Nk is the total number of such k. 2.3. Lattice Gaussian-Fluctuation (LGF) Calculations. 2.3.1. LGF Theory. Here we start from the same model system as in Sec. 2.1 and insert the identity Πrδρ̂(r),ρ0 = Πr∫ [dω(r)/2π] · exp{−Σriω(r)[ρ̂(r) − ρ0]} into eq 1, where ω(r) is the conjugate field enforcing the incompressibility constraint and i ≡ √−1. This gives

r

∫

⎡ ⎤ dω(r) ·exp⎢ρ0 ∑ iω(r) + n ln Q ⎥ ⎢⎣ ⎥⎦ 2π r

where the single-chain partition N N Q ≡ ∏s = 1 ∑R · exp[−βhC − i ∑s = 1 ω(R s)]/. , s

(10)

η=0

⎡

∫ dω2π(r) ·exp⎢⎢ρ0 ∑ iω(r) ⎣

r

⎤ + n ln Q (η)⎥ ⎥⎦

(6)

{b}

A = .n ∏

1 ∂ ln A(η) ∂η n

where (5)

with B(k) ≡

(9)

where ρc = ρ0/N denotes the chain number density. 2.3.2. General Method of Calculating Single-Chain Properties in LGF Theory. In order to calculate the ensemble average of a single-chain property A, we extend the method of Wang;56 that is, we introduce an auxiliary parameter η and define βhk C [η] ≡ βhk C + ηAk , where Ak denotes the property A of the kth chain. We then have

⟨(R k , s + Ns − 1 − R k , s) ⟩

s=1

(3)

S1(k) =

(8)

k≠0

N − Ns + 1

∑ ∑

∑ P(k)δω̂ (k)δω̂ (−k) k≠0

where δω(r) ≡ ω(r) − ω*(r) and P(k) = S1,0(k)/N. Substituting the above into eq 7 and evaluating the Gaussian integral over δω̂ (k), we can approximate A by

where bk,s ≡ Rk,s+1 − Rk,s is the s bond of the kth chain, u ≡ Ns/ N with Ns being the number of segments between bonds bk,s and bk,s+Ns, and b the effective bond length (equal to the lattice spacing for HEX lattice); th

Re2(Ns)

N2 2V 2

(11)

with Q (η) ≡ ∏sN= 1 ∑R s · exp⎡⎣−βhC − ηA − i ∑sN= 1 ω(R s)⎤⎦ .(η) N

and .(η) ≡ ∏s = 1 ∑R · exp( −βhC − ηA). s LGF theory then gives ALGF(η) = .n(η) ∏ [2πρc N 2P(k, η)]−1/2

(12)

k≠0

where P(k,η) is the Fourier transform of {[V/N2Q(η)][δ2Q(η)/δiω(r)δiω(r′)]}|ω*, and ⟨A⟩LGF = ⟨A⟩0 +

1 2n

∑ k≠0

1 ∂P(k, η) P(k) ∂η

(13)

η=0

with ⟨⟩0 denoting the ensemble average in an ideal-chain system. In the thermodynamic limit (i.e., V → ∞ at constant ρ0), the summation over k is replaced by an integral (i.e., Σk≠0 → [V/(2π)d] ∫ dk) where the integration limits in each dimension are from 0 to 2π, and we have

(7)

function with

⟨A⟩LGF = ⟨A⟩0 +

N

. ≡ ∏s = 1 ∑R · exp( −βhC ). The self-consistent field solution s of this homogeneous system is given by ω*(r) = 0 (which can be shifted by an arbitrary constant) and ρ*(r) = ρ0.

1 2ρc

∫ (2dπk)d P(1k) ∂P(∂kη, η)

η= 0

(14)

In eqs 13 and 14,

12061

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

∂P(k, η) ∂η

η=0

Article

1 = 2 N

=

1 N2

N

C ⎡ N−1 ⎤ ∂ ⎢ ∏s ′= 1 ∑ bs′ ·exp[−βh − ηA − i k·(R t − R s)] ⎥ N−1 ⎥ ∂η ⎢⎣ ∏s ′= 1 ∑ b ·exp(− βhC − ηA) ⎦

N

∑∑ s=1 t=1 N

s′

∑ ∑ {⟨A⟩0 ⟨exp[−i k·(R t − R s)]⟩0 − ⟨Aexp[−i k·(R t − R s)]⟩0 } (15)

s=1 t=1

where bs ≡ Rs+1 − Rs is the sth bond on the chain. Note that eq 14 gives the same results as the first-order perturbation theory (for the same model).50 In the following, we derive [∂P(k,η)/ ∂η]|η=0 for the calculation of various single-chain properties in LGF theory. 2.3.3. Bond−Bond Correlation Function. We first consider Cb(e1, e2 − 1) ≡ ⟨be1 · be2−1⟩/b2 with 1 ≤ e1 < e2 ≤ N. According to the above, for A = be1 · be2−1/b2, we have ∂P(k, η) ∂η

η= 0

and Σmd = 1[∂2BHEX(k)/∂k2m]{kn≠m} = −BHEX(k). Substituting eq 18 into eq 13, we obtain R2LGF(e1,e2) with R20(e1,e2) = (e2 − e1)b2, which is then averaged along the chain contour to obtain R2e,LGF(Ns) as Re2,LGF(Ns) =

(16)

C b,LGF(s , s + Ns)

s=1

(17)

∂P(k′, η) ∂η

In Section 2 of the Supporting Information, we derive the Gaussian-fluctuation (GF) prediction Cb,GF(u) for incompressible melts of the discrete and continuous Gaussian chains in continuum. 2.3.4. Mean-Square Partial Chain End-to-End Distance. We first consider the mean-square distance between segments e1 and e2, R2(e1, e2) ≡ ⟨(Re2 − Re1)2⟩, with 1 ≤ e1 < e2 ≤ N. For A = (Re2 − Re1)2, we have ∂P(k, η) ∂η

η= 0

+

e1− 1

+

+

+

+

∑ s = e1

k ′≠ 0

1 ∂P(k′, η) P(k′) ∂η

η=0

⎧ e − 1 e2 2 ⎪1 ⎨ ∑ ∑ [Bt − s (k′)Be2 − e1(k) 2⎪ N ⎩ s = 1 t = e1+ 1

N

[Bt − s (k′)Be2 − e1(k) − Be1− s + t − e2 (k′)Be2 − e1(k′ + k)]

e2

[Bt − s (k′)Be2 − e1(k) − B s − e1+ e2 − t (k)Bt − s (k′ + k)]

s = e1 t = s + 1

+

N

∑ ∑

[Bt − s (k′)Be2 − e1(k) − Bt − e2 (k′)B s − e1(k)Be2 − s (k′

s = e1 t = e 2 + 1

⎫ ⎪ + k)]⎬ ⎪ ⎭

P0(k, t − s , e 2 − e1)

(22)

Substituting eq 22 into eq 21, we obtain ω̂ LGF(e1, e2; k) with ω̂ 0(e1, e2; k) = Be2−e1(k), from which S1,LGF(k) is calculated as

P0(k, t − s , t − s)

s = e1 t = s + 1 e2 − 1

η= 0

∑ ∑ e2 − 1

e2

∑ ∑

∑

s = 1 t = e2 + 1

s = 1 t = e2 + 1 e2 − 1

=

∑ ∑ e2 − 1

N

∑ ∑

1 2n

− Be1− s (k′)Be2 − t (k)Bt − e1(k′ + k)]

⎡ e − 1 e2 2 1 = 2 ⎢ ∑ ∑ P0(k, t − s , t − e1) N ⎢⎣ s = 1 t = e1+ 1 e1− 1

s=1

We first consider the correlation function ω̂ (e1, e2; k) ≡ ⟨exp[ik · (Re2 − Re1)]⟩ between two segments e1 and e2; since ω̂ (e1, e2; k) = ω̂ (e2, e1; k) and ω̂ (e1, e1; k) = 1, we only consider 1 ≤ e1 < e2 ≤ N. For A(k) = exp[ik · (Re2 − Re1)], we have

N − Ns − 1

∑

2 RLGF (s , s + Ns − 1)

(21)

Substituting eq 16 into eq 13, we obtain Cb,LGF(e1, e2 − 1) with Cb,0(e1,e2 − 1) = 0, which is then averaged along the chain contour to obtain Cb,LGF(u) as N N − Ns − 1

∑

as ⟨A(k)⟩LGF = ⟨A(k)⟩0 +

n ≠ m}

N − Ns + 1

2.3.5. Single-Chain Structure Factor. Here we rewrite eq 13

d

C b,LGF(u) =

1 N − Ns + 1

(20)

2[Be1(k) − 1][B N + 1(k) − Be2 (k)] = b2N 2Be1+ 2 (k)[1 − B(k)]2 ⎛ ∂B(k) ⎞2 × ∑⎜ ⎟ m = 1 ⎝ ∂k m ⎠{k

η=0

N

⎤ ∑ P0(k, t − s , e2 − s)⎥⎥ t = e2 + 1 ⎦ N

S1,LGF(k) = (18)

⎛

d

∑ ⎜ ∂ B(2k) ⎟

m = 1 ⎝ ∂k m ⎠{k d

+ t(t − 1)B s − 2 (k)

⎞

2

⎛

⎞

∑ ⎜ ∂B(k) ⎟ m=1

n ≠ m}

n ≠ m}

N

s = 1 s ′= 1

(23)

3. RESULTS AND DISCUSSION With the chain length N = 40, the total number of polymer segments, ρ0, at each lattice site is the only parameter in our model system of incompressible homopolymer melts in the thermodynamic limit, where the system size L → ∞. When ρ0 = 1, our system recovers the self- and mutual-avoiding walk

2

⎝ ∂km ⎠{k

N

∑ ∑ ω̂ LGF(s , s′; k)

Finally, we calculate S1,LGF(k) ≡ Σ|k|=k S1,LGF(k)/Nk as in FLMC simulations.

where P0(k, s , t ) = tb2B s (k) + tB s − 1(k)

1 N

(19) 12062

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

Article

(SMAW). For finite ρ0 > 1, it is very close to that studied by Semenov and Johner,17 except that the latter is in continuum and weakly compressible [i.e., |ρ̂(r) − ρ0| ≪ ρ0, where ρ̂(r) denotes the total number of segments at spatial position r]. In the limit of ρ0 → ∞, chains adopt the ideal conformations and the self-consistent field (SCF) theory becomes exact. In the following, we first examine the effects of finite L, then present the effects of ρ0 on the bond−bond correlation function Cb(u) (defined in eq 2), the mean-square end-to-end distance R2e (Ns) of a partial chain having Ns segments, and the single-chain structure factor S1(k) with k being the wave vector length, obtained from both FLMC simulations and LGF calculations. Since LGF theory only includes Gaussian fluctuations, the differences between FLMC and LGF results unambiguously quantify the effects of non-Gaussian fluctuations. 3.1. Finite-Size Effects. Finite-size effects (FSE) can be present in both LGF calculations and FLMC simulations. As shown in Sec. 3.2 below, ρ0Cb(u) is the largest at u = 1/N = 0.025; we therefore expect that it exhibits the largest FSE at u = 0.025. Similarly, we use ρ0R2e (N) and ρ0S1(k*Re,0 = 4.53) to examine FSE, where Re,0 = (N − 1)1/2b with b denoting the lattice spacing. Figure 1 shows FSE in LGF calculations, quantified by |ALGF(L) − ALGF(L → ∞)|, where ALGF(L) at finite L is given

exception is Figure 4a below, where LGF calculations are performed in the thermodynamic limit (i.e., L → ∞). 3.2. Bond−Bond Correlation Function. Since chains are ideal in the mean-field limit, SCF theory gives Cb(u) = 0. As shown in eq 14, Gaussian-fluctuation (GF) theory then predicts Cb,GF(u) to be inversely proportional to ρc. Figure 2a therefore

Figure 1. Logarithmic plot of the finite-size effects on A = ρ0 [R2e (N)/ R2e,0 − 1], ρ0[S1,0(k*) − S1(k*)] with k*Re,0 = 4.53, and ρ0Cb(u = 1/N) in LGF calculations. N = 40 is used here.

by eq 13 and ALGF(L → ∞) by eq 14; here, A is either ρ0[R2e (N)/R2e,0−1], ρ0Cb(1/N), or ρ0[S1,0(k*) − S1(k*)]. We see that R2e exhibits the largest FSE, but in any case |ALGF(L) − ALGF(L → ∞)| ≲ 0.003 for L ≥ 20. To examine FSE in our simulations, we use ρ0 = 1, where the largest FSE are expected. Table 1 lists AFLMC(L) for L = 20 and 40, where we see that for all the quantities the difference between these two system sizes is within our sampling error. Since this is also found for |ALGF(L ≥ 20) − ALGF(L → ∞)| given above, FSE in both FLMC simulations and LGF calculations are negligible when L ≥ 20. Hereafter we compare FLMC results with LGF predictions both at L = 20; the only

Figure 2. (a) Semilogarithmic plot of the bond−bond correlation function Cb,GF(u) in the thermodynamic limit predicted by Gaussian fluctuation (GF) theory for incompressible melts of chains on the hexagonal (HEX) lattice, discrete Gaussian chains (DGC), continuous Gaussian chains (CGC), and CGC with Padé approximation, where the inset shows the logarithmic plot. (b) Cb(u) obtained from FLMC simulations at various ρ0 and LGF calculations, where the inset shows the logarithmic plot. (c) Logarithmic plot of the deviation of Cb(u) in FLMC simulations from its LGF prediction, ΔCb defined in eq 24, where the straight line has a slope of −1. N = 40 and L = 20 are used in both parts (b) and (c).

Table 1. Finite-Size Effects on R2e (N)/R2e,0 −1, Cb(1/N), and S1,0(k*) − S1(k*) with k*Re,0 = 4.53 in Our Simulationsa R2e (N)/R2e,0 L = 20 L = 40 a

−1

0.627 ± 0.006 0.626 ± 0.005

Cb(1/N)

S1,0(k*) − S1(k*)

7.597 ± 0.004 7.593 ± 0.005

6.623 ± 0.004 6.624 ± 0.004

shows ρcR2e,0Cb,GF(u) for incompressible homopolymer melts of various chain models, including chains on the hexagonal (HEX) lattice, and continuous Gaussian chains (CGC) and discrete Gaussian chains (DGC) in continuum. In addition, an analytical expression, eq 14 in the Supporting Information, is obtained for CGC model under Padé approximation57 and also

N = 40 and ρ0 = 1 are used here. 12063

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

Article

shown in Figure 2a. In all cases, Cb,GF(u) decreases monotonically with increasing u as expected. Furthermore, at u = 0.025, HEX model has the largest Cb,GF, while the DGC model has the smallest, and Padé approximation gives slightly larger Cb,GF than CGC model. The differences are solely due to the different chain models and become smaller with increasing u, indicating the universal behavior of Cb,GF at large u. On the other hand, eq 14 in the Supporting Information gives Cb,GF(u) ∝ 1/u at small u for CGC model under Padé approximation, and the same behavior is also found for the CGC model (without Padé approximation) as shown in the inset of Figure 2a. Figure 2b compares ρ0 C b (u) obtained from FLMC simulations with LGF prediction (i.e., ρ0Cb,GF(u) on HEX lattice). In FLMC simulations, Cb(u) exhibits anticorrelations (i.e., Cb < 0) for 0.2 ≤ u ≤ 0.75 at ρ0 = 1 (i.e., in SMAW model); similar results are found by Meyer et al.44,45 and Wen et al.58 When ρ0 ≥ 2, Cb decreases monotonically with increasing u and no anticorrelations are found. Figure 2b also shows that ρ0Cb from FLMC simulations approaches LGF prediction with increasing ρ0. To quantify non-Gaussian fluctuation effects on Cb, we define N−2

ΔC b ≡

∑uN = 1 [C b,FLMC(u) − C b,LGF(u)]2 N−2

(24)

and examine its dependence on ρ0. Figure 2c shows ρ0ΔCb ≈ C1ρ−1 0 with C1 = 0.47 ± 0.01 obtained from the unweighted linear least-squares fitting. This can be understood as follows: in general, a physical quantity A at finite ρ0 can be expanded as a power series of 1/ρ0, i.e., A(ρ0 ) = ASCF + A1/ρ0 + A 2 /ρ02 + ···

Figure 3. (a) Difference between the mean-square end-to-end distance of a partial chain having Ns segments, R2e (Ns), and its mean-field (SCF) prediction, R2e,0(Ns) = Ns −1, obtained from FLMC simulations at various ρ0 and LGF theory. (b) Logarithmic plot of the deviation of R2e (Ns) in FLMC simulations from its LGF prediction, ΔR2e defined in eq 26, where the straight line has a slope of −1. N = 40 and L = 20 are used here.

(25)

where the zeroth-order term ASCF is the SCF prediction, the first-order term A1/ρ0 is the GF contribution which can be obtained from GF theory, and the rest is non-Gaussian fluctuation contribution which, to the leading order, is inversely proportional to ρ20; we therefore have non-Gaussian fluctuation effects revealed approximately as ρ0(A − ALGF) ∝ ρ−1 0 . 3.3. Mean-Square Partial Chain End-to-End Distance. SCF theory gives the ideal-chain result R2e,0(Ns) = (Ns − 1)b2. Figure 3a shows ρ0[R2e (Ns) − R2e,0(Ns)] obtained from LGF theory and FLMC simulations. We see that SCF theory underestimates the chain dimensions. On the other hand, LGF theory predicts that ρ0[R2e (Ns) − R2e,0(Ns)] is independent of ρ0 and increases monotonically with increasing Ns. In FLMC simulations, however, ρ0[R2e (Ns) − R2e,0(Ns)] exhibits a small maximum at Ns = 39 when ρ0 = 1 (i.e., in SMAW model). In addition, when ρ0 ≥ 2, FLMC results are always larger than LGF prediction and approach the latter monotonically with increasing ρ0; at ρ0 = 1, however, ρ0[R2e (Ns) − R2e,0(Ns)] obtained from FLMC simulations is larger than LGF prediction for 4 ≤ Ns ≤ 25 while the opposite occurs for Ns = 3 and Ns ≥ 26. To quantify non-Gaussian fluctuation effects on R2e , similar to ΔCb we define

of non-Gaussian fluctuation effects is inversely proportional to ρ20. 3.4. Single-Chain Structure Factor. SCF theory gives the ideal-chain result S1,0(k) (see eq 5), shown in the inset of Figure 4a for various chain models. The differences at large k are solely due to the different chain models. For example, the oscillation of S1,0(k) on HEX lattice is due to the spatial discretization of the lattice, and the slightly larger value of DGC model than CGC is due to the finite chain discretization. On the other hand, at small kRe,0 S1,0(k) of the different chain models collapse together, indicating their universal behavior. In addition, Padé approximation57 introduces deviation from CGC model at intermediate k with the largest difference, SCGC 1,0 (k)/ N−SPadé (k)/N ≈ 0.072, occurring at kR ≈ 2.343. 1,0 e,0 Figure 4a shows ρ0[S1,0(k) − S1(k)] obtained from LGF theory and FLMC simulations. We see that SCF theory overestimates S1(k) for all k except at k = 0 [note that S1(k = 0) = N]. On the other hand, LGF theory predicts that ρ0[S1,0(k) − S1(k)] is independent of ρ0 and exhibits a maximum at kRe,0 ≈ 3.645. A maximum is also found in FLMC simulations, where kRe,0 can only take discrete values commensurate with L = 20. In addition, we see that LGF theory overestimates ρ0[S1,0(k)− S1(k)] at all these kRe,0 values, except at the smallest kRe,0 ≈ 2.265 when ρ0 ≥ 2, and that FLMC results approach LGF prediction with increasing ρ0. To quantify non-Gaussian fluctuation effects on S1, we define

N

ΔR e2

≡

2 2 ∑ N = 3 [R e,FLMC (Ns) − R e,LGF (Ns)]2 s

N−2

(26)

Figure 3b shows ρ0ΔR2e ≈ C2ρ−1 0 with C2 = 6.1 ± 0.2 obtained from the unweighted linear least-squares fitting, where the data point at ρ0 = 1 is excluded. We again see that the leading order 12064

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

ΔS1 ≡

∑kR

e,0 > 0

Article

tions neglected by the latter. This is the fifth paper in our series of studies1−4 with similar purpose. Our FLMC simulations are performed in a canonical ensemble. We have generalized the cooperative motion algorithm (CMA),7 as well as the related vacancy diffusion algorithm,8 originally proposed for the self- and mutualavoiding walk (SMAW, where ρ0 = 1) to the case of multiple occupancy of lattice sites (MOLS, where ρ0 > 1). For incompressible homopolymer melts with MOLS, our generalized CMA is highly efficient (i.e., nearly rejection-free) and gives the accurate simulation results needed for quantitatively revealing non-Gaussian fluctuation effects. On the other hand, only single-chain properties are of interest in our simple model of incompressible melts, and we have extended the method of Wang56 to calculate in LGF theory the bond−bond correlation function Cb(u) (defined in eq 2), the mean-square end-to-end distance R2e (Ns) of a partial chain having Ns segments, and the single-chain structure factor S1(k) with k being the wave vector length. With L = 20, where the finite-size effects in both FLMC simulations and LGF calculations are found to be negligible compared to our sampling errors, we find that FLMC results approach LGF predictions with increasing ρ0, and that the leading order of non-Gaussian fluctuation effects is inversely 2 proportional to ρ20; in particular, ΔCb = (0.47 ± 0.01)ρ−2 0 , ΔRe −2 −2 = (6.1 ± 0.2)ρ0 (for ρ0 ≥ 2), and ΔS1 = (1.028 ± 0.006)ρ0 , where ΔCb, ΔR2e , and ΔS1 defined in eqs 24, 26, and 27, respectively, measure the deviation in the corresponding property between FLMC and LGF results. Given the closeness of the line (corresponding to the above relation) to the points (representing FLMC results) shown in Figures 2c, 3b, and 4b, our work suggests that theories capturing the first-order nonGaussian fluctuation effects may give quantitative agreement with FLMC simulations of incompressible homopolymer melts at ρ0 ≥ 2 in 2D and 3D. This work extends FLMC simulations to incompressible polymer melts and paves the way to FLMC simulations of more interesting and complicated systems such as incompressible polymer blends and block copolymers commonly studied with polymer field theories. It also provides accurate simulation data for quantitative test of advanced theories capturing nonGaussian fluctuation effects in polymeric systems and promotes the development of such theories.

[S1,FLMC(k) − S1,LGF(k)]2 Mk

(27)

where the summation is over all the discrete k values commensurate with L = 20 and Mk the total number of such values. Figure 4b shows ρ0ΔS1 ≈ C3ρ−1 0 with C3 = 1.028 ±

Figure 4. (a) Logarithmic plot of the difference between the singlechain structure factor S1(k) in incompressible melts and its mean-field (SCF) prediction S1,0(k) obtained from FLMC simulations at L = 20 and various ρ0, as well as from LGF calculations in the thermodynamic limit, where the inset shows the logarithmic plot of S1,0(k) for chains on the hexagonal (HEX) lattice, discrete Gaussian chains (DGC), continuous Gaussian chains (CGC), and CGC with Padé approximation. (b) Logarithmic plot of the deviation of S1(k) in FLMC simulations from its LGF prediction, ΔS1 defined in eq 27, where the straight line has a slope of −1, and N = 40 and L = 20 are used.

■

ASSOCIATED CONTENT

S Supporting Information *

0.006; the leading order of non-Gaussian fluctuation effects on S1 is therefore inversely proportional to ρ20, consistent with our ΔCb and ΔR2e data.

Generalized vacancy diffusion algorithm for incompressible homopolymer solution, as well as the Gaussian-fluctuation predictions of the bond−bond correlation function for incompressible melts of discrete and continuous Gaussian chains in continuum. This material is available free of charge via the Internet at http://pubs.acs.org.

4. CONCLUSIONS Using fast lattice Monte Carlo (FLMC) simulations5 and the corresponding polymer lattice field theories, including the lattice self-consistent field (LSCF) and Gaussian-fluctuation (LGF) theories, we have studied a model system of incompressible homopolymer melts on a hexagonal lattice. This system has a few parameters: the chain length N = 40, the system size L, and the number of polymer segments at each lattice site ρ0 controlling the system fluctuations and correlations. Direct comparisons between FLMC and LGF results, both of which are based on the same Hamiltonian (thus without any parameter-fitting between them), unambiguously and quantitatively reveal the effects of non-Gaussian fluctua-

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by NSF CAREER Award CBET0847016, which is gratefully acknowledged. 12065

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

■

Article

(27) Dickman, R.; Hall, C. Equation of State for Athermal Lattice Chains. J. Chem. Phys. 1986, 85, 3023. (28) Dickman, R.; Hall, C. Equation of State for Chain Molecules: Continuous-Space Analog of Flory Theory. J. Chem. Phys. 1986, 85, 4108. (29) Dickman, R. New Simulation Method for the Equation of State of Lattice Chains. J. Chem. Phys. 1987, 87, 2246. (30) Dickman, R. On the Equation of State of Athermal Lattice Chains: Test of Mean-Field and Scaling Theories in Two Dimensions. J. Chem. Phys. 1989, 91, 454. (31) Carmesin, I.; Kremer, K. Static and Dynamic Properties of TwoDimensional Polymer Melts. J. Phys. (Paris) 1990, 51, 915. (32) Dickman, R. Equation of State of Two Dimensional Lattice Chains at the Theta Point. J. Chem. Phys. 1992, 96, 1516. (33) Gauger, A.; Weyersberg, A.; Pakula, T. Monte Carlo Studies of Static Properties of Interacting Lattice Polymers with the CooperativeMotion Algorithm. Makromol. Chem., Theory Simul. 1993, 2, 531. (34) Reiter, J.; Duering, E. Monte Carlo Simulation of Freely Jointed Chains with Variable Bond Lengths. Macromol. Theory Simul. 1995, 4, 667. (35) Nelson, P.; Hatton, T.; Rutledge, G. General Reptation and Scaling of 2d Athermal Polymers on Close-Packed Lattices. J. Chem. Phys. 1997, 107, 1269. (36) Ostrovsky, B.; Smith, M.; Bar-Yam, Y. Simulations of Polymer Interpenetration in 2D Melts. Int. J. Mod. Phys. C 1997, 8, 931. (37) Teraoka, I.; Wang, Y. Crossover from Two- to ThreeDimensional Contraction of Polymer Chains in Semidilute Solutions Confined to a Narrow Slit. Macromolecules 2000, 33, 6901. (38) Polanowski, P.; Pakula, T. Studies of Polymer Conformation and Dynamics in Two Dimensions using Simulations Based on the Dynamic Lattice Liquid (DLL) Model. J. Chem. Phys. 2002, 117, 4022. (39) Balabaev, N.; Darinskii, A.; Neelov, I.; Lukasheva, N.; Emri, I. Molecular Dynamics Simulation of a Two-Dimensional Polymer Melt. Polym. Sci., Ser. A 2002, 44, 781. (40) Yethiraj, A. Computer Simulation Study of Two-Dimensional Polymer Solutions. Macromolecules 2003, 36, 5854. (41) Hehmeyer, O.; Arya, G.; Panagiotopoulos, A. Phase Transitions of Confined Lattice Homopolymers. J. Phys. Chem. B 2004, 108, 6809. (42) Cavallo, A.; Müller, M.; Wittmer, J.; Johner, A.; Binder, K. Single Chain Structure in Thin Polymer Films: Corrections to Flory’s and Silberberg’s Hypotheses. J. Phys.: Condens. Matter 2005, 17, S1697. (43) Polanowski, P.; Jeszka, J. Microphase Separation in TwoDimensional Athermal Polymer Solutions on a Triangular Lattice. Langmuir 2007, 23, 8678. (44) Meyer, H.; Kreer, T.; Aichele, M.; Cavallo, A.; Johner, A.; Baschnagel, J.; Wittmer, J. Perimeter Length and Form Factor in TwoDimensional Polymer Melts. Phys. Rev. E 2009, 79, 050802. (45) Meyer, H.; Wittmer, J.; Kreer, T.; Johner, A.; Baschnagel, J. Static Properties of Polymer Melts in Two Dimensions. J. Chem. Phys. 2010, 132, 184904. (46) Meyer, H.; Schulmann, N.; Zabel, J.; Wittmer, J. The Structure Factor of Dense Two-Dimensional Polymer Solutions. Comput. Phys. Commun. 2011, 182, 1949. (47) Schulmann, N.; Meyer, H.; Wittmer, J.; Johner, A.; Baschnagel, J. Interchain Monomer Contact Probability in Two-Dimensional Polymer Solutions. Macromolecules 2012, 45, 1646. (48) Schulmann, N.; Xu, H.; Meyer, H.; Polińska, P.; Baschnagel, J.; Wittmer, J. Strictly Two-Dimensional Self-Avoiding Walks: Thermodynamic Properties Revisited. Eur. Phys. J. E 2012, 35, 93. (49) Schulmann, N.; Meyer, H.; Kreer, T.; Cavallo, A.; Johner, A.; Baschnagel, J.; Wittmer, J. Strictly Two-Dimensional Self-Avoiding Walks: Density Crossover Scaling. Polym. Sci., Ser. C 2013, 55, 181. (50) Edwards, S. The Size of a Polymer Molecule in a Strong Solution. J. Phys. A: Math. Gen. 1975, 8, 1670. Muthukumar, M.; Edwards, S. Extrapolation Formulas for Polymer Solution Properties. J. Chem. Phys. 1982, 76, 2720. (51) Flory, P. The Configuration of Real Polymer Chains. J. Chem. Phys. 1949, 17, 303.

REFERENCES

(1) Zhang, P.; Zhang, X.; Li, B.; Wang, Q. Quantitative Study of Fluctuation Effects by Fast Lattice Monte Carlo Simulations. I. Compressible Homopolymer Melts. Soft Matter 2011, 7, 4461. (2) Zhang, P.; Li, B.; Wang, Q. Quantitative Study of Fluctuation Effects by Fast Lattice Monte Carlo Simulations. 2. Homopolymer Brushes in an Implicit, Good Solvent. Macromolecules 2011, 44, 7837. (3) Zhang, P.; Li, B.; Wang, Q. Quantitative Study of Fluctuation Effects by Fast Lattice Monte Carlo Simulations. III. Homopolymer Brushes in an Explicit Solvent. Macromolecules 2012, 45, 2537. (4) Zhang, P.; Wang, Q. Quantitative Study of Fluctuation Effects by Fast Lattice Monte Carlo Simulations: Compression of Grafted Homopolymers. J. Chem. Phys. 2014, 140, 044904. (5) Wang, Q. Studying Soft Matter with “Soft” Potentials: Fast Lattice Monte Carlo Simulations and Corresponding Lattice SelfConsistent Field Calculations. Soft Matter 2009, 5, 4564; ibid. 2010, 6, 6206. (6) Fredrickson, G. H. The Equilibrium Theory of Inhomogeneous Polymers; Oxford University Press: Oxford, 2006. (7) Pakula, T. Cooperative Relaxations in Condensed Macromolecular Systems. 1. A Model for Computer Simulation. Macromolecules 1987, 20, 679. (8) Reiter, J.; Edling, T.; Pakula, T. Monte Carlo Simulation of Lattice Models for Macromolecules at High Densities. J. Chem. Phys. 1990, 93, 837. (9) de Gennes, P. G. Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca, 1978. (10) Douglas, J.; Cherayil, B.; Freed, K. Polymers in Two Dimensions: Renormalization Group Study using Three-Parameter Model. Macromolecules 1985, 18, 2455. (11) Duplantier, B. Exact Critical Exponents for Two-Dimensional Dense Polymers. J. Phys. A: Math. Gen. 1986, 19, L1009. (12) Duplantier, B.; Saleur, H. Exact Critical Properties of TwoDimensional Dense Self-Avoiding Walks. Nucl. Phys. B 1987, 290, 291. (13) Petschek, R.; Pfeuty, P. Two-Dimensional Polymer Melts. Phys. Rev. Lett. 1987, 58, 1096. (14) Honnell, K.; Hall, C. A New Equation of State for Athermal Chains. J. Chem. Phys. 1989, 90, 1841. (15) Duplantier, B. Statistical Mechanics of Polymer Networks of Any Topology. J. Stat. Phys. 1989, 54, 581. (16) Kushwaha, K.; Khanna, K. SAFT-D Theory for Hard Sphere and Hard Disc Chain Fluids. Mol. Phys. 1999, 97, 907. (17) Semenov, A.; Johner, A. Theoretical Notes on Dense Polymers in Two Dimensions. Eur. Phys. J. E 2003, 12, 469. (18) Yethiraj, A.; Sung, B.; Lado, F. Integral Equation Theory for Two-Dimensional Polymer Melts. J. Chem. Phys. 2005, 122, 094910. (19) Okamoto, H. Monte Carlo Study of Systems of Linear Oligomers in Two-Dimensional Spaces. J. Chem. Phys. 1976, 64, 2686. (20) Wall, F.; Chin, J.; Mandel, F. Configurations of Macromolecules Subject to Intermolecular Volume Exclusions. J. Chem. Phys. 1977, 66, 3143. (21) Wall, F.; Seitz, W. Simulation of Polymers by Self-Avoiding, Nonintersecting Random Chains at Various Concentrations. J. Chem. Phys. 1977, 67, 3722. (22) Bishop, M.; Ceperley, D.; Frisch, H.; Kalos, M. Investigations of Static Properties of Two-Dimensional Bulk Polymer Systems. J. Chem. Phys. 1981, 75, 5538. (23) Baumgärtner, A. Excluded Volume Effect on Polymer Films. Polymer 1982, 23, 334. (24) Okamoto, H.; Itoh, K.; Araki, T. Scaling Relations of TwoDimensional Athermal Multichain Systems by Computer Simulation. J. Chem. Phys. 1983, 78, 975. (25) Bishop, M.; Kalos, M.; Sokal, A.; Frisch, H. Scaling in Multichain Polymer Systems in Two and Three Dimensions. J. Chem. Phys. 1983, 79, 3496. (26) Khalatur, P.; Pletneva, S.; Papulov, Y. Monte Carlo Simulation of Multiple Chain Systems: The Concentration Dependence of Osmotic Pressure in Two- and Three-Dimensional Spaces. Chem. Phys. 1984, 83, 97. 12066

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067

The Journal of Physical Chemistry B

Article

(52) Wang, Q. On the Anisotropy of Lattice Polymers. J. Chem. Phys. 2009, 131, 234903. (53) For STAR, BFM2 and BFM3 lattices, if different bonding potential energies are used as in ref 52, they need to be taken into account in the acceptance criterion of MC trial moves. (54) Jacucci, G.; Rahman, A. Comparing the Efficiency of Metropolis Monte Carlo and Molecular-Dynamics Methods for Configuration Space Sampling. Il Nuovo Cimento 1984, 4D, 341. (55) Wang, Q.; Yin, Y. Fast Off-Lattice Monte Carlo Simulations with “Soft” Repulsive Potentials. J. Chem. Phys. 2009, 130, 104903. (56) Wang, Z.-G. Chain Dimensions in Amorphous Polymer Melts. Macromolecules 1995, 28, 570. (57) See Eq (2.83) in Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Clarendon Press: Oxford, 1986. (58) Wen, P.; Zheng, N.; Li, L.; Li, H.; Sun, G.; Shi, Q. Polymerlike Statistical Characterization of Two-Dimensional Granular Chains. Phys. Rev. E 2012, 85, 031301.

12067

dx.doi.org/10.1021/jp507391j | J. Phys. Chem. B 2014, 118, 12059−12067