Ultra-coarse-graining of homopolymers in inhomogeneous systems

The microscopically resolved (MR) reference systems consisted of N_p homopolymers dispersed in an implicit solvent. Polymers were represented by a generic bead-spring model with N = 20 monomers per chain. Non-bonded interactions between the monomers were modeled using the Lennard-Jones (LJ) potential:

$$\begin{equation}{U}^{\text{LJ}}\left({r}_{ij}\right)=\begin{cases}4\varepsilon \left[{\left(\frac{\sigma }{{r}_{ij}}\right)}^{12}-{\left(\frac{\sigma }{{r}_{ij}}\right)}^{6}\right],\quad \hfill & \quad {r}_{ij}{< }{r}_{\text{c}}\hfill \\ 0,\quad \hfill & \quad {r}_{ij}{\geqslant}{r}_{\text{c}}\hfill \end{cases}\end{equation} \tag{ 1 }$$

with radial distance r_ij between particles i and j, interaction strength ɛ = k_B T, and bead diameter σ. The potential was truncated at the cutoff radius of r_c = 5σ.

Polymer bonds were modeled through the finitely extensible nonlinear elastic (FENE) potential: [50]

$$\begin{equation}{U}^{\text{FENE}}\left({r}_{ij}\right)=\begin{cases}-\frac{1}{2}k{r}_{0}^{2}\enspace \mathrm{ln}\left[1-{\left({r}_{ij}/{r}_{0}\right)}^{2}\right],\quad \hfill & \quad {r}_{ij}{< }{r}_{0}\hfill \\ \infty ,\quad \hfill & \quad {r}_{ij}{\geqslant}{r}_{0}\hfill \end{cases}\end{equation} \tag{ 2 }$$

with spring constant k = 30k_B T/σ² and maximum bond extension r₀ = 1.5σ [51].

All MD simulations were conducted in the $\mathcal{N}VT$ ensemble (unless stated otherwise explicitly), $\mathcal{N}={N}_{\text{p}}N$ being the total number of monomers in the system. The temperature was kept constant at T = 1.0ɛ/k_B through a Langevin thermostat with friction coefficient ξ = 1.0m/τ, where m is the monomer mass, and $\tau =\sqrt{m{\sigma }^{2}/\left({k}_{\text{B}}T\right)}$ is the unit of time. The equations of motion were integrated using a Verlet scheme with time step Δt = 0.001τ. All simulations were conducted in cubic simulation boxes with edge length L and periodic boundary conditions in all three Cartesian directions.

To characterize the polymer properties in homogeneous systems, we simulated a bulk polymer melt with N_p = 2000 chains. The melt was first equilibrated in the $\mathcal{N}PT$ ensemble at a pressure of P = 0k_B T/σ³, resulting in an average monomer number density of ρ_b = 0.921σ⁻³ (all quantities extracted from bulk simulations will be denoted by a subscript 'b', while the subscript '0' marks quantities obtained or parameterized at infinite dilution). Once the system reached equilibrium, we switched to the $\mathcal{N}VT$ ensemble using a cubic simulation box with L = 35.15σ, and simulated for 2 × 10⁷ time steps, saving configurations and taking measurements every 1000 time steps. We characterized the conformation of the polymers through the average radius of gyration tensor:

$$\begin{equation}\mathbf{G}=\frac{1}{N}\sum\limits _{i=1}^{N}\left\langle {\Delta}{\mathbf{r}}_{i}\otimes {\Delta}{\mathbf{r}}_{i}\right\rangle \end{equation} \tag{ 3 }$$

where Δr_i is the vector from the polymer center of mass to monomer i, and ⊗ indicates the dyadic product. The radius of gyration is then ${R}_{\text{g}}=\sqrt{{G}_{xx}+{G}_{yy}+{G}_{zz}}$, with R_g,b ≈ 2.2σ in the bulk simulations.

