A coarse-grained polymer model for studying the glass transition

TL;DR

This paper introduces a new coarse-grained polymer model that effectively simulates the glass transition by capturing key physical behaviors like density changes and chain mobility slowdown, with minimal temperature dependence of chain structure.

Contribution

A novel coarse-grained polymer model with a new bond bending potential that accurately reproduces glass transition phenomena and allows easy switching between existing models.

Findings

Kuhn length and internal chain distances are weakly temperature-dependent.

Glass transition characterized by density changes and non-Arrhenius mobility slowdown.

Model enables fast switching between different polymer simulation models.

Abstract

To study the cooling behavior and the glass transition of polymer melts in bulk and with free surfaces a coarse-grained weakly semi-flexible polymer model is developed. Based on a standard bead spring model with purely repulsive interactions an attractive potential between non-bonded monomers is added, such that the pressure of polymer melts is tuned to zero. Additionally, the commonly used bond bending potential [Everaers et al., Science 303, 823 (2004)] controlling the chain stiffness is replaced by a new bond bending potential. For this model, we show that the Kuhn length and the internal distances along the chains in the melt only very weakly depend on temperature, just as for typical experimental systems. The glass transition is observed by the temperature dependency of the melt density and the characteristic non-Arrhenius slowing down of the chain mobility. The model is set to allow for a fast switch between models, for which a wealth of data already exists.

Equations

