Equilibration of polymer melts containing highly entangled long polymer chains in confinement or with free surfaces is a challenge for computer simulations. We approach this problem by first studying polymer melts based on the soft-sphere coarse-grained model confined between two walls with periodic boundary conditions in two directions parallel to the walls. Then, we insert the microscopic details of the underlying bead-spring model. Tuning the strength of the wall potential, the monomer density of confined polymer melts in equilibrium is kept at the bulk density even near the walls. In a weak confining regime, we observe the same conformational properties of chains as in the bulk melt showing that our confined polymer melts have reached their equilibrated state. Our methodology provides an efficient way of equilibrating large polymer films with different thicknesses and is not confined to a specific underlying microscopic model. Switching off the wall potential in the direction perpendicular to the walls enables to study free-standing highly entangled polymer films or polymer films with one supporting substrate.

## I. INTRODUCTION

Polymer confinement plays an important role for many aspects of adhesion, wetting, lubrication, and friction of complex fluids from both theoretical and technological points of view. For more than two decades, both theoretical and experimental works^{1–8} have shown that the dynamic and structural properties of polymers subject to confinement may deviate from that in the bulk as the interaction between polymers and confining surfaces becomes non-negligible. For example, the conformations of polymer chains in a melt near the wall in the direction perpendicular to the wall shrink remarkably compared to bulk chains, while they only extend slightly parallel to the wall.^{3,9,10} Long chain mobility in confined melts is also affected by entanglement effects, while, in turn, the conformational deviations from bulk chains have influence on the distribution of entanglements.^{11–15} Furthermore, many studies have also focused on the dependency between the glass transition temperature *T*_{g} and the nature of the confinement effects.^{16–18} Therefore, it is important to understand the mechanical properties of confined polymer melts and how confinement impacts both viscous and elastic properties of amorphous polymer films with different surface substrates or even free surfaces.

Computer simulations provide a powerful method to mimic the behavior of polymers under well-defined external conditions covering the range from atomic to coarse-grained (CG) scales.^{19,20} However, the cost of computing time rises dramatically as the size and complexity of systems increase. To avoid this, applying appropriate coarse-grained models that keep the global thermodynamic properties and the local mechanical and chemical properties is still an important subject.^{21–30} One of the successful monomer-based models, namely, the bead-spring (BS) model^{22,23} together with an additional bond-bending potential,^{31–34} has been successfully employed to provide a better understanding of generic scaling properties of polymer melts in bulk. For such a model, static and dynamic properties of highly entangled polymer melts in bulk have been extensively studied in our previous work.^{35,36} We have verified the crossover scaling behavior of the mean square displacement of monomers between several characteristic time scales as predicted by the Rouse model^{37} and reptation theory^{38–41} over several orders of magnitude in the time. For weakly semiflexible polymer chains of sizes of up to *N* = 2000 monomers, we have also confirmed that they behave as ideal chains to a very good approximation. For these chains, the entanglement length *N*_{e} = 28 in the melt, estimated through the primitive path analysis and confirmed by the plateau modulus,^{34,35,42} results in polymer chains of *N* = 2000 ≈ 72*N*_{e}. Thus, we here focus on this model and a related variant for the study of the equilibration of polymer melts confined between two repulsive walls as supporting films, and free-standing films after walls are removed. Each film contains *n*_{c} = 1000 chains of *N* = 2000 monomers at the bulk melt density *ρ*_{0} = 0.85*σ*^{−3}. Directly equilibrating such large and highly entangled chains in bulk or confinement is not feasible within reasonably accessible computing time.

A novel and very efficient methodology has recently been developed^{28,43} for equilibrating large and highly entangled polymer melts in bulk described by the bead-spring model.^{22,23} Through a hierarchical backmapping of CG chains described by the soft-sphere CG model^{27,28} from low resolution to high resolution and a reinserting of microscopic details of bead-spring chains, finally, highly entangled polymer melts in bulk are equilibrated by molecular dynamics (MD) simulations using the package ESPResSO++.^{44,45} To first order, the required computing time depends only on the overall system size and becomes independent of the chain length. Similar methodologies have also been used to equilibrate high-molecular-weight polymer blends^{46} and polystyrene melts.^{47} In this paper, we extend the application of the soft-sphere approach to confined polymer melts and subsequently free-standing films. As polymer chains are described at a lower resolution, the number of degrees of freedom becomes smaller. Here, we adapt this hierarchical approach to equilibrate polymer melts confined between two walls in detail. Moreover, we apply our newly developed, related model^{48,49} to prepare polymer films with one or two free surfaces. Different from Refs. 27 and 43, we take the bending elasticity of bead-spring chains in a bulk melt into account for the parameterization of the soft-sphere CG model. Namely, the underlying microscopic bead-spring chains are weakly semiflexible (bending constant *k*_{θ} = 1.5*ε*) instead of fully flexible (*k*_{θ} = 0.0*ε*).

The outline of this paper is as follows: In Sec. II, we introduce the main features of the microscopic bead-spring model and soft-sphere coarse-grained model used for studying the confined polymer melts. The application of the soft-sphere CG model for confined coarse-graining melts and conformational properties of fully equilibrated confined CG melts are addressed in Sec. III. In Sec. IV, we reinsert the microscopic details of confined CG melts and discuss the equilibration procedures. In Sec. V, we show how to prepare films with one or two free surfaces by switching to another variant of the bead-spring model. Finally, our conclusion is given in Sec. VI.