Thin films were constructed by initially placing the polymers in a slab at the center of the simulation box, so that the normal vectors of the polymer–solvent interfaces lied parallel to the z-axis. The polymer droplets were created similarly by initializing the polymers close to each other in the box center. The size of the simulation boxes was chosen sufficiently large to prevent nonphysical self-interactions and the coalescence of the films/droplets with their periodic images (see table 1). Simulations were then performed in the $\mathcal{N}VT$ ensemble at T = 1.0ɛ/k_B. This state point falls inside the two-phase coexistence region of the system [52], so that the polymers separated into a high and low density phase that coexisted in the same simulation box. Simulation snapshots of typical films and droplets from the MR simulations are shown in figure 1. We determined the film thicknesses H and droplet diameters 2R from the full width at half maximum of the monomer density profiles (see figure 2(a)). The resulting H and R values are summarized in table 1, and they are in good agreement with the estimates $H=\mathcal{N}/\left({L}^{2}{\rho }_{\text{b}}\right)$ and $R={\left[3\mathcal{N}/\left(4\pi {\rho }_{\text{b}}\right)\right]}^{1/3}$, respectively, which assume a homogeneous monomer density inside the films/droplets and a perfectly sharp polymer–solvent interface.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2(a) shows the monomers number density along the film normal, ρ(z), for the system containing N_p = 1245 homopolymers (the data for the thinner films look qualitatively similar and have been omitted for brevity). Since the simulated temperature was far below the critical temperature of this system (T_c ≈ 2.65ɛ/k_B [52]), all polymers were part of the thin film (the density in the 'vapor' phase was strictly zero in our simulations). The density profile has a flat plateau with density ρ = 0.922σ⁻³, which is in excellent agreement with the value measured in the bulk systems (ρ_b = 0.921σ⁻³). Measurement uncertainties of the density profiles were estimated from the standard error of the mean between ten subdivided blocks of the data. Note that we did not determine the (apparent) width of the polymer–solvent interfaces, as this measurement would only be meaningful when the dependence on the lateral box dimensions is analyzed to account for the broadening due to capillary waves [54].

The density profile of the polymer centers of mass, ρ_p(z), is also flat in the bulk-like interior of the film, but has two distinct peaks at z = ±8σ, indicating a surplus of polymers near the polymer–solvent interfaces (figure 2(a)). This local excess of polymers can be rationalized by considering the conformation of the individual chains, which we have characterized via the diagonal components of the radius of gyration tensor (figure 2(b)). The polymers are essentially isotropic in the central region of the film, but they assume an oblate ellipsoidal shape near the edges of the film, with G_zz < G_yy = G_xx and a maximum aspect ratio of about G_xx /G_zz ≈ 4.

Finally, we determined the surface tension of the thin films from the anisotropy of the pressure tensor:

$$\begin{equation}\gamma =\frac{{L}_{z}}{2}\left\langle {P}_{zz}-\frac{{P}_{xx}+{P}_{yy}}{2}\right\rangle \end{equation} \tag{ 4 }$$

where P_αα denotes the diagonal components of the instantaneous pressure tensor, which was computed from the Clausius virial equation. The factor of 1/2 in equation (4) is due to the presence of two interfaces (see figure 1(a)). Here, we find a value of γ = 1.42k_B T/σ² for all three investigated film thicknesses.

In our CG model, an entire homopolymer chain is represented by a single spherical particle, as depicted in figure 3. We tested three different approaches, using (i) only two-body interactions, (ii) a combination of two-body and three-body interactions, and (iii) a combination of two-body and (mean-field) many-body interactions. In the following, we briefly discuss these interactions and how they have been parameterized. The resulting CG potentials are then presented in section 3.

Zoom In Zoom Out Reset image size

2.2.1. Two-body interaction

In all models, the pairwise interactions were determined using FM [18, 24, 25], where the force on each CG bead is calculated as the sum of the forces on the monomers of the corresponding polymer. We determined the two-body interaction potentials between the CG polymers using FM from bulk simulations, ${U}_{\text{b}}^{2\text{b}}$, as well as from simulations containing two isolated polymers, ${U}_{0}^{2\text{b}}$. The first approach results in a pair potential that (implicitly) includes many-body effects due to the surrounding polymers, while the second approach provides the pairwise interaction of only two polymer chains.

We introduced a cutoff radius, ${r}_{\text{c}}^{2\text{b}}$, beyond which the CG particles do not interact with each other. The two-body forces in the CG model, f^2b, were parameterized using cubic splines with a uniform grid spacing Δr^2b. Thus, there are $K\equiv {r}_{\text{c}}^{2\text{b}}/{\Delta}{r}^{2\text{b}}$ grid points and basis functions which depend linearly on 2K spline coefficients λ_i . Imposing continuity and smoothness conditions, reduces the number of free coefficients to K. These coefficients were then determined by matching the forces in the CG representation to the ones from the MR simulations. To determine the K coefficients, we first divided the MR trajectory in blocks containing N_s snapshots each, and then solved the resulting set of N_s N_p linear equations with a constrained least-squares algorithm for each block. Finally, we tabulated f^2b with a grid spacing of ${\Delta}{r}_{\text{tab}}^{2\text{b}}=0.01\sigma$, averaging over the blockwise results. To ensure a smooth decay to zero at ${r}_{\text{c}}^{2\text{b}}$, we multiplied the tabulated forces with a smoothing function f^sm for all distances r > r^sm:

$$\begin{equation}{f}^{\text{sm}}\left(r\right)=\mathrm{cos}\left(\frac{\pi }{2}\frac{r-{r}^{\text{sm}}}{{r}_{\text{c}}^{2\text{b}}-{r}^{\text{sm}}}\right).\end{equation} \tag{ 5 }$$