$$\displaystyle U_{\rm ATT}(r)=\left\{\begin{array}[]{ll}\alpha\left[\cos(\pi\left(\frac{r}{r_{\rm cut}}\right)^{2})\right],&{r_{\rm cut}}\leq r<r^{a}_{c}\\ &\\ 0,&{\rm otherwise}\end{array}\right.\,,$$

$$U_{\mathrm{BEND}}(\theta )=-a_{\theta }{\sin }^2(b_{\theta }\theta ),0<\theta <{\theta }_c$$

Full text

Computer simulations play an important role in investigating the structure and molecular motion (viscosity) of polymeric systems under a variety of different conditions. For studying glassy polymers, both atomistic and coarse-grained models are widely used in the literature Binder and Kob (2005); Barrat et al. (2010). The structure and thermal behavior of fluid mixtures can also be analyzed by tuning relative resolution in a recently developed hybrid model combing the fine-grained and coarse-grained models Chaimovich et al. (2015). Our aim is to eventually study generic properties of large and highly entangled polymer melts in bulk, in confinement and with free surfaces as a function of temperature within accessible computing times. For this we adopt a highly efficient coarse-grained model Kremer and Grest (1990). Usually in these models the excluded volume interaction is taken care of by a purely repulsive Lennard-Jones (LJ) potential, the Weeks-Chandler-Andersen (WCA) potential Kremer and Grest (1990), which prevents the study of surfaces Dünweg et al. (1997); Kopf et al. (1997) and displays a rather high pressure ($P\approx 5.0\epsilon/\sigma^{3}$, $T=1.0\epsilon/k_{B}$, density $\rho=0.85\sigma^{-3}$ in standard Lennard-Jones (LJ) units of energy and length, and $k_{B}$ being the Boltzmann factor). To reduce the pressure the cut-off of the WCA potential for non-bonded pairs of monomers is often doubled from $r_{\rm cut}=2^{1/6}\sigma$ to $r_{c}=2r_{\rm cut}$, resulting $P=1.0\epsilon/\sigma^{3}$ Bennemann et al. (1998); Binder (1999); Buchholz et al. (2002); Binder et al. (2003); Schnell et al. (2011); Frey et al. (2015). The two main shortages of this setting are: (1) There is a small discontinuity in the force at the cut-off making microcanonical runs impossible and (2) the pressure is still not very close to zero. Furthermore, chain stiffness usually is taken into account by a bond bending potential Everaers et al. (2004); Svaneborg and Everaers (2018); Svaneborg et al. (2018), which tends to stretch the chains out with decreasing temperatures Grest (2016). As will be shown below, this leads to rather artificial chain conformations upon cooling, while in experiment chain conformations only very weakly depend on temperature Wind et al. (2003); Fetters et al. (2007). Our new coarse-grained model is set to overcome these shortages.

Our starting point is the standard bead spring model (BSM) Kremer and Grest (1990) with a weak bending elasticity Everaers et al. (2004) (the bending strength $k_{\theta}=1.5\epsilon$) for which a huge body of data already exists (see e.g. Zhang et al. (2014); Moreira et al. (2015); Hsu and Kremer (2016, 2018a, 2018b); Svaneborg and Everaers (2018)). While focusing on $k_{\theta}=1.5\epsilon$, our approach easily applies to other bending constants as well. At the standard melt density of $0.85\sigma^{-3}$ ($\sigma$ being the unit of length) the weak bending elasticity combined with the chain packing result in an entanglement length of only $N_{e}=28$ monomers. $N_{e}=28$ is small enough to allow for extremely efficient simulations of highly entangled, huge polymeric systems, while at the same time the subchain of length $N_{e}$ is already well described by a Gaussian chain. The purpose of this communication is to replace/extend the WCA excluded volume interaction potential to arrive at a pressure of $P=0.0\epsilon/\sigma^{3}$, which allows to study free surfaces in interaction with gases, liquids, and particles for example, and to replace the standard bending potential $U_{\rm BEND}^{\rm(old)}(\theta)=k_{\theta}(1-\cos\theta)$ by a new modified $U_{\rm BEND}(\theta)$, which should lead to the typical very weak temperature dependence of chain conformations in melts. The close resemblance to the standard semiflexible bead spring model will allow to switch “on the fly” between the models and to make use of the already broadly available data.

In a first step we add an attractive well to the WCA excluded volume in order to reduce the pressure in the system from $P=5.0\epsilon/\sigma^{3}$ to $P=0.0\epsilon/\sigma^{3}$. For this we add $U_{\rm ATT}(r)$ (see Figure 1a),

[TABLE]

between all non-bonded monomers. $U_{\rm ATT}(r)$ is set to not alter the local bead packing. It is chosen to have zero force at the cut-off as well as at the contact point between the two parts of the potential at $r_{c}=2^{1/6}\sigma$, which is needed in the case microcanonical simulations are performed. As illustrated in Figure 2a, adding this term to the standard model equilibrates and reduces the pressure to zero in less than $5\tau$ ($\tau$ being the standard LJ unit of time). This time corresponds to a small, local bead displacement of about $1\sigma$, for which the characteristic time is Hsu and Kremer (2016) $\tau_{0}\approx 2.89\tau$. Furthermore, since the number of particles $Z$ in the interaction range $r_{c}^{a}=1.5874\sigma$ is $\approx 15$ instead of $\approx 45$ at $r_{c}=2.25\sigma$ $(P=1.0\epsilon/\sigma^{3})$ or $\approx 60$ at $r_{c}=2.5\sigma$ $(P=0.0\epsilon/\sigma^{3})$ using the standard LJ potential, the present model is computationally significantly more efficient. In the next step we replace the standard bond bending potential $U^{\rm(old)}_{\rm BEND}(\theta)=k_{\theta}(1-\cos\theta)$ which would lead to a rod-like chain in the ground state at $T=0.0\epsilon/k_{B}$ by a new bending potential $U_{\rm BEND}(\theta)$ with the goal to (1) match the chain conformations at $T=1.0\epsilon/k_{B}$ and (2) to approximately preserve them upon cooling. Thus it should satisfy the condition that the mean square end-to-end distance of chains, $\langle R_{e}^{2}\rangle$, does not (preferably) or only very weakly depend on the temperature $T$. The new bond bending potential $U_{\rm BEND}(\theta)$ (see Figure 1b) is chosen as

[TABLE]

with the bond angle $\theta$ defined by $\theta=\cos^{-1}\left(\frac{\vec{b}_{j}\cdot\vec{b}_{j+1}}{\mid\vec{b}_{j}\mid\mid\vec{b}_{j+1}\mid}\right)$ where $\vec{b}_{j}=\vec{r}_{j}-\vec{r}_{j-1}$ is the bond vector between monomers $j$ and $(j-1)$ along the chain. The fitting parameters $a_{\theta}$ and $b_{\theta}$, and the cut-off $\theta_{c}=\pi/{b_{\theta}}$ where the force $\mid\vec{F}(\theta=\theta_{c})\mid=0$ are adjusted such that the estimates of the mean square internal distance $\langle R^{2}(s)\rangle$ for all chemical distance $s$ between two monomers along the same chain follow the same curve as obtained from the model using $U^{\rm(old)}_{\rm BEND}(\theta)$ with $k_{\theta}=1.5\epsilon$. Comparing to the reference data for a polymer melt of $n_{c}=2000$, $N=50$ shown in Figure 2b, we find that $a_{\theta}=4.5\epsilon$, $b_{\theta}=1.5$. leads to an almost perfect match of the two systems. Our data are also in perfect agreement with the theoretical prediction described by a freely rotating chain (FRC) model Rubinstein and Colby (2003); Hsu and Kremer (2016).

Compared to the original model, the profiles of the pair distribution function $g(r)$ of all, inter, and intra pairs of monomers for polymer melts show that the two potentials $U_{\rm BEND}(\theta)$ and $U_{\rm ATT}(r)$ only have very small effects on the local packing of monomers (Figure 2c). Results of the collective structure factor $S(q)$ also show that using the new model, the occurrence of the first peak remains at $q=q*\approx 6.9\sigma^{-1}$ indicating the same mean distance between monomers in the first neighbor shell of the polymer melt. The peak itself is slightly higher, indicating a slightly more structured local environment, in agreement with the observed weakly enhanced bead friction.

We now turn to the temperature dependency and compare melts of the new model to the standard semiflexible polymer model. For that we perform molecular dynamics (MD) simulations (Hoover Barostat with Langevin thermostat Martyna et al. (1994); Quigley and Probert (2004) implemented in ESPResSo++ Halverson et al. (2013)) at constant temperature $T$ by a stepwise cooling Buchholz et al. (2002), and constant pressure $P=0.0\epsilon/\sigma^{3}$ ($P=5.0\epsilon/\sigma^{3}$ for the old model), i.e. in the isothermal-isobaric ensemble (NPT), for two polymer melts of $n_{c}=2000$, $N=50$, and $n_{c}=1000$, $N=500$, respectively. The temperature is reduced in steps of $\Delta T=0.05\epsilon/k_{B}$ with a relaxation time between each step of $\Delta t=60000\tau$ resulting in a cooling rate of $\Gamma=\Delta T/\Delta t=8.3\times 10^{-7}\epsilon/(k_{B}\tau)$. $\Delta t$ corresponds to $\approx 8.3\tau_{R,N=50}\approx 0.083\tau_{R,N=500}$ ($\tau_{R,N}$ being the Rouse time of the chains at $T=1.0\epsilon/k_{B}$ for the old model). Results of the mean square internal distances $\langle R^{2}(s)\rangle$ and the bond angle probability distribution, $P(\theta)$, are shown in Figures 3, 4. First let us focus on the standard weakly semiflexible model. As temperature decreases the chains stretch out as displayed in Figure 3a for $N=50$. While for $N=50$ the cooling rate is slow enough to allow for equilibration over a wide temperature range, for longer chains ($N=500$, Figure 3c) the system cannot equilibrate anymore even on short length scales ($s\leq 50$), leading to a characteristic maximum in $\langle R^{2}(s)\rangle/s$. For long chain simulations, it will not be possible to avoid this artefact. Also the strong increase of $\langle R^{2}(s)\rangle$ of the standard semiflexible polymer model, is an artefact of the model when compared to experiments. This increase in chain stiffness is related to the shift of the probability distribution $P(\theta)$ towards smaller angles as revealed in Figure 4a and which directly connects to the shape of the standard bending potential Schnell et al. (2011). In contrast, the new excluded volume and bending potential not only leads to a conformational very close match with the old one at $T=1.0\epsilon/k_{B}$, but it also avoids a significant temperature shift. Figure 4b demonstrates for $N=50$ that $\langle R^{2}(s)\rangle/s$ becomes independent of $T$ within the error bars. As a consequence we also do not observe the maximum in $\langle R^{2}(s)\rangle$ for $N=500$ as a function of temperature (Figure 3d). These observations fit to the $T$ dependence of the distribution $P(\theta)$, Figure 4b, which only becomes somewhat sharper but does not reveal any shift of the maximum.

Finally we report some preliminary results for our new model in the glass transition region. As we are not interested here in details of the transition itself, we focus on $N=50$ ($n_{c}=2000$) and one cooling rate ($\Gamma=8.3\times 10^{-7}\epsilon/(k_{B}\tau))$, which, however, allows for a full relaxation of the system up to the region very close to $T_{g}$, the observed glass transition temperature. $T_{g}$ can be determined from the change of density $\rho$ or volume $V$ as a function of temperature Buchholz et al. (2002). The intersection of linear extrapolation of $\ln V(T)$ between the liquid branch ($\ln V_{\rm liquid}=a_{\rm liquid}+\alpha_{\rm liquid}T$) and glass branch ($\ln V_{\rm glass}=a_{\rm glass}+\alpha_{\rm glass}T$) gives a good estimate of $T_{g}$. Here $\alpha_{\rm liquid}$ and $\alpha_{\rm glass}$ are thermal expansion coefficients for polymer melts in the liquid state and the glass state, respectively. Results of $\ln V$ plotted versus $T$ are shown in Figure 5a. The glass transition occurs around $T_{g}=0.64\epsilon/k_{B}$. To investigate the mobility of chains at $T>T_{g}$, we perform additional NVT MD simulations with a weak coupling Langevin thermostat for polymer melts at $k_{B}T/\epsilon=1.0$, $0.95$, $0.90$, $0.85$, $0.80$, $0.75$, $0.70$, and $0.65$. The initial configuration and volume of the polymer melt at each temperature $T$ are taken from the last configuration of the NPT run in the cooling process. According to the Rouse model Rouse (1953), the mean square displacement (MSD) of monomers, $g_{1}(t)$, is expressed in terms of the Rouse rate $W=12k_{B}T/(\pi\zeta\sigma^{2})$ as $g_{1}(t)=\sigma^{2}(Wt)^{1/2}$. Here $\zeta(\propto D^{-1}\propto\eta)$ being the monomeric friction coefficient is related to the self-diffusion coefficient $D=k_{B}T/(N\zeta)$ and the viscosity $\eta$ using the Stokes-Einstein relation. Results of $g_{1}(t)$ taking from the average MSD of inner $12$ monomers are shown in Figure 5b. We also include the data at $T=1.0\epsilon/k_{B}$ for the old model for comparison. The Rouse rate $W$ depending on the temperature is determined by the best fit of a straight line with slope $1/2$ going through our data on log-log scales. At $T=1.0\epsilon/k_{B}$, the Rouse rate for the old model ($W=0.20\tau^{-1}$) is faster than the new model ($W=0.09\tau^{-1}$). From the well-known Vogel-Fulcher-Tammann (VFT) equation Vogel (1921); Fulcher (1925); Tammann et al. (1926), $\log_{10}\eta=A+\frac{B}{T-T_{0}}$, where $A$, $B$, and $T_{0}$ are constants and $T$ is the absolute temperature, Angell Angell (1988, 1991) has proposed that the fragility parameter $m$, defined by Nascimento and Aparicio (2007): $m=d(\log_{10}\eta)/d(T_{g}/T)\mid_{T=T_{g}}$. Thus, plotting $\log_{10}(1/W)$ versus $T_{g}/T$ in Figure 5c, we obtain the characteristic behavior of a polymer approaching the glass transition.

In summary, based on the standard BSM, we have introduced a new non-bonded short range attractive potential $U_{\rm ATT}(r)$ and bond bending potential $U_{\rm BEND}(\theta)$ for studying polymer melts subject to cooling. The functional form of these two new interaction potentials also is directly applicable to other standard BSM models with different stiffness Svaneborg and Everaers (2018) just by adjusting the coefficients. By keeping $\alpha=0.5145\epsilon$, which results in a density of $0.85\sigma^{-3}$ for all longest ($N=2000$) systems within the error bars, we get $a_{\theta}=4.5\epsilon$ for $0\leq k_{\theta}/\epsilon\leq 2.0$, and $b_{\theta}=1.32$, $1.40$, $1.50$, and $1.70$ for $k_{\theta}/\epsilon=0.5$, $1.0$, $1.5$, and $2.0$, respectively. The new coarse-grained model captures the major features of glass-forming polymers, and preserves the Kuhn length as well as internal distances and can also be used to study systems with free surfaces. By construction it can directly take advantage of available simulation data of standard BSM models at $T=1.0\epsilon/k_{B}$ and can be applied to available large deformed polymer melts Hsu and Kremer (2018a, b) and for understanding the viscoelastic behavior of these polymeric systems.

Acknowledgement: We are grateful to B. Dünweg for a critical reading of the manuscript. This work has been supported by European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 340906-MOLPROCOMP. We also gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JUROPA at Jülich Supercomputing Centre (JSC), and the Max Planck Computing and Data Facility (MPCDF).