## II. MODELS

### A. Generic microscopic bead-spring models

In the microscopic bead-spring model,^{22,23} all monomers at a distance *r* interact via a shifted Lennard-Jones (LJ) potential *U*_{LJ}(*r*),

with

where *ε* is the energy strength of the pairwise interaction and *r*_{cut} is the cutoff in the minimum of the potential such that force and potential are zero at *r*_{cut} = 2^{1/6}*σ*. These LJ units also provide a natural time definition via $\tau =\sigma m/\epsilon $, with *m* = 1 being the mass of the monomers. The temperature is set to *T* = 1.0*ε*/*k*_{B}, with *k*_{B} being the Boltzmann factor, which is set to one. Along the backbone of the chains, a finitely extensible nonlinear elastic (FENE) binding potential *U*_{FENE}(*r*) is added,

where *k* = 30*ε*/*σ*^{2} is the force constant and *R*_{0} = 1.5*σ* is the maximum value of bond length. For controlling the bending elasticity, i.e., chain stiffness, the standard bond-bending potential^{31–34} is given by

with the bond angle *θ* defined by $\theta =cos\u22121b\u2192j\u22c5b\u2192j+1|b\u2192j\Vert b\u2192j+1|$, where $b\u2192j=r\u2192j+1\u2212r\u2192j$ is the bond vector between the (*j* + 1)th monomer and the *j*th monomer along the identical chain. The bending factor *k*_{θ} is set to 1.5*ε*, and the melt density is set to the widely used value of 0.85*σ*^{−3} throughout the whole paper.

For studying polymer melts under confinement, we first consider the simpler example of polymer melts that are confined between two planar, structureless repulsive walls. The walls placed at *z* = 0 and *z* = *L*_{z} are described by the 10-4 Lennard-Jones planar wall potential,^{3,50}

with

Here, *ε*_{w} is the interaction strength between monomers and the walls, and *z* and *L*_{z} − *z* are the vertical distances of a monomer from the two walls, respectively.

For preparing free-standing polymer films with one or two free surfaces, i.e., mimicking free-standing films in vacuum or later on studying them in presence of a liquid or a gas, we first have to stabilize the system at zero pressure, particularly, in the direction perpendicular to the walls. Only then, we can switch off the wall potential and prevent system instability. Therefore, a short-range attractive potential to reduce the pressure to zero is added^{48,49} with an additional shift term,

such that *U*_{ATT}(*r*) = *U*_{LJ}(*r*) at *r* = *r*_{cut}. Here, *α* = 0.5145*ε* is the strength parameter, and $rca=2rcut\u22481.5874\sigma $ is the upper cutoff such that it has zero force at *r* = *r*_{cut} and $r=rca$. Note that this additional potential does not alter the characteristic conformations at *T* = 1*ε*/*k*_{B} so that we can switch between these different models as needed.^{22,23,48,49,51}

Note that using $UBEND(old)(\theta )$, bead-spring chains tend to stretch out as the temperature decreases. To avoid such an artificial chain stretching, Eq. (4) can be replaced by^{48,49}

if one is interested in studying polymer melts under cooling. The fitting parameters *a*_{θ} = 4.5*ε* and *b*_{θ} = 1.5 are determined so that the local conformations of chains remain essentially unchanged compared to those with *k*_{θ} = 1.5*ε* using Eq. (4) at temperature *T* = 1.0*ε*/*k*_{B}.

### B. Soft-sphere coarse-grained model