The CG pair potential U^2b was then obtained by numerical integration of f^2b. This coarse-graining procedure is implemented in the VOTCA-CSG package [55].

The pair potential ${U}_{\text{b}}^{2\text{b}}$ was determined from the MR melt simulations (see section 2.1) using a cutoff distance of ${r}_{\text{c}}^{\mathrm{2}\mathrm{b}}=8\enspace \sigma$ and a grid spacing of Δr^2b = 0.1σ. Each block contained 300 frames, and the final potential was computed by taking the average of the 67 blocks. For computing ${U}_{0}^{2\text{b}}$, we performed 64 independent simulations of two isolated polymers in a large simulation box (L = 60σ), so that they did not interact with their periodic images. The polymers were initialized at large distances and then approached each other due to the monomer–monomer attraction (see section 2.1). By evaluating the trajectories of these multiple runs, there were enough data for capturing the forces at all relevant distances. For the FM procedure for calculating ${U}_{0}^{2\text{b}}$, we used ${r}_{\text{c}}^{\mathrm{2}\mathrm{b}}=8\enspace \sigma$ and Δr^2b = 0.2σ, averaging over 14 blocks each containing 750 000 frames. In both cases, we chose r^sm = 7.5σ.

2.2.2. Three-body interaction

Three-body interactions were taken into account using the Stillinger–Weber (SW) potential with flexible angular dependence, as implemented in the VOTCA-CSG package [29] and available under https://gitlab.mpcdf.mpg.de/votca

$$\begin{align}\hfill {U}^{\text{SW}}\left({r}_{IJ},{r}_{IK},{\theta }_{IJK}\right)& ={f}^{\text{SW}}\left({\theta }_{IJK}\right)\mathrm{exp}\left(\frac{\eta }{{r}_{IJ}-{r}_{\text{c}}^{\text{SW}}}\right)\enspace \hfill \\ \hfill & \quad {\times}\mathrm{exp}\left(\frac{\eta }{{r}_{IK}-{r}_{\text{c}}^{\text{SW}}}\right),\hfill \end{align} \tag{ 6 }$$

with angle θ_IJK between the three particles (I being the central one), and angular interaction term f^SW(θ_IJK ). The cutoff distance ${r}_{\text{c}}^{\text{SW}}$ determines how many triplets are included into the local environment of each CG polymer, and U^SW smoothly decays to zero when one of the two inter-particle distances r_IJ and r_IK reaches ${r}_{\text{c}}^{\text{SW}}$. The parameter η controls the steepness of this decay with small η corresponding to a steeper transition.

The angular dependence of U^SW allows for capturing anisotropic interactions, which occur near the polymer–solvent interfaces in the MR simulations (see figure 2). The angular interaction term ${f}^{\text{SW}}\left({\theta }_{IJK}\right)$ was fitted to the residual forces of the MR simulations Δf via the FM procedure [29]. We determined Δf acting on each CG polymer chain by subtracting the CG two-body force from the total reference force from the MR simulations, f_MR:

$$\begin{equation}{\Delta}\mathbf{f}={\mathbf{f}}_{\text{MR}}-{\mathbf{f}}^{2\text{b}}.\end{equation} \tag{ 7 }$$

In practice, we determined Δf by recomputing the forces acting on the centers of mass of the polymers in the MR trajectories using the pairwise forces f^2b. For f^2b, we either used the forces due to ${U}_{\text{b}}^{2\text{b}}$ or ${U}_{0}^{2\text{b}}$ (see section 2.2.1). To distinguish the resulting three-body potentials, we will refer to them as ${U}_{{\Delta}\mathrm{b}}^{\text{SW}}$ and ${U}_{{\Delta}0}^{\text{SW}}$, respectively. The subscript 'Δ' indicates that the SW potential was parameterized using the residual forces Δf (equation (7)). Parameterization of U^SW according to Δf ensures the orthogonality of the two-body and three-body terms [29]. The fit parameters then depend on the choice of the pair potential, so that the additional three-body contribution can be regarded as a correction term to either ${U}_{\text{b}}^{2\text{b}}$ or ${U}_{0}^{2\text{b}}$.

All trajectory reruns with pair potentials were carried out with the GROMACS simulation package [56]. All reruns with three-body potentials were carried out with the LAMMPS simulation package [57] with the user pair style sw/table [29] as available under https://gitlab.mpcdf.mpg.de/votca/lammps.

