Correlations in Hard- and Soft-Core Generic Polymer Models

Generic polymer models capturing the chain connectivity and the non-bonded excluded-volume interactions between polymer segments can be classified into hard- and soft-core models depending on their non-bonded pair potential. Here we compared the correlation effects on the structural and thermodynamic properties of the hard- and soft-core models given by the polymer reference interaction site model (PRISM) theory, and found different behaviors of the soft-core models at large invariant degree of polymerization (IDP) depending on how IDP is varied. We also proposed an efficient numerical approach, which enables us to accurately solve the PRISM theory for chain lengths as large as 106.

While they do not correspond to any chemically specific polymer, generic polymer models are widely used in theoretical and simulation studies in the field of polymer physics as they capture two essential features of all polymers: chain connectivity and non-bonded excluded-volume interactions. Compared to atomistic models that can represent specific polymers used in experiments, molecular simulations of generic models can reach much larger length scales and much longer time scales, and theoretical studies of generic models can also be performed. Depending on whether or not the excluded-volume interactions in generic models prevent complete overlapping of polymer segments, they can be classified into hard-core models (such as those based on the hard-sphere chain model, the Kremer-Grest model [1], and the various self-and mutual-avoiding walk models on a lattice) and soft-core models (such as those used in the dissipative particle dynamics (DPD) simulation [2], fast Monte Carlo simulations [3][4][5][6][7], field-theoretic simulation (FTS) [8], variants of FTS under the partial saddle-point approximation [9], single-chain-in-mean-field simulation [10] and hybrid particle field molecular dynamics simulation [11] both under the quasi-instantaneous field approximation [10]). Taking the study of polymer melts as an example, while hard-core models have been used in conventional molecular simulations for a long time, they have the disadvantage that their chain lengths N used in such simulations (as limited by the computational cost) are too short compared to those in typical experiments; in other words, such conventional simulations significantly exaggerate the fluctuations in polymer melts compared to experiments [6,7,12]. In contrast, simulations of the more recently proposed soft-core models can readily reach the extent of fluctuations in typical experiments by increasing the chain number density (or equivalently the segment number density ρ at finite N) instead of N [6,7,12].
In this Letter we focus on a simple but important class of generic models for compressible homopolymer melts (or equivalently homopolymer solutions in an implicit solvent) in the continuum, with the excluded-volume interaction between polymer segments described by a short-range, isotropic and purely repulsive pair potential βu nb (r), where β ≡ 1/k B T with k B being the Boltzmann constant and T the thermodynamic temperature of the system; this is the basis of more complicated polymer models having attractions and/or more species. The hard-and soft-core models can then be classified according to whether drβu nb (|r|) diverges or not. This classification becomes clear after we write the total dimensionless non-bonded interaction energy for a system of n chains each having N segments in volume V under the commonly used pairwise additivity as where r i denotes the spatial position of the ith segment in the system, φ(r) ≡ ∑ nN i=1 δ(r − r i )/ρ is the segment volume fraction at r, and the last term deducting the self-interaction of segments gives an unimportant constant; while molecular simulations of this system can be performed at finite ρ ≡ nN/V for any βu nb (r) (along with a chain-connectivity model), for a homogeneous system (i.e., φ(r) = 1) the widely used polymer self-consistent field (SCF) theory [13] gives the dimensionless internal energy per chain due to the non-bonded interaction βU nb /n = (Nρ/2) drβu nb (|r|) − (N/2)βu nb (0) (due to its mean-field approximation that neglects the system fluctuations and correlations), which diverges if drβu nb (|r|) does (i.e., for the hard-core models). It is then clear that the SCF theory can only be applied to soft-core models, where one can define the dimensionless excludedvolume interaction parameter ε > 0 via u nb (r) = εu 0 (r) with the normalized pair potential u 0 (r) satisfying drβu 0 (|r|) = 1. Another necessary condition for applying the SCF theory (i.e., having finite βU nb /n) is ε ∝ ρ −1 .
Here we compare the correlation effects on the structural and thermodynamic properties of hard-and soft-core generic polymer models, which has rarely been reported [14], to further reveal their differences. For this purpose we choose the polymer reference interaction site model (PRISM) theory proposed by Schweizer and Curro [15], which has been applied to many polymeric systems, including homopolymer melts, solutions, blends, block copolymers, nanocomposites, polyelectrolytes, etc. [16][17][18][19] It can be considered as the most successful molecular-level theory to date for studying the correlations in homogeneous polymeric systems.
For the above homopolymer melts, the PRISM equation is given bŷ where h(r) is the interchain total segment pair correlation function (PCF) with r ≡ r/σ and σ the segment diameter (i.e., the range of βu nb (r)), ω(r) is the normalized (i.e., 4π ∞ 0 drr 2 ω(r) = 1) intrachain segment PCF, c(r) is the interchain direct segment PCF, f ≡ (4π/q) ∞ 0 drr f (r) sin(qr) denotes the 3D Fourier transform of a radial function f (r) with q being the wavenumber (in units of 1/σ), and ρ ≡ nNσ 3 /V is the dimensionless segment number density. For given N, ρ and ω, to solve for both h and c, a closure providing an approximate relation between them is needed; here we take the commonly used Percus-Yevick (PY) closure [20] which works well for our class of generic models where βu nb (r) is short-range and purely repulsive.
To be more specific, we consider two commonly used generic polymer models: the tangent hard-sphere chain (THSC) and the DPD models; the former is a hard-core model that uses exp −βu b (r) = δ(r − 1)/4π with βu b (r) specifying the dimensionless bonded potential between two adjacent segments on the same chain and the hard-sphere (HS) potential βu HS (r) → ∞ for r < 1 and 0 otherwise as βu nb (r), and the latter is a soft-core model that uses βu b (r) = 2r 2 and the DPD potential βu DPD (r) = (a/2)(1 − r) 2 for r < 1 and 0 otherwise as βu nb (r) with the dimensionless interaction parameter a = 15ε/π = 75/ρ chosen to mimic the compressibility of water [2]. In the thermodynamic limit, the structural and thermodynamic properties of these two models are controlled only by N and ρ; typically, molecular simulations of the DPD model uses ρ = 3 or 5.
Finally, we note thatω is needed as an input for PRISM calculations. While in general the chain conformations characterized byω depend on ρ, for simplicity here we use the ideal-chain for the THSC model and exp −q 2 /8 for the DPD model, where r ≡ r/σ, q is the wave vector for the 3D Fourier transform, and q=|q|.
We calculate h(r i ) = γ(r i ) + c(r i ) (i = 0, . . . , m-1 for the DPD model and i = 0, . . . , m for the THSC model), then use the residual errors of Equation (2) (which becomes h(r i ) = −1 for the THSC model) to converge the independent variables via the Anderson mixing [28], which has the computational complexity of O(m) and can quickly converge a large set of nonlinear equations to a high accuracy.
We use the convergence criterion of ε c < 10 −10 with ε c denoting the maximum absolute value of the residual errors of the PY closure over all r i (i = 0, . . . , m-1 for the DPD model and i = 0, . . . , m for the THSC model), and choose the values of m (=4096 for the THSC model and 512 for the DPD model) and r c (≈ 10 √ N if N < 100 and 2 √ N otherwise, rounded to the nearest integer, to capture the correlation-hole effect [29]) such that the discretization errors are comparable to ε c . Our numerical approach has the least number of independent variables to be iteratively solved, greatly reduces m (thus M) both by analytically treating the discontinuities in the THSC model and by taking the inverse Fourier transform only forγ (which decays toward 0 with increasing q much faster than bothĉ andĥ), and is essential for us to accurately solve the PRISM-PY theory for N as large as 10 6 (where for the THSC model M is about 8.2 × 10 6 !). To the best of our knowledge, analytically treating the discontinuities caused by the HS potential has not been reported in numerical calculations of even the widely studied Ornstein-Zernike (OZ) equation [30] (to which Equation (1) reduces for N = 1); in Supplemental Material we show that our numerical approach gives several orders of magnitude more accurate results than pyPRISM [19], a recently developed Python-based open-source framework for PRISM calculations.
In the limit of N→∞ and σ→0 at finite root-mean-square end-to-end distance of the ideal chain R e,0 ≡ √ N − 1σ, the THSC model becomes the hard-core Gaussian thread model [31] (HC CGC-δ, where R e,0 is taken as the unit of length); to compare the PRISM-PY results of these two models, we define two dimensionless parameters: C 0 ≡ N 2ĉ 0 σ 3 /R e,0 3 and the invariant degree of polymerization [32] N ≡ nR e,0 3 /V 2 ; N controls the fluctuations in polymer melts, and for the THSC model it is easy to show that N ∝ N at large N. Figure 1a shows how C 0 varies with N for the THSC and HC CGC-δ models; for the latter model, N is the only parameter, the PRISM-PY equation is given by Equation (18) in our previous work [14] and the corresponding numerical results for N ≥ 100 are shown in figure 8b there. We see that, while −C 0 increases monotonically with increasing N for the HC CGC-δ model, it exhibits a minimum for the THSC model. At given N due to its N→∞ the HC CGC-δ model corresponds to the limit of ρ = √ N N/(N − 1) 3/2 → 0 of the THSC model as implied in Figure 1a. At large N , we see that −C 0 ∝ √ N in all cases. This is in accordance with an asymptotic value ofĉ 0 < 0 at given ρ for the THSC model, whilê c 0 → 0 for the HC CGC-δ model. Figure 1a also shows that the DPD model at ρ = 3 gives qualitatively the same behavior of C 0 vs. N as that for the THSC model.  With the normalized isothermal compressibility given by the compressibility equation, where ρc ≡ n/V is the chain number density and   , is the isothermal compressibility with P denoting the system pressure, Figure 1b presents essentially the same data as in Figure 1a, but in a way that With the normalized isothermal compressibility κ T ≡ ρ c κ T /β = 1/ 1 − √ N C 0 given by the compressibility equation, where ρ c ≡ n/V is the chain number density and κ T ≡ −(∂V/∂P) n,β /V is the isothermal compressibility with P denoting the system pressure, Figure 1b presents essentially the same data as in Figure 1a, but in a way that can be compared with real polymers used in experiments. As shown in figure 2 of our previous work [14], κ T N ≈ 1.38 for polyethylene (at 180 • C) and 0.119 for polystyrene (at 280 • C), independent of their N ≥10 3 ; this is consistent with −C 0 ∝ √ N at large N shown in Figure 1a. On the other hand, while κ T N ∝ N is expected for very small N , the smallest N (given by N = 2) is about 0.0025, 0.090 and 0.95, respectively, for the THSC model at ρ = 0.1 and 0.6 and the DPD model at ρ = 3. Clearly, both hard-and soft-core models can be used to describe the excluded-volume interactions in real polymers, and experimental values of κ T can be achieved by adjusting ρ, for example, in the THSC and DPD models. We attribute the largest κ T at the same N given by the HC CGC-δ model to its σ→0, and note that the DPD model at ρ = 3 is actually "harder" (i.e., more difficult to compress) than the hard-core models studied here.   Figure 1c shows that at large N , the dimensionless excess (virial) pressure due to the interchain repulsion βR e,0 3 P ex = (2π/3) N N 2 /(R e,0 /σ) 3 (h(r = 1) + 1) scales with N 3/2 for the THSC model; this is due to the same scaling of R e,0 3 with N and also found for the HC CGC-δ model (where βR e,0 3 P ex = −C 0 N /2). We also see that the HC CGC-δ model gives a much smaller βR e,0 3 P ex than the THSC model at the same N , again due to its σ→0. Figure 1c further shows that at large N , the DPD model at ρ = 3 gives the same scaling of βR e,0 3 P ex = −(2π/3) N N 2 /(R e,0 /σ) 3 1 0 drr 3 (h(r) + 1)(dβu DPD (r)/dr) with N as the hard-core models; at the same N , it has even the largest βR e,0 3 P ex , in accordance with its smallest κ T shown in Figure 1b.
Note that for both the THSC and DPD models, N is varied by changing N at fixed ρ in Figure 1, which makes N and N to be approximately on the same order; it is therefore very difficult, if possible at all, to reach via this way in molecular simulations even a relatively small N -value of 10 4 used in experiments. As aforementioned, molecular simulations of soft-core models can readily reach N -values used in experiments by increasing ρ at fixed N [6,7,12]. For the DPD model at large ρ, βu DPD (r) = (75/2ρ)(1 − r) 2 ≈ 0 and the PY closure approaches the random-phase approximation (RPA) closure [33,34] c RPA (r) = −βu DPD (r), which gives c RPA 0 = −75/2ρ andĉ RPA 0 = −5π/ρ independent of N.
We then obtain κ RPA T = 1/(1 + 5πN) from the compressibility equation. Figure 2a shows κ T vs. 1/ρ obtained via the compressibility equation from our PRISM-PY calculations of the DPD model at various N, where each curve exhibits a minimum with its location shifting to smaller 1/ρ with increasing N and the intercept of each (extrapolated) curve with the left axis (i.e., in the limit of ρ → ∞ ) gives the corresponding κ RPA T . Clearly, the difference between κ T and κ RPA T is entirely due to that between the PY and RPA closures. A C 0 vs. N plot (not shown) can be obtained from Figure 2a for the DPD model. In particular, the RPA closure gives −C RPA at large N ; this is in clear contrast to −C 0 ∝ N 1/2 for the hard-core models and the DPD model at ρ = 3 shown in Figure 1a, but consistent with the soft-core Gaussian thread (SC CGC-δ) model (which is equivalent to the well-known Edwards model [35]) shown in figure 8a of our previous work [14], where N→∞ and σ→0 at finite R e,0 and βu nb (r) = κ/N 2 ρ c δ(r) is used with a finite dimensionless parameter κ > 0 controlling the repulsion strength between polymer segments. The behavior of soft-core models at large N , therefore, depends on how N is varied, i.e., whether by changing N at fixed ρ (thus the excluded-volume interaction parameter ε is fixed) or by changing ρ at fixed N (thus ε is also varied as ∝ ρ −1 ); in the former case correlations exist even in the limit of N→∞, while in the latter case the SCF theory becomes exact in the limit of ρ → ∞ (at finite N) where neither fluctuations nor correlations exist. As aforementioned, with increasing ρ at fixed N, the PY closure approaches the RPA closure, which givesĉ RPA = −(5π/ρ)βû 0 thusĥ RPA = −5πN 2 ω DPD id 2 βû 0 /ρ 1 + 5πNω DPD id βû 0 according to Equation (1). In the limit of ρ → ∞ , we have c RPA (r) → 0 and h RPA (r) → 0 , thus the SCF results of βσ 3 P SCF ex /ρ = 5π/2 and βu SCF c,ex /N = 5π/2 independent of N, where βu c,ex = 75πN 1 0 drr 2 (h(r) + 1)(1 − r) 2 denotes the dimensionless excess internal energy per chain due to the interchain repulsion. On the other hand, the differences between the SCF and RPA results as given by βσ 3 Figure 2b shows that βσ 3 P ex − P RPA ex /aρ ∝ ρ −1 at large ρ; note that P ex > P RPA ex at large ρ while P ex < P RPA ex at small ρ, which leads to the cusp of each curve shown in the figure with its location (i.e., the ρ-value at which P ex = P RPA ex ) increasing with increasing N (the cusp at N = 1 is located around ρ = 2.6). We also note that βσ 3 P SCF ex − P RPA ex /aρ ≈ 0.0327, 0.119 and 0.182 for N = 1, 10 and 100. Therefore, with increasing ρ, both βσ 3 P ex /ρ and βσ 3 P RPA ex /ρ approach βσ 3 P SCF ex /ρ. Similar results are found for β u c,ex − u RPA c,ex /aN as shown in Figure 2c, where u c,ex < u RPA c,ex at large ρ while u c,ex > u RPA c,ex at small ρ (with the cusp at N = 1 located around ρ = 2.2); also note that β u SCF c,ex − u RPA c,ex /aN ≈ 0.144, 0.253 and 0.322 for N=1, 10 and 100. In particular, the PRISM-RPA theory withω DPD id is equivalent to the Gaussian-fluctuation theory neglecting non-Gaussian fluctuations in the system and gives a correction ∝ ρ −1 to the SCF result, while the PRISM-PY theory captures non-Gaussian fluctuations in an approximate way and gives a leading-order correction ∝ ρ −2 to the Gaussian-fluctuation result. These are consistent with our previous study of compressible [36] and incompressible [37] homopolymer melts using fast lattice Monte Carlo simulations [6,7]. Given this and the agreement of our Figure 1b with experimental results at large N , we do not expect that the use of more accurateω can qualitatively change our PRISM-PY results here.
To summarize, we have compared the correlation effects on the structural and thermodynamic properties of hard-core models (i.e., the THSC model and its limit of N→∞ at finite R e,0 (or equivalently ρ → 0 at given N ), the HC CGC-δ model [31]) and soft-core models (i.e., the DPD model and its limit of N→∞ at finite R e,0 , the Edwards model [35]) for compressible homopolymer melts (or equivalently homopolymer solutions in an implicit solvent) given by the PRISM-PY theory. The behavior of soft-core models at large N depends on how N is varied, i.e., whether by changing N at fixed ρ (thus ε is fixed) or by changing ρ at fixed N (thus ε is also varied as being inversely proportional to ρ). In the former case, correlations exist even in the limit of N→∞, and both the hard-core and the DPD models give −C 0 ∝ N 1/2 at large N , consistent with real polymers used in experiments; it is, however, very difficult to reach via this way in molecular simulations even a relatively small N -value of 10 4 used in experiments. This problem is solved in the latter case, where the widely used polymer SCF theory becomes exact in the limit of ρ → ∞ (at finite N), the Gaussian-fluctuation theory gives a correction ∝ ρ −1 to the SCF result, and the PRISM-PY theory captures non-Gaussian fluctuations in the system in an approximate way and gives a leading-order correction ∝ ρ −2 to the Gaussian-fluctuation result, consistent with our previous simulations [36,37]. The soft-core models, however, give −C 0 ∝ N −1/2 at large N , suggesting that it would be difficult, if possible at all, for the various recently proposed simulation methods [3][4][5][6][7][8][9][10][11] to capture both the fluctuations and correlations in experimental systems. We also proposed an efficient numerical approach, which enables us to accurately solve the PRISM-PY theory for N as large as 10 6 ; numerical calculations of such theories can, therefore, capture both the fluctuations and correlations in experimental systems.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/polym15051180/s1, Figure S1: Logarithmic plot of the numerical errors given by pyPRISM and our approach for the HS model; Table S1: List of variables used in the main text.
Funding: This research and the APC were funded by the American Chemical Society Petroleum Research Fund (grant number 66122-ND6). Acknowledgment is made to the donors of The American Chemical Society Petroleum Research Fund for partial support of this research.

Data Availability Statement:
The in-house code written in C and the data presented in this study are available upon reasonable request from the corresponding author.

Conflicts of Interest:
The author declares no conflict of interest.