In the soft-sphere approach,^{27} each bead-spring polymer chain in a melt is represented by a chain of *N*_{CG} = *N*/*N*_{b} fluctuating soft spheres in a CG view, as shown in Fig. 1. The coordinate of the center and the radius of the *s*th sphere in the *i*th chain, $r\u2192i(s)$ and *σ*_{i}(*s*), are described by

and

respectively. Here, $r\u2192i,k$ is the coordinate of the *k*th monomer in chain *i*, and $R\u2192CM,i(s)$ and *R*_{g,i}(*s*) are the center of mass (CM) and the radius of gyration of *N*_{b} monomers in its subchains, respectively. Since we extend the application of the soft-sphere CG model to semiflexible chains and follow the slightly different strategies of Refs. 28 and 43, we list the details of soft-sphere potentials as follows. The spheres are connected by a harmonic bond potential

and an angular potential

with

where $d\u2192i(s)=r\u2192i(s+1)\u2212r\u2192i(s)$ is the bond vector and *b*_{CG} denotes the effective bond length. The radius of the *s*th sphere in chain *i* and its fluctuation are controlled by the following two potentials:^{27,28}

and

respectively. The non-bonded interactions between any two different spheres due to excluded volume interactions according to Flory’s theory^{21,38,52} are taken care of by the potential

Here, the parameters *k*_{bend}, *a*_{1}, *a*_{2}, *a*_{3}, *b*_{CG}, and *ε*_{1} depending on *N*_{b} are determined by a numerical approach via curve fitting.

Similarly, as shown in Eq. (5), the two soft repulsive planar walls located at *z* = 0 and *z* = *L*_{z} depending on the radius of each sphere are given by

with

where $\epsilon w(CG)$ is the interaction strength between soft spheres and the walls and *z* and *L*_{z} − *z* are the vertical distances from the two walls, respectively.

For the parameterization of the soft-sphere CG model, we take 15 independent and fully equilibrated bulk polymer melts of bead-spring polymer chains with *k*_{θ} = 1.5*ε* obtained from the previous works^{35,42,43} as our reference systems. Using Eqs. (9) and (10) with *N*_{b} = 25, each melt containing *n*_{c} = 1000 bead-spring chains of *N* = 2000 = *N*_{CG}*N*_{b} monomers in a cubic simulation box of size *V* = *L*^{3} (*L* ≈ 133*σ* with periodic boundary conditions in *x*-, *y*-, and *z*- directions) at the melt density *ρ*_{0} = 0.85*σ*^{−3} is mapped onto a CG melt containing *n*_{c} = 1000 soft-sphere chains of *N*_{CG} = 80 spheres at the CG melt density $\rho 0(CG)(NCG)=\rho 0NCG/N=0.034\sigma \u22123$. The parameters *a*_{1} = 6.3444 × 10^{−3}*σ*^{6}, *a*_{2} = 4.1674*σ*^{−2}, *k*_{bend} = 1.3229*ε*, *a*_{3} = 22.0*εσ*^{3}, *b*_{CG} = 6.68*σε*^{−1/2}, and *ε*_{1} = 290.0*εσ*^{3} are determined^{27,28} such that the average conformational properties of the reference melt systems in a CG representation are reproduced by fully equilibrated CG melts of soft-sphere chains. Quantitatively, the conformational properties are characterized by the probability distributions of the radius of soft spheres, *P*(*σ*, *N*_{b}), the bond length connecting two successive soft spheres, *P*(*d*, *N*_{b}), and the bond angle between two successive bonds, *P*(*θ*, *N*_{b}), the average mean square internal distance between the *j*th soft sphere and the (*j* + *s*)th soft sphere along the identical chain, ⟨*R*^{2}(*s*, *N*_{b})⟩, and the pair distribution of all pairs of soft spheres, *g*(*r*, *N*_{b}) (see Figs. 2 and 3 discussed in Sec. III).

Since the excluded volume effect between *N*_{b} monomers in each subchain is ignored in the parameterization of the soft-sphere approach and self-entanglements on this length scale are negligible as *N*_{b} < *N*_{e}, subchains behave as ideal chains (alternatively, one could include excluded volume for semi dilute solutions via a Flory term). For *N*_{b} = 25(<*N*_{e} = 28), the chain to a good approximation can be described as a Gaussian chain, while this is not the case for significantly shorter chains. In this case, we can simplify several steps of hierarchical backmapping^{43} to only one step of fine-graining to introduce microscopic details of subchains once a CG melt reaches its equilibrated state.

## III. EQUILIBRATION OF SOFT-SPHERE CHAINS IN A CONFINED CG MELT

In this section, we extend the application of the soft-sphere CG model for polymer melts in bulk to polymer melts confined between two repulsive walls [Eq. (16)], first focusing on the CG melt containing *n*_{c} = 1000 chains of *N*_{CG} = 80 spheres at the CG bulk melt density $\rho 0(CG)=0.034\sigma \u22123$. For the comparison to our reference systems, we set the distance between two walls compatible with the bulk melt. Thus, we locate two walls at *z* = 0*σ* and *z* = *L*_{z} ≈ 133*σ* while keeping the periodic boundary conditions along the *x*- and *y*- directions with the lateral linear dimensions *L*_{x} = *L*_{y} = 133*σ*. Of course, one can adjust *L*_{z} and extend/reduce *L*_{x} and *L*_{y} for keeping the bulk melt density as needed.

The initial configurations of soft-sphere chains in terms of {*σ*_{i}(*s*), *d*_{i}(*s*), *θ*_{i}(*s*)} are randomly generated according to their corresponding Boltzmann weights $exp[\u2212\beta Usphere(CG)(\sigma i(s))]$, $exp[\u2212\beta Ubond(CG)(di(s))]$, and $exp[\u2212\beta Uang(CG)(\theta i(s))]$, respectively, where *β* = 1/(*k*_{B}*T*) = 1.0*ε*^{−1}. Additionally, we set 1*σ* < *σ*_{i}(*s*) < *σ*_{max} = 8*σ* and 0*σ* < *d*_{i}(*s*) < *d*_{max} = 21*σ* but restrict the coordinates of centers of spheres, $r\u2192i(s)$, satisfying the condition $\sigma i(s)<r\u2192i(s)<(Lz\u2212\sigma i(s))$. It is computationally more efficient to perform Monte Carlo simulations to equilibrate confined CG melts. Similar to Ref. 28, our simulation algorithm including three types of MC moves at each step is as follows: (i) For a local move, one of the *n*_{c}*N*_{CG} spheres is randomly selected, e.g., the *s*th sphere in the *i*th chain, the sphere of radius *σ*_{i}(*s*) at $r\u2192i(s)=(ri,x(s),ri,y(s),ri,z(s))$ is allowed to move within the range −*σ*_{i}(*s*) < Δ*r*_{i,x}(*s*), Δ*r*_{i,y}(*s*), Δ*r*_{i,z}(*s*) < *σ*_{i}(*s*). The trial move is therefore accepted if $exp[\u2212\beta (\Delta Unb(CG)+\Delta Ubond(CG)+\Delta Uang(CG)+\Delta Uwall(CG))]>\eta $, where *η* is a random number and *η* ∈ [0, 1). (ii) For a snake-slithering move, one end of the *n*_{c} chains is randomly selected, and $\sigma i(new)(s)$, $di(new)(s)$, and $\theta i(new)(s)$ of the selected sphere are randomly generated according to their corresponding Boltzmann weights, respectively. The trial move is accepted if $exp[\u2212\beta (\Delta Unb(CG)+\Delta Uself(CG)+\Delta Uwall(CG))]>\eta $. (iii) For a trial change in the sphere size, one of the *n*_{c}*N*_{CG} spheres is randomly selected, and $\sigma i(new)(s)$ of the selected sphere is randomly generated according to its Boltzmann weight. Finally, the trial move is accepted if $exp[\u2212\beta (\Delta Unb(CG)+\Delta Uself(CG))]>\eta $. A cutoff at $rcut(CG)=20\sigma $ for calculating the non-bonded interactions between two different spheres, $Unb(CG)(ri(s),\sigma i(s),ri\u2032(s\u2032),\sigma i\u2032(s\u2032))$, is also introduced^{27} since the contributions for *r* > 20*σ* are negligible. Nevertheless, there is no influence on measurements of any physical observable while it speeds up the simulations by a factor of four. Applying a linked-cell algorithm with the cell size *L*_{c} = 2.66*σ* (*L*_{x}/*L*_{c} = *L*_{y}/*L*_{c} is very close to an integer), smaller than the cutoff value $rcut(CG)$, speeds up the simulation even more by an additional factor of 2.5, i.e., all together, it speeds up by a factor of ten. It takes about 10 h central processing unit (CPU) time on an Intel 3.60 GHz PC for a confined CG melt to reach its equilibrated state (after 2 × 10^{7} MC steps are performed). The acceptance ratio is about 73% for a trial change of the sphere size, 45% for a local move, and 41% for a snake-slithering move.

Choosing the wall strength $\epsilon w(CG)\u2248O(0.1\u22121)\epsilon $, there is no detectable influence on the probability distributions *P*(*σ*, *N*_{b}), *P*(*d*, *N*_{b}), and *P*(*θ*) comparing to that for an equilibrated CG melt in bulk, as shown in Fig. 2. Since soft spheres are allowed to penetrate each other in CG melts, we observe that the distributions *P*(*σ*, *N*_{b}) and *P*(*d*, *N*_{b}) are slightly narrower for both equilibrated CG melts in bulk and in confinement than that for the reference systems while the average values of radius and bond length remain the same. We have also compared the estimates of ⟨*R*^{2}(*s*, *N*_{b})⟩ and *g*(*r*, *N*_{b}) for an equilibrated confined CG melt to an equilibrated CG melt in bulk and the reference data in Fig. 3. The curve of ⟨*R*^{2}(*s*, *N*_{b})⟩ taken the average over all *n*_{c} = 1000 chains for a confined CG melt deviates from the bulk behavior for *s* > 10 due to the confinement effect, while for *s* < 6, ⟨*R*^{2}(*s*, *N*_{b})⟩ is a bit smaller compared to the bulk value. It is due to the artifact of the soft-sphere CG model, where the excluded volume effect within the size of spheres is not considered. After monomers are reinserted into soft-sphere chains, local excluded volume and the corresponding correlation hole effect automatically correct for these deviations. However, the discrepancy for *s* < 6 is still within fluctuations observed in bulk. When monomers are reinserted into soft-sphere chains, the estimate of *g*(*r*, *N*_{b}) starts to increase at *r* ≈ 3*σ* and then decrease at *r* ≈ 5*σ*. It indicates that near the walls, the distance between any two spheres decreases due to the confinement effect.

To investigate the confinement effect on packing and conformations of a polymer melt, we determine the soft-sphere density profile between two walls as follows:

Figure 4 shows that the soft-sphere density profiles with bin size 1*σ* for three different values of the interaction strength between the soft spheres and walls, $\epsilon w(CG)/\epsilon =1.0$, 0.5, and 0.1, are the same within small fluctuation. The bulk melt density persists, i.e., $\rho (CG))(z,Nb)=\rho 0(CG)$, between *z* = 5 and *z* = *L*_{z} − 5. *ρ*^{(CG)}(*z*, *N*_{b}) increases and reaches a maximum value at *z* = 3*σ* and then approaches zero next to the walls. It indicates that the confinement effect is weak for spheres sitting in the middle regime between two walls. The change in *ρ*^{(CG)}(*z*, *N*_{b}) is related to *P*(*σ*, *N*_{b}) since *P*(*σ*, *N*_{b}) has its maximum at *z* = ⟨*σ*_{i}(*s*)⟩ ≈ 3.06*σ* (see Fig. 2). The two components of the mean square radius of gyration depending on the *z*-component of CMs of chains are defined as follows:

where $r\u2192i,\alpha =\Vert (s)=r\u2192i,x(s)+r\u2192i,y(s)$ and $r\u2192i,\alpha =\u22a5(s)=r\u2192i,z(s)$. Figure 4(b) shows that the linear dimensions of confined chains having their CM in the regime $0.3Lz\u2264r\u2192i,z(CM)\u22640.7Lz$ are the same as in a bulk melt within fluctuation. For polymer chains of *N* = 2000 monomers in a bulk melt, the mean square radius of gyration is $\u27e8Rg2\u27e90=\u27e8Rg,\Vert 2\u27e9+\u27e8Rg,\u22a52\u27e9\u2248909\sigma 2$. With a decrease in the distance *z* from the walls, $\u27e8Rg,\Vert 2(z,Nb)\u27e9$ increases moderately with larger fluctuations, while $\u27e8Rg,\u22a52(z,Nb)\u27e9$ decreases gradually even already in the regime where the monomer density is $\rho (CG)(z,NCG)\u2248\rho 0(CG)$. Note that none of chains has its CM next to the wall. From the results shown in Fig. 4(b), we should expect that ⟨*R*^{2}(*s*, *N*_{b})⟩ and *g*(*r*, *N*_{b}) follow the bulk behavior if we count those chains sitting in the middle regime ($0.3Lz\u2264r\u2192i,z(CM)\u22640.7Lz$) between two walls. It is indeed seen in Figs. 3(b) and 3(d). Similar behavior has been observed for shorter chains confined between two walls.^{3,9,10,53}

## IV. EQUILIBRATING BEAD-SPRING CHAINS IN A CONFINED MELT

### A. Backmapping procedure

After equilibration of the CG melt, we apply the similar backmapping strategy developed in Ref. 43 to reinsert the microscopic details of the bead-spring model described in Sec. II A (see Fig. 5). In this strategy, two monomers along the chains are bonded via the FENE potential [Eq. (3)] and the shifted LJ potential [Eq. (1)]. The non-bonded and bond-bending interactions are excluded at this step. The confinement effect is introduced by the soft repulsive wall potential [Eq. (5)]. Each soft-sphere CG chain of *N*_{CG} = 80 spheres is now replaced by a bead-spring chain of *N* = 2000 monomers. To preserve the relationship between a soft sphere and a subchain of *N*_{b} = 25 monomers given in Eqs. (9) and (10), two pseudopotentials^{43} for the *s*th soft sphere in the *i*th chain are implemented as follows:

and

where *k*_{cm} and *k*_{g} determine the coupling strength. The forces derived from these two potentials can drive the center of mass and the radius of gyration of each subchain to the center and the radius of the corresponding soft sphere, respectively. Namely, each bead-spring chain is then sitting on top of its corresponding soft-sphere CG chain (see Fig. 1). During this backmapping procedure, it is more practical to perform MD simulations in the NVT ensemble with a weak coupling Langevin thermostat at *T* = 1.0*ε*/*k*_{B} by setting the friction constant Γ = 0.5*τ*^{−1}. Choosing *k*_{cm} = 50.0*ε* and *k*_{g} = 5.0*ε*, the integration time step is set to Δ*t* = 0.005*τ*. At this stage, all 1000 soft-sphere CG chains can be mapped into 1000 bead-spring chains confined between two walls simultaneously since there is no interaction between different chains. Snapshots of the configurations of the fully equilibrated confined CG melt and the backmapped confined melt of bead-spring chains are shown in Fig. 5. Here, the strength of the wall potential *U*_{wall}(*z*) is set to *ε*_{w} = 0.005*ε*. To keep the bulk melt density *ρ*_{0} = 0.85*σ*^{−3} in the middle regime between two walls (the weak confinement regime) in a microscopic representation, we set *L*_{z} = 134*σ* instead of *L*_{z} = 133*σ* taking the repulsive potential of the wall, which is a steep but smooth function, into account. The reinsertion MD time is about 80*τ* and the CPU time is about 2.8 h on a single processor in an Intel 3.60 GHz PC.