We used a cubic spline representation for ${f}^{\text{SW}}\left({\theta }_{ijk}\right)$ with K^SW = 31 grid points and a linear dependence on the 2K^SW spline coefficients. Treating the two exponential terms in equation (6) as prefactors led to a linear set of equations which was solved by a constrained least-squares solver. For the thinnest films with H ≈ 4R_g,b, we averaged over 200 blocks containing each 200 frames, while for the remaining systems we averaged over 800 blocks each containing 50 frames. We systematically varied the cutoff radius of ${U}_{{\Delta}\mathrm{b}}^{\text{SW}}$ and ${U}_{{\Delta}0}^{\text{SW}}$ from ${r}_{\text{c}}^{\text{SW}}=4\enspace \sigma$ to ${r}_{\text{c}}^{\text{SW}}=10\enspace \sigma$, and the parameter η from η = σ to η = 4σ.

2.2.3. Local density-dependent interaction

Alternatively, many-body interactions were included in our CG simulations by supplementing the two-body pair potential with a local density-dependent potential (LDP). In this representation, each CG polymer carries a cloud, which effectively describes its monomers that were integrated out during the coarse-graining procedure (see figure 3). For simplicity, we assumed that the overlap of these clouds can be described by radial weight functions ω(r_IJ ), which only depend on the distance r_IJ between the CG particles I and J. In this work, we constructed the weight function ω(r_IJ ) by computing the pairwise overlap integral between the average monomer density clouds around the polymers centers of mass, ρ_cloud(r). To determine ω(r_IJ ), we first computed ρ_cloud(r) from the MR bulk simulations:

$$\begin{equation}{\rho }_{\text{cloud}}\left(r\right)=\left\langle \frac{1}{N4\pi {r}^{2}}\sum\limits _{i=1}^{N}\delta \left(\left\vert {\mathbf{r}}_{i}-{\mathbf{r}}_{\text{com}}\right\vert -r\right)\right\rangle ,\end{equation} \tag{ 8 }$$

where r_com is the center of mass position of the polymer. Assuming a pairwise overlap of the density clouds, the weight function ω(r_IJ ) is then given by:

$$\begin{equation}\omega \left({r}_{IJ}\right)={\int }_{V}{\rho }_{\text{cloud}}\left(\left\vert {\mathbf{r}}_{\text{com}}^{I}-\mathbf{r}\right\vert \right){\rho }_{\text{cloud}}\left(\left\vert {\mathbf{r}}_{\text{com}}^{J}-\mathbf{r}\right\vert \right)\enspace {\mathrm{d}}^{3}\mathbf{r}.\end{equation} \tag{ 9 }$$

Thus, ω(r_IJ ) can be interpreted as the number density of monomer pairs between a pair of (CG) polymers at center-to-center distance r_IJ . Figure 4 shows ρ_cloud(r) as well as ω(r), which are both bounded and monotonically decrease with increasing r. The monomer density decays almost completely at r ≈ R_g,b while ω(r) vanishes almost entirely at r ≈ 2R_g,b, as expected.

Zoom In Zoom Out Reset image size

The potential energy of CG particle I due to the LDP is then:

$$\begin{equation}{U}_{I}^{\text{LDP}}=G\left({\varphi }_{I}\right)\end{equation} \tag{ 10 }$$

where G is an embedding function, and φ_I is the local density of (fictitious) monomer pairs at the position of CG particle I due to the other CG particles. We posit that φ_I can be expressed as the linear superposition of the pair density clouds of the other particles, that is:

$$\begin{equation}{\varphi }_{I}=\sum\limits _{J\ne I}\omega \left({r}_{IJ}\right).\end{equation} \tag{ 11 }$$

In this context, the embedding function G is a mean-field representation of many-body effects [58], and local density gradients enter the force calculation via the gradient of the embedding function. Thus, many-body effects are related to the curvature of G: if G is a linear function of φ_I , then the total potential reduces to a two-body pair potential. If G has a positive curvature, i.e., d² G/dφ² > 0, then the LDP becomes more repulsive as the local density increases. Due to the assumed pairwise additivity of the weight functions ω(r_IJ ), the LDP can be interpreted as a generalized pair potential, which is why in practice LDPs are typically implemented as pair potentials [59–61]. Hence, the pressure tensor can be calculated as usual via the standard Clausius virial equation [62].

The embedding function was implemented as a cardinal B-spline function, with equally spaced nodes λ_i at distances Δφ = 0.5. This representation has the advantage that the derivatives with respect to the amplitude of the nodes, ∂U^LDP/∂λ_i , are linear, which facilitates the calculation in the update step (see equation (15)). The number of nodes was adapted during the coarse-graining procedure to cover the required density range. We determined the embedding function G from our MR simulations using relative entropy minimization (REM) [20, 63]. The goal of this technique is to optimize the CG potential in such a way that the difference between the probability distribution of the MR configurations in the CG and the MR representation is minimized. This difference can be measured by the relative entropy, also known as the Kullback–Leibler divergence [64], which can be interpreted as the difference in information conveyed by the probability distributions. So in contrast to other coarse-graining methods like FM or IBI, which are designed to reproduce the potential of mean force or the radial distribution function, respectively, REM aims to minimize the difference in the configurational probability distributions. As a consequence, the resulting CG potentials will differ, because there is in general no unique parameterization which optimizes all quantities of the target system.

In the canonical ensemble, the relative entropy can be written as:

$$\begin{equation}{S}_{\text{rel}}=\beta {\left\langle {U}_{\text{CG}}-{U}_{\text{MR}}\right\rangle }_{\text{MR}}-\beta \left({A}_{\text{CG}}-{A}_{\text{MR}}\right)+{\left\langle {S}_{\text{map}}\right\rangle }_{\text{MR}}\end{equation} \tag{ 12 }$$

where U_MR and U_CG denote the potential functions, and A_MR and A_CG the configurational parts of the Helmholtz free energy for the MR and CG system, respectively. The term ${\left\langle {S}_{\text{map}}\right\rangle }_{\text{MR}}$ is the (unavoidable) contribution to the relative entropy due to the mapping. The relative entropy S_rel can be optimized using approximate and/or iterative methods to find a (locally) optimal set of parameters of the potential. In this work, we used a Newton–Raphson update rule to find the root of ∂S_rel/∂λ, that is

$$\begin{equation}{\lambda }^{j+1}={\lambda }^{j}-\left(\frac{\partial {S}_{\text{rel}}}{\partial \lambda }\right)/\left(\frac{{\partial }^{2}{S}_{\text{rel}}}{\partial {\lambda }^{2}}\right)\end{equation} \tag{ 13 }$$

where j denotes the iteration. Here, the parameters to be optimized are the amplitudes of the nodes λ_i of the cardinal B-spline representation of the embedding function G. To improve the stability of the updates and the convergence of the optimization, we multiplied the change of the parameter Δλ_i with a constant factor α = 0.01 and clipped it at ±1k_B T. The resulting update rule for the coefficients is then:

$$\begin{equation}{\lambda }_{i}^{j+1}={\lambda }_{i}^{j}-{\Delta}{\lambda }_{i}^{j}\end{equation} \tag{ 14 }$$

$$\begin{align}\hfill {\Delta}{\lambda }_{i}^{j}& =\mathrm{min}\left(-1,\mathrm{max}\left(1,\alpha \left[{\left\langle \frac{\partial {U}^{\text{LDP}}}{\partial {\lambda }_{i}}\right\rangle }_{\text{MR}}-{\left\langle \frac{\partial {U}^{\text{LDP}}}{\partial {\lambda }_{i}}\right\rangle }_{\text{CG}}\right]\right.\right.\hfill \\ \hfill & \left.\left.\quad {\times}\left[{\left\langle \frac{{\partial }^{2}{U}^{\text{LDP}}}{\partial {\left({\lambda }_{i}\right)}^{2}}\right\rangle }_{\text{MR}}-{\left\langle \frac{{\partial }^{2}{U}^{\text{LDP}}}{\partial {\left({\lambda }_{i}\right)}^{2}}\right\rangle }_{\text{CG}}+\beta {\left\langle {\left(\frac{\partial {U}^{\text{LDP}}}{\partial {\lambda }_{i}}\right)}^{2}\right\rangle }_{\text{CG}}\right.\right.\right.\hfill \\ \hfill & \left.\left.{\left.\quad -\beta {\left\langle \frac{\partial {U}^{\text{LDP}}}{\partial {\lambda }_{i}}\right\rangle }_{\text{CG}}^{2}\right]}^{-1}\right)\right)\hfill \end{align} \tag{ 15 }$$

where ${\left\langle \cdot \right\rangle }_{\text{MR}}$ and ${\left\langle \cdot \right\rangle }_{\text{CG}}$ indicate averaging in the mapped MR and CG representation, respectively. For all investigated systems, we used 500 iterations with 200 000 simulation timesteps per iteration. All simulations with LDPs were carried out on graphics processing units with the HOOMD-blue simulation package (v. 2.4.2) [59–61].

In a bulk polymer melt, every monomer is isotropically surrounded by other monomers, both from the same chain as from the other ones. Consequently, the intra- and inter-chain interactions compensate each other, and the polymers behave as dispersed in a Θ-solvent [2, 65, 66]. ³ This does, however, not mean that the effective polymer–polymer interactions vanish, but that they are purely entropic instead. Figure 5 shows the effective pair potential between CG polymers in a melt, ${U}_{\text{b}}^{2\text{b}}$, which is repulsive and bounded, allowing a partial overlap of the CG polymers. This behavior is due to the fractal and open nature of the polymer coils, which lets the centers of mass of two coils be at the same place while each chain can fluctuate without intersecting the other (see figure 3(a)).

Zoom In Zoom Out Reset image size