### B. Equilibration procedure

In the next step, the full excluded volume interaction, as listed in Sec. II A, has to be introduced. To avoid the “explosion” of the system due to the overlap of monomers, we have to switch on the excluded volume interactions between the non-bonded pairs of monomers in a quasi-static way (slow push-off).^{42,51} Therefore, the shifted LJ potential for each non-bonded pair of monomers at a distance *r*, [Eq. (1)], is first replaced by a force-capped LJ potential

where *r*_{fc} is an adjustable cutoff distance in this warm-up procedure. $rfc=rfco$ decreases monotonically at each cycle from 2^{1/6}*σ* to 0.8*σ* in the warm-up procedure for the non-bonded pairs except the next-nearest neighboring (nn) pairs along identical chains. For the nn pairs, $rfc=rfcnn$ is set to 1.1*σ* initially and tuned according to the following cost function:^{42}

at the end of each cycle in the warm-up procedure with the restriction that $0.8\sigma \u2264rfcnn\u22641.1\sigma $. Master curve refers to fully equilibrated polymer melts in bulk (the reference systems). Here, we have assumed that the mean square internal distance, ⟨*R*^{2}(*s*)⟩, for the confined chains in a melt finally should coincident with the master curve for *s* < 50 (up to two soft spheres in a CG representation) at least. For a polymer melt under strong confinement effect, the bulk behavior may no longer be valid even for *s* < 50 ≈ 2*N*_{e}. However, it is more important to first correct chain distortion due to the absence of excluded volume interactions between monomers [see Fig. 6(b)]. Once we remove the criterion of the cost function, confined chains will relax very fast on a short length scale dominated by both the entanglement and confinement effects. Note that this final equilibration step only affects subchain lengths of up to the order of *N*_{e}. In this process, the nn-excluded volume is adjusted according to the value of *C*, namely, reduced (increase $rfcnn$ by 0.01*σ*) if *C* < 0 and enhanced (decrease $rfcnn$ by 0.01*σ*) if *C* > 0. If |*C*| < 0.0001*σ*^{2}, $rfcnn$ remains unchanged.