Previous consideration within renormalization-group theory [70] predicted that the effective potential can be approximated rather accurately by a Gaussian function ${U}_{\text{b}}^{2\text{b}}\propto \mathrm{exp}\left[-{\left({r}_{IJ}/{\sigma }_{\text{p}}\right)}^{2}\right]$, with characteristic length-scale of the interaction σ_p. Our simulation data can be fitted perfectly to this functional form with σ_p = 3.4σ. Such Gaussian polymer–polymer interactions have also been reported in previous lattice [71] and off-lattice [1, 70] simulations of athermal polymer solutions. In particular, Louis et al varied the polymer concentration over a wide range, from dilute solutions up to almost five times the overlap concentration, finding that ${U}_{\text{b}}^{2\text{b}}$ barely changed [1]. Further, they showed that ${U}_{\text{b}}^{2\text{b}}$ reproduces rather accurately the equation of state of a polymer solution. However, a CG description of the polymers only in terms of ${U}_{\text{b}}^{2\text{b}}$ is insufficient for modeling polymer films and droplets, as the CG polymers will repel each other rather than form a stable film, unless additional external constraints are applied.

The pair potential between two isolated polymers, ${U}_{0}^{2\text{b}}$, is attractive (see figure 5), resembling a soft square-well potential with a well depth of about −13k_B T and a well width of roughly 2σ. This strong polymer–polymer attraction is due to the attractive tail of the LJ interaction between the monomers in the MR model (see equation (1)). At such low polymer concentrations, the employed MR model effectively describes polymers in a poor (implicit) solvent [52, 72, 73], which have collapsed into compact globules with a distinctly smaller radius of gyration compared to their coil-like analogs in the melt, i.e., R_g,0 ≈ 1.5σ vs R_g,b ≈ 2.2σ. Once two (or more) of these globules approach each other in the MR simulations, they coalesce into a small polymer droplet that remains stable over the investigated simulation time. However, simulations containing many CG polymers interacting only through ${U}_{0}^{2\text{b}}$ behave in a nonphysical manner: due to the strong attraction for r ≲ 5σ and the lack of repulsive excluded volume interactions at short distances, the CG polymers merge into a single small aggregate, with a diameter that even decreases with increasing aggregation number. Alternatively, one could try to parameterize the pair potential at the Θ-temperature, where the solvent quality is exactly poor enough to cancel the polymer-swelling due to excluded volume effects, so that the polymer scales like an ideal chain. However, the resulting pair potential between two isolated polymers will be zero (apart from statistical fluctuations), leading to incorrect thermodynamic properties in the bulk and confined systems.

We determined the SW three-body potentials for the three different film systems, as well as for the droplet with radius R ≈ 6R_g,b. For all investigated systems, the forces of the reference systems, f_MR, were reproduced most faithfully when the residual forces Δf were fitted using ${U}_{0}^{2\text{b}}$ with a cutoff distance of ${r}_{\text{c}}^{\text{SW}}=10\enspace \sigma$. The dependence on η was weak and we chose η = σ (see the supporting information (https://stacks.iop.org/JPCM/33/254002/mmedia) for a comparison of all different parameter combinations, as well as the results when fitting to the residual forces using ${U}_{\text{b}}^{2\text{b}}$).

Figure 6(a) shows the average force on the CG polymers perpendicular to the film surface, F_z (z), in the film with H ≈ 8R_g,b, where measurements were taken every 100τ. (The results for the other systems are qualitatively similar and are included in the supporting information.) In the MR simulations, the average net force on the center of mass of the polymers vanishes in the central bulk-like region of the film, |z| < 6σ due to the film's symmetry along the z-axis. Further away from the film center, the net force on the polymers points inwards with its magnitude increasing the further the chain is located away from the film center. This net attraction reflects the cohesive forces between the monomers in the MR simulations.

Zoom In Zoom Out Reset image size

Using the particle positions from the MR trajectories (see section 2.2.2 for details), we computed the force profiles in the CG model and show the resulting profiles also in figure 6(a). The resulting net force has a similar shape as the reference force, although the flat plateau in the film center is replaced by a weak oscillation. The contribution from the two-body potential ${U}_{0}^{2\text{b}}$ leads to a strong net attraction between the CG polymers, which is especially pronounced near the film surfaces because of the inhomogeneous polymer density distribution (see figure 2). The pairwise attraction in the CG model extends deep into the bulk region of the film, which is (partially) compensated by the repulsive three-body SW interactions. It should be noted, that these force profiles have been created using the CG model in trajectories from the MR reference simulations, and not from dedicated simulations using the CG model. Thus, the good agreement of the force profiles shown in figure 6(a) does not guarantee the stability of the films and droplets in the CG representation, which we will test in section 3.4.

In figure 6(b), we plot the angular part of the SW potential f^SW(θ_IJK ) for the three different films as well as for the droplet with radius R ≈ 6R_g,b. The fitted functions f^SW(θ_IJK ) have similar shapes in all investigated cases, indicating a good transferability of the three-body potential. As the two exponential terms of the SW potential are strictly positive (see equation (6)), a positive (negative) sign of f^SW(θ) results in a net repulsion (attraction) between the central CG particle I of a triplet configuration and the two CG particles J and K. The SW potential is purely repulsive for small angles θ_IJK < 30°, where the three CG polymers I, J, and K are lined up with J at a smaller distance and K at a larger distance behind chain J (see figure 6(c)). For larger angles θ_IJK > 30°, f^SW(θ) oscillates around zero with a rather small amplitude.

Using ${U}_{0}^{2\text{b}}$ as the two-body potential between the CG polymers, we determined ${U}_{0}^{\text{LDP}}$ for all films and droplets via REM, as described in section 2.2.3. To verify that the coarse-graining procedure converged, we computed the probability density function of the local density, P(φ), in both the mapped MR reference simulations as well as in the CG simulations. Figure 7(a) shows the corresponding results for the film with thickness H ≈ 8R_g,b and the droplet with radius R ≈ 4R_g,b. The data for the CG and MR simulations are in perfect agreement, which indicates that the LDPs have been parameterized accurately. Small φ values correspond to polymers located near the polymer–solvent interface, while large φ values are associated with polymers in the central bulk-like region of the systems. All P(φ) have a maximum near φ ≈ 15, which indicates that there are more polymers in the bulk-like region than at the surfaces. This maximum becomes more pronounced as the film thickness H (droplet radius R) increases, because the surface-to-volume ratio decreases as 2/H (3/R). Note also, that a droplet with diameter 2R = H has a larger surface-to-volume ratio than a film with the same thickness, which is reflected by the larger P(φ) values at small φ for droplets compared to films (see figure 7(a)).

Zoom In Zoom Out Reset image size

Figure 7(b) shows the fitted embedding functions G(φ) for all CG systems, which look nearly identical except for a slight vertical shift in the region φ ≳ 4. This good agreement indicates that the obtained LDPs are transferable across the different systems, that is, the embedding function G(φ) parameterized for a small droplet is also applicable to a thick film. G(φ) increases linearly with φ for 5 ⩽ φ ⩽ 22, which corresponds to the range of local densities typically observed in the films and droplets (see figure 7(a)). For such linear embedding functions G(φ) = kφ + C, the potential energy of particle I due to the LDP can be written as:

$$\begin{equation}{U}_{I}^{\text{LDP}}=C+k{\varphi }_{I}=C+\sum\limits _{I\ne J}k\omega \left({r}_{IJ}\right)=C+\sum\limits _{I\ne J}{U}^{\text{LDP,}2\text{b}}\left({r}_{IJ}\right).\end{equation} \tag{ 16 }$$

Hence, in this case, the LDP behaves like an additive pair potential U^LDP,2b(r), which restores the excluded volume interactions between the CG polymers. Figure 7(c) shows the resulting effective pair interaction between CG particles in the bulk-like region of the films and droplets, ${U}_{0}^{2\text{b}}+{U}^{\text{LDP,}2\text{b}}$, which is repulsive at short distances to prevent the CG particles from collapsing onto each other, and has a minimum at r ≈ R_g,b = 2.2σ to achieve the desired inter-particle spacing in the bulk-like region. Outside the density region 5 ⩽ φ ⩽ 22, the derivative d² G/dφ² is nonzero (see inset of figure 7(b)), so that many-body effects become relevant. For smaller φ < 5, the derivative d² G/dφ² is negative, so that ${U}_{0}^{\text{LDP}}$ becomes less repulsive as the local density increases. In contrast, d² G/dφ² > 0 for φ > 22, which ensures the incompressibility of the polymer film.

To assess the computational performance of the CG model, we performed bulk simulations with N_p = 2000 polymers in a cubic box with edge length L = 35.15σ. Here, we found that the CG simulations using only the two-body interactions ${U}_{\text{b}}^{2\text{b}}$ achieved about 40 times more timesteps per second than the MR simulations, while the CG simulations with both ${U}_{0}^{2\text{b}}$ and ${U}_{0}^{\text{LDP}}$ achieved about 10 times more timesteps per second than the MR simulations. The actual speedup with respect to the physical timescales is likely even higher due to the accelerated dynamics in the soft CG models.

To investigate whether the CG representations using the SW potential or the LDP are able to reproduce stable films and droplets, we performed MD simulations using the CG models starting from the final snapshots of the MR reference simulations. Figure 8 shows typical snapshots of the CG simulations of the films with thickness H ≈ 8R_g,b at different simulation times. The starting configuration is shown in figure 8(a), while the results from the runs with the SW potential and LDP are shown in figures 8(b), (c) and (d), (e), respectively.

Zoom In Zoom Out Reset image size

The films in the CG simulations with ${U}_{0}^{2\text{b}}$ and ${U}_{{\Delta}0}^{\text{SW}}$ start to become unstable already after about 2τ, as the CG polymers collapse into small spherical aggregates, which are stabilized by a long range repulsion between them. The inability of this combination of two- and three-body interactions to reproduce the MR simulations even qualitatively stems likely from the restricted functional form of ${U}_{{\Delta}0}^{\text{SW}}$. The residual force fit of the SW potential (see section 2.2.2) most probably only captures a local minimum of the potential of mean force, whereas the global minimum of the fitted CG model corresponds to this collapsed state. This transition can be better understood by considering the strong attraction of ${U}_{0}^{2\text{b}}$ at short distances (see figure 5), whereas the three-body contribution ${U}_{{\Delta}0}^{\text{SW}}$ is repulsive only for small angles θ_IJK < 30°, where three CG beads are lined up as shown in figure 6(b). Thus, small deviations in θ_IJK are already sufficient to overcome the repulsive barrier between neighboring particles. This interpretation is corroborated by the distribution of angles θ_IJK in the collapsed droplets, where we find almost no triplets with angles in the range of θ_IJK < 30°. CG simulations with ${U}_{\text{b}}^{2\text{b}}$ and ${U}_{{\Delta}\mathrm{b}}^{\text{SW}}$ do not reproduce stable films or droplets either, as in these cases the (attractive) three-body contribution ${U}_{{\Delta}\mathrm{b}}^{\text{SW}}$ is not sufficient to compensate the repulsive pair potential ${U}_{\text{b}}^{2\text{b}}$ (see the supporting information for the fitting results of ${U}_{{\Delta}\mathrm{b}}^{\text{SW}}$). In addition, including SW interactions with a relatively large cutoff (${r}_{\text{c}}^{\text{SW}}=10\enspace \sigma$) significantly slows down the CG simulations due to the large number of triplets included into the force calculation. In fact, the CG simulations became less efficient than the original MR ones.

In contrast, the CG simulations with the pair potential ${U}_{0}^{2\text{b}}$ and LDP ${U}_{0}^{\text{LDP}}$ lead to stable film and droplet configurations, even after long simulation times (see figure 8). (Note that simulations using the effective pair interactions ${U}_{0}^{2\text{b}}$ and U^LDP,2b, see section 3.3, exhibit a similar instability as the simulations with ${U}_{0}^{2\text{b}}$ and ${U}_{{\Delta}0}^{\text{SW}}$.) To investigate the resulting configurations in more detail, we computed the density profiles along the z-axis for the films and in the radial direction for the droplets. The corresponding results for the film with thickness H ≈ 8R_g,b and the droplet with radius R ≈ 4R_g,b are plotted in figure 9. Again the measurement uncertainties of the density profiles are estimated from the standard error of the mean between ten subdivided blocks of the data. Overall, the density profiles in the CG simulations are reasonably close to the ones from the MR simulations, with relative deviations of H and R below 4%, and similar polymer densities in the central region of the films and droplets. However, the polymer–solvent interfaces are slightly broader in the CG simulations, which is also reflected by the smaller surface tension measured in the CG simulations, that is γ = 0.31 ± 0.01k_B T/σ² vs γ = 1.42k_B T/σ² for the planar films. Further, the peaks in the density profiles near the polymer–solvent interfaces are slightly more pronounced and more narrow in the MR simulations compared to the CG simulations. This difference can be understood if one considers the oblate ellipsoidal shape of the polymers at the interfaces in the MR simulations (see figure 2(b)). This shape anisotropy near the polymer–solvent interface can, however, not be captured by our CG model, where the polymers are modeled as soft spheres.

Zoom In Zoom Out Reset image size

In addition, we computed the force profiles along the z-axis in the films and in the radial direction r for the droplets (see figure 10). The force profiles were calculated as the mean forces of particles at the given z or r position with data taken every 100τ. In both the MR and CG simulations, the net force is essentially zero in the central region of the film and droplet, due to the symmetry of the systems. Polymers near the surface of the film/droplet experience attractive forces toward the center of the film/droplet. These cohesive forces are stronger and act on a much narrower region in the MR simulations compared to the CG simulations. If one decomposes the net force into the contributions from the two- and many-body interactions, one sees that ${U}_{0}^{2\text{b}}$ and ${U}_{0}^{\text{LDP}}$ have opposite signs and are each significantly larger than the resulting net force. Both forces also extend much further into the central region of the films and droplets than expected from the width of the bumps at the extremities of the force profiles. Evidently, the stability of the films and droplets is the result of a delicate balance between two- and many-body interactions in our CG model.

Zoom In Zoom Out Reset image size