We perform MD simulations in the NVT ensemble with a Langevin thermostat at the temperature *T* = 1.0*ε*/*k*_{B} using the package ESPResSo++^{44,45} for equilibrating a confined polymer melt containing *n*_{c} = 1000 chains of *N* = 2000 monomers under three procedures as follows: (a) In the warm-up procedure, 120 cycles of 1.5 × 10^{5} MD steps per cycle in the first 80 cycles and 5 × 10^{4} MD steps per cycle in the rest 40 cycles are performed with a larger friction constant Γ = 1.0*τ*^{−1} and a small time step Δ*t* = 0.0002*τ*. Differently from Ref. 42, more MD steps and a slower rate of decreasing the cutoff value of $rfco$ for non-bonded monomer pairs except nn pairs are required for the confined polymer melt system due to the competition between the excluded volume effect and the confinement effect. In the first 100 cycles, $rfcnn$ is updated at each cycle associated with the cost function. In the rest 20 cycles, $rfcnn$ decreases by 0.01*σ* before reaching the minimum value 0.8*σ* [see Fig. 6(a)]. (b) In the relaxation procedure, the shifted LJ potential *U*_{LJ}(*r*) is restored. We first perform 10^{5} MD steps with Γ = 0.5*τ*^{−1} and Δ*t* = 0.001*τ* and then another 2 × 10^{6} MD steps with Δ*t* = 0.005*τ* to ensure that the confined melt reaches its equilibrated state. Afterward, we can set the time step to its standard value, Δ*t* = 0.01*τ*, for the further study of the confined polymer melt in equilibrium.

During warming up the confined polymer melt, the local monomer density profile *ρ*(*z*) near the walls varies with the interaction strength *ε*_{w} between monomers and walls, as shown in Fig. 6(a). For *ε*_{w} = 0.005*ε*, the bulk melt density *ρ*_{0} is conserved all the way up to the wall, where the bin size is set to 0.5*σ*. At the same time, after the confined polymer melt is warmed up, the curve of ⟨*R*^{2}(*s*)⟩ coincides with the master curve for *s* < 50, as shown in Fig. 6(b). A typical variation of cutoff distances, $rfcnn$ and $rfco$, for the non-bonded pairs of monomers in the warm-up procedure are shown in Fig. 6(c). At the end of the warm-up procedure, *U*_{FC-LJ}(*r*) approaches to *U*_{LJ}(*r*). Thus, there is no problem to switch back to *U*_{LJ}(*r*) for the further relaxation of the confined polymer melt.

First, we examine the difference between the equilibrated confined melt and the reference bulk melt, as shown in Fig. 7. We compare the whole system in terms of the mean square internal distance ⟨*R*^{2}(*s*)⟩, the single chain structure factor *S*_{c}(*q*), and the collective structure $S\u2225$(*q*) in the direction parallel to the walls. We see that the curve of ⟨*R*^{2}(*s*)⟩ estimated only for those chains having their CMs satisfying $0.3Lz\u2264r\u2192i,z(CM)\u22640.7Lz$ follows the master curve, while the curve obtained by taking the average over all 1000 chains starts to deviate from bulk behavior at *s* = 250. To detect any anisotropy of chain conformations under confinement, we distinguish *S*_{c,⊥}(*q* = *q*_{z}), where the wave vector $q\u2192$ is oriented in the *z*-direction perpendicular to the wall and $Sc,\Vert (q=(qx2+qy2)1/2)$. $Sc,\u2225$(*q*) follows the same chain structure as in the bulk melt, as shown in Fig. 7(b). The shift of *S*_{c,⊥}(*q*) toward to a slightly larger value of *q* from the bulk curve also indicates that the estimate of $\u27e8Rg,\u22a52\u27e9$ for all 1000 chains is smaller. Nevertheless, chains still behave as ideal chains. On comparing the collective structure factor of the whole melt between the confined polymer and the bulk melt in the parallel direction to the walls, we see that there is no difference as the distance between two walls is compatible to the bulk melt [see Fig. 7(c)]. The local packing of monomers characterized by the pair distribution *g*(*r*) for both inter- and intra-pairs of monomers is also in perfect agreement with the bulk melt, as shown in Fig. 7(d).

Finally, we go back to the outset and analyze the conformational properties of the confined equilibrated melt in a CG representation by mapping bead-spring chains to soft-sphere chains using Eqs. (9) and (10). As shown in Figs. 3 and 4(b), all the results obtained by taking the average over 14 configurations within 6*τ*_{e} are consistent with the data obtained from the MC simulations of a confined CG melt within fluctuation. Since in the microscopic model, the excluded volume interactions between monomers are properly taken into account, the curves of ⟨*R*^{2}(*s*, *N*_{b})⟩ *g*(*r*, *N*_{b}) at short length scales do not deviate from the curves for the reference systems in a CG representation. This shows that the confined polymer melt indeed reaches the equilibrium on all length scales.

## V. PREPARATION OF SUPPORTED AND FREE-STANDING FILMS AT ZERO PRESSURE

In order to study free-standing or supported films at pressure *P* = 0.0*ε*/*σ*^{3}, stabilizing non-bonded attractive monomer monomer interactions are required. For this, we switch to our recently developed CG model^{48,49} based on the bead-spring model^{22,23} by simply turning on the attractive potential *U*_{ATT}(*r*) [Eq. (7)] for non-bonded monomer pairs and replacing the bending potential $UBEND(old)(\theta )$ with *U*_{BEND}(*θ*) [Eq. (8)]. This choice of interaction has the additional advantage that it allows us to study glassy films as well. Starting from a fully equilibrated polymer melt confined between two walls obtained in the last section, we perform MD simulations in the NVT ensemble with a Langevin thermostat at the temperature *T* = 1.0*ε*/*k*_{B} using the package ESPResSo++,^{44,45} keeping the short range repulsion from the walls. Figure 8(a) shows that the three diagonal terms of the pressure tensor *P*_{αβ}(*t*) first drop from 4.9*ε*/*σ*^{3} (5.0*ε*/*σ*^{3} for bulk melts) to 1.4*ε*/*σ*^{3} right after *U*_{ATT}(*r*) is switched on and then to 0.1*ε*/*σ*^{3} in a very short time about 20*τ*. We further relax the confined polymer film for 30 000*τ* ≈ 13*τ*_{e}, with *τ*_{e} ≈ 2266*τ* being the entanglement time,^{35} by performing MD simulations in the NPT ensemble (Hoover barostat with Langevin thermostat^{54,55} implemented in ESPResSo++^{44,45}) at temperature *T* = 1.0*ε*/*k*_{B} and pressure *P* = (*P*_{xx} + *P*_{yy} + *P*_{zz})/3 = 0.0*ε*/*σ*^{3} to finally adjust the pressure from 0.1*ε*/*σ*^{3} to 0.0*ε*/*σ*^{3}. Under this circumstance, an equilibrated free-standing film is generated after removing two walls by turning off the wall potential at *z* = 0*σ* and *z* = *L*_{z} = 134*σ*. If we only remove one of the walls, we get a polymer film with one supporting substrate, where one, of course, can introduce appropriate adhesion interactions.

The overall conformations of all chains and inner chains ($0.3Lz\u2264ri,z(CM)\u22640.7Lz$) as characterized by ⟨*R*^{2}(*s*)⟩ between two walls for a confined polymer melt based on two variants of bead-spring model and after relaxing for 30 000*τ* (NPT MD simulations) are preserved within fluctuation, as shown in Figs. 8(b) and 8(c). The monomer density profile *ρ*(*z*) in the direction perpendicular to the wall compared to the density in the interior of confined melt is also preserved. as shown in Fig. 8(d). The thickness of films, *D*, is determined according to the concept of the Gibbs dividing surface that has been applied to identify the interface between two different phases,^{56–58} e.g., liquid and vapor, and polymer and vacuum, based on the density profile in the direction perpendicular to the interfaces. The locations of the Gibbs dividing surfaces (planar surfaces) corresponding to the upper and lower bounds of films,

and

are obtained by the requirement of equal areas that

and

respectively. Here, *z*_{min} and *z*_{max} are the two limits where *ρ*(*z*) approaches to zero, and *z*_{c} = (*z*_{min} + *z*_{max})/2. The mean monomer density is given by

where the value of Δ*z* is chosen such that the monomer density profile *ρ*(*z*) in the interval [*z*_{c} −Δ*z*, *z*_{c} + Δ*z*] reaches a plateau value within small fluctuation. Choosing Δ = 2.5*σ*, the thickness of confined film, $D=zG(upper)\u2212zG(lower)$, at *T* = 1.0*ε*/*k*_{B} reduces from 133.4*σ* at *P* = 4.9*ε*/*σ*^{3} to 131.6*σ* at *P* = 0.1*ε*/*σ*^{3} due to the short-range attractive interaction between non-bonded monomers, and finally is stabilized at 130.3*σ* at *P* = 0.0*ε*/*σ*^{3} where the lateral dimensions of film increase slightly from *L*_{x} = *L*_{y} ≈ 133.0*σ* to *L*_{x} = *L*_{y} ≈ 134.0*σ*.

To further relax the free-standing film, we have also performed MD simulations in the NVT ensemble at the temperature *T* = 1.0*ε*/*k*_{B} for 30 000*τ*, where the resulting pressure is *P* = 0.0*ε*/*σ*^{3}. The perpendicular simulation box size is set to *L*_{z} = 194*σ* for preventing any interaction between monomers and the lateral surfaces of the box. Snapshots of the configurations of confined and free-standing films after relaxing for 30 000*τ* are shown in Fig. 9. The estimate of ⟨*R*^{2}(*s*)⟩ for all chains and chains in the middle part of the free-standing film, and *ρ*(*z*) with bin size 0.25*σ* are shown in Figs. 8(b), 8(c), and 10, respectively. We see that after relaxing polymer chains in a free-standing film for 30 000*τ* ≈ 13*τ*_{e}, the thickness of the free-standing film, 130.5*σ*, is still compatible with that of the confined polymer film, 130.3*σ*. The average conformations of all chains and inner chains also remain the same, while the tails of the monomer density profile become longer, indicating that the surface becomes a bit rougher.

Finally, the density of chain ends near surfaces for both confined and free-standing films at *T* = 1.0*ε*/*k*_{B} and *P* = 0.0*ε*/*σ*^{3} is examined. For this, we compare the normalized density of all monomers, *ϕ*(*z*) = *ρ*(*z*)/*ρ*_{0}, the density of end monomers, *ϕ*_{e}(*z*) = *ρ*_{e}(*z*)/(2*ρ*_{0}/*N*), and the relative excess of end monomers, $\varphi e(excess)(z)=[\rho e(z)\u2212(2/N)\rho (z)]/(2\rho 0/N)$, along the perpendicular direction of interfaces, averaged over 30 configurations within 6*τ*_{e}, which are shown in Fig. 11. To illustrate the dependency on the chain length or bead size, we have also mapped the bead-spring chains of *N* = 2000 monomers of size *σ* = 1 onto underlying soft-sphere chains of *N*_{CG} = 80 spheres of approximate size 6*σ* [average radius ⟨*R*_{g}(*N*_{b} = 25)⟩ ≈ 3.0*σ*]. For films, we observe a weak enrichment of end monomers at the surfaces due to a potential gain of entropy.^{59,60} The related depletion zone for chain ends near the surface, as predicted by self-consistent field theories,^{59,60} however, turns out to be too small to be resolved within the fluctuations of our data. The coarse-graining slightly smears out this enrichment effect due to the overlap of soft spheres. For the free-standing film, the enrichment effect of end monomers near the interfaces is slightly less pronounced, while at the same time, the interface widens. This indicates a weak roughening of the free surface compared to the confined film. Nevertheless, the indicated $\varphi e(excess)(z)$ near the surfaces in both cases is very small and levels off on the scale of the typical bulk density correlation length given, e.g., in Fig. 7.

## VI. CONCLUSION

In this paper, we have developed an efficient methodology to equilibrate long chain polymer films and applied this method to a polymer film where 1000 chains of 2000 ≈ 72*N*_{e} monomers are confined between two repulsive walls at the bulk melt density *ρ* = 0.85*σ*^{−3}. Starting from a confined CG melt of 1000 chains of 80 soft spheres^{27,28} at rather high resolution such that each sphere corresponds to 25 monomers, it takes only 12.8 h CPU time on a single processor in Intel 3.6 GHz PC to prepare an initial configuration based on the bead-spring model (10 h for equilibrating the confined CG melt using a MC simulation and 2.8 h for reinserting monomers into soft spheres using a MD simulation). By gradually switching on the excluded volume interactions between two monomers, overlapping monomers are pushed away slowly in a warm-up procedure. Finally, the confined polymer melt is relaxed with full standard potentials. This takes about 182 h CPU time using 48 cores (2.7 GHz) on Dual Intel Xeon Platinum 8168 (155 h for the warm-up procedure and 27 h for the relaxation procedure). Similarly, as found in the previous studies,^{43,46} the required MD time for equilibrating confined polymer melts based on the bead-spring model is only about *t* = 12 900*τ* ≈ 5.69*τ*_{e}.

Following the same strategy, one can easily equilibrate highly entangled polymer melts confined between two walls at distances ranging from thick films to thin films (in which the distance between two walls is smaller than the radius of gyration of chains) within easily manageable computing time. Our work opens ample possibilities to study static and dynamic properties of highly entangled polymer chains in large polymer films, including, e.g., entanglement distributions. Varying the interaction potential between walls and monomers, or even replacing the wall potential by other potentials, only requires local short relaxation runs starting from a fully equilibrated polymer melt confined between two repulsive walls. Switching to our recently developed coarse-grained model for studying polymer melts under cooling,^{48,49} both fully equilibrated confined and free-standing films at the temperature *T* = 1.0*ε*/*k*_{B} and pressure *P* = 0.0*ε*/*σ*^{3} are also obtained in this work. This provides a direct route to further investigate the relation between the glass transition temperature and the thickness of films of highly entangled polymer chains at the zero pressure. Beyond that, it is also interesting to analyze the rheological properties and local morphology of deformed films by stretching or shearing.

## ACKNOWLEDGMENTS

H.-P.H. thanks T. Ohkuma and T. Stuehn for helpful discussions. The authors also thank K. Ch. Daoulas for carefully reading our paper. This work was supported by the European Research Council under the European Union’s Seventh Framework Programme (Grant No. FP7/2007-2013)/ERC Grant Agreement No. 340906-MOLPROCOMP. The authors also gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JUWELS at the Jülich Supercomputing Centre (JSC).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.