A coarse-grained polymer model for studying the glass transition
aa r X i v : . [ c ond - m a t . s o f t ] S e p A coarse-grained polymer model for studying the glass transition
Hsiao-Ping Hsu ∗ and Kurt Kremer † Max-Planck-Institut für Polymerforschung, Ackermannweg 10, 55128, Mainz, Germany
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 purelyrepulsive interactions an attractive potential between non-bonded monomers is added, such that the pressureof polymer melts is tuned to zero. Additionally, the commonly used bond bending potential controlling thechain stiffness is replaced by a new bond bending potential. For this model, we show that the Kuhn length andthe internal distances along the chains in the melt only very weakly depend on temperature, just as for typicalexperimental systems. The glass transition is observed by the temperature dependency of the melt density andthe characteristic non-Arrhenius slowing down of the chain mobility. The new model is set to allow for a fastswitch between models, for which a wealth of data already exists.
Polymer materials are omnipresent in our daily life with ap-plications in medicine, technology as well as as ‘simple’ com-modities to name a few. Very often these materials are in theglassy state [1]. In the liquid more rubbery state the viscositydramatically increases close to the glass transition tempera-ture T g in a non-Arrhenius way [2–5]. This slowing down ofthe chain mobility is of both high scientific and technologicalinterest. Experimentally, T g of polymers can be determinedas such by observing the change in the heat capacity of poly-mers using differential scanning calorimetry (DSC) [6], or bymeasuring the thermal expansion coefficient using thermo me-chanical analysis (TMA) [7]. However, the nature of the glasstransition is still not fully understood [8–13]. It is the purposeof this communication to present a most simple, efficient beadspring model, which allows to study these effects and whichcan make contact to the huge body of simulation data avail-able in the literature.Computer simulations play an important role in investigat-ing the structure and molecular motion (viscosity) of poly-meric systems under a variety of different conditions. Forstudying glassy polymers, both atomistic and coarse-grainedmodels are widely used in the literature [10, 11]. The struc-ture and thermal behavior of fluid mixtures can also be an-alyzed by tuning relative resolution in a recently developedhybrid model combing the fine-grained and coarse-grainedmodels [14]. Our aim is to eventually study generic prop-erties of large and highly entangled polymer melts in bulk, inconfinement and with free surfaces as a function of tempera-ture within accessible computing times. For this we adopt ahighly efficient coarse-grained model [15]. Usually in thesemodels the excluded volume interaction is taken care of bya purely repulsive Lennard-Jones (LJ) potential, the Weeks-Chandler-Andersen (WCA) potential [15], which prevents thestudy of surfaces [16, 17] and displays a rather high pressure( P ≈ . ε / σ , T = . ε / k B , density ρ = . σ − in stan-dard 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 monomersis often doubled from r cut = / σ to r c = r cut , resulting ∗ [email protected] † [email protected] P = . ε / σ [18–23]. The two main shortages of this set-ting are: (1) There is a small discontinuity in the force at thecut-off making microcanonical runs impossible and (2) thepressure is still not very close to zero. Furthermore, chainstiffness usually is taken into account by a bond bending po-tential [24–26], which tends to stretch the chains out with de-creasing temperatures [27]. As will be shown below, this leadsto rather artificial chain conformations upon cooling, while inexperiment chain conformations only very weakly depend ontemperature [28, 29]. Our new coarse-grained model is set toovercome these shortages.Our starting point is the standard bead spring model(BSM) [15] with a weak bending elasticity [24] (the bendingstrength k θ = . ε ) for which a huge body of data already ex-ists (see e.g. [25, 30–34]). While focusing on k θ = . ε , ourapproach easily applies to other bending constants as well.At the standard melt density of 0 . σ − ( σ being the unit oflength) the weak bending elasticity combined with the chainpacking result in an entanglement length of only N e = N e =
28 is small enough to allow for extremelyefficient simulations of highly entangled, huge polymeric sys-tems, while at the same time the subchain of length N e is al-ready well described by a Gaussian chain. The purpose ofthis communication is to replace/extend the WCA excludedvolume interaction potential to arrive at a pressure of P = . ε / σ , which allows to study free surfaces in interactionwith gases, liquids, and particles for example, and to replacethe standard bending potential U ( old ) BEND ( θ ) = k θ ( − cos θ ) bya new modified U BEND ( θ ) , which should lead to the typicalvery weak temperature dependence of chain conformationsin melts. The close resemblance to the standard semiflexiblebead spring model will allow to switch “on the fly” betweenthe models and to make use of the already broadly availabledata.In a first step we add an attractive well to the WCA ex-cluded volume in order to reduce the pressure in the systemfrom P = . ε / σ to P = . ε / σ . For this we add U ATT ( r ) (see Figure 1a), U ATT ( r ) = α (cid:20) cos ( π (cid:16) rr cut (cid:17) ) (cid:21) , r cut ≤ r < r ac , otherwise , (1) (a) -2 0 2 4 6 0.7 1.1225 1.5874r cut r ca P o t en t i a l [ ε ] r [ σ ] U WCA (r) U
ATT (r) (b) -4-2 0 2 4 0 π /4 π /2 3 π /4 πθ c P o t en t i a l [ ε ] θ U BEND(old) ( θ )U BEND ( θ ) FIG. 1. (a) Non-bonded and short-range repulsive potential U WCA ( r ) and attractive potential U ATT ( r ) with α = . ε {Eq. (1)} plotted asa function of distance r . (b) Standard and new bond bending potentials, U ( old ) BEND ( θ ) with k θ = . ε and U BEND ( θ ) with a θ = . ε , b θ = . θ . In (a)(b), the cut-off values are pointed by arrows. between all non-bonded monomers. U ATT ( r ) is set to not alterthe local bead packing. It is chosen to have zero force at thecut-off as well as at the contact point between the two partsof the potential at r c = / σ , which is needed in the casemicrocanonical simulations are performed. As illustrated inFigure 2a, adding this term to the standard model equilibratesand reduces the pressure to zero in less than 5 τ ( τ being thestandard LJ unit of time). This time corresponds to a small,local bead displacement of about 1 σ , for which the character-istic time is [32] τ ≈ . τ . Furthermore, since the num-ber of particles Z in the interaction range r ac = . σ is ≈
15 instead of ≈
45 at r c = . σ ( P = . ε / σ ) or ≈
60 at r c = . σ ( P = . ε / σ ) using the standard LJ potential, thepresent model is computationally significantly more efficient.In the next step we replace the standard bond bending poten-tial U ( old ) BEND ( θ ) = k θ ( − cos θ ) which would lead to a rod-likechain in the ground state at T = . ε / k B by a new bendingpotential U BEND ( θ ) with the goal to (1) match the chain con-formations at T = . ε / k B and (2) to approximately preservethem upon cooling. Thus it should satisfy the condition thatthe mean square end-to-end distance of chains, h R e i , does not(preferably) or only very weakly depend on the temperature T . The new bond bending potential U BEND ( θ ) (see Figure 1b)is chosen as U BEND ( θ ) = − a θ sin ( b θ θ ) , < θ < θ c (2)with the bond angle θ defined by θ = cos − (cid:18) ~ b j · ~ b j + | ~ b j || ~ b j + | (cid:19) where ~ b j = ~ r j − ~ r j − is the bond vector between monomers j and ( j − ) along the chain. The fitting parameters a θ and b θ ,and the cut-off θ c = π / b θ where the force | ~ F ( θ = θ c ) | = h R ( s ) i for all chemical distance s between twomonomers along the same chain follow the same curve as ob-tained from the model using U ( old ) BEND ( θ ) with k θ = . ε . Com-paring to the reference data for a polymer melt of n c = N =
50 shown in Figure 2b, we find that a θ = . ε , b θ = .
5. leads to an almost perfect match of the two systems. Our dataare also in perfect agreement with the theoretical predictiondescribed by a freely rotating chain (FRC) model [32, 35].Compared to the original model, the profiles of the pairdistribution function g ( r ) of all, inter, and intra pairs ofmonomers for polymer melts show that the two potentials U BEND ( θ ) and U ATT ( r ) only have very small effects on thelocal packing of monomers (Figure 2c). Results of the collec-tive structure factor S ( q ) also show that using the new model,the occurrence of the first peak remains at q = q ∗ ≈ . σ − indicating the same mean distance between monomers in thefirst neighbor shell of the polymer melt. The peak itself isslightly higher, indicating a slightly more structured local en-vironment, in agreement with the observed weakly enhancedbead friction.We now turn to the temperature dependency and comparemelts of the new model to the standard semiflexible polymermodel. For that we perform molecular dynamics (MD) sim-ulations (Hoover Barostat with Langevin thermostat [36, 37]implemented in ESPResSo++ [38]) at constant temperature T by a stepwise cooling [20], and constant pressure P = . ε / σ ( P = . ε / σ for the old model), i.e. in the isothermal-isobaric ensemble (NPT), for two polymer melts of n c = N =
50, and n c = N = ∆ T = . ε / k B with a relaxationtime between each step of ∆ t = τ resulting in a cool-ing rate of Γ = ∆ T / ∆ t = . × − ε / ( k B τ ) . ∆ t correspondsto ≈ . τ R , N = ≈ . τ R , N = ( τ R , N being the Rouse timeof the chains at T = . ε / k B for the old model). Results ofthe mean square internal distances h R ( s ) i and the bond an-gle probability distribution, P ( θ ) , are shown in Figures 3, 4.First let us focus on the standard weakly semiflexible model.As temperature decreases the chains stretch out as displayedin Figure 3a for N =
50. While for N =
50 the cooling rateis slow enough to allow for equilibration over a wide temper-ature range, for longer chains ( N = s ≤ h R ( s ) i / s . (a) c = 1000, N = 2000 α =0.5145 ε P [ ε / σ ] t [ τ ] U (old)BEND ( θ )U (old)BEND ( θ )+U ATT (r)U
BEND ( θ )+U ATT (r) (b) c = 2000 < R ( s ) > / s [ σ ] s U (old)BEND ( θ )U BEND ( θ )+U ATT (r) (c) c = 2000, N = 50 g (r) r [ σ ]U BEND ( θ )+U ATT (r)U (old)BEND ( θ ) (d) c = 2000 S ( q ) q [ σ -1 ]U BEND ( θ )+U ATT (r)U (old)BEND ( θ ) FIG. 2. (a) Pressure P plotted versus the relaxation time t . (b) Rescaled mean square internal distance, h R ( s ) i / s , plotted versus the chemicaldistance s between two monomers along the same chain. (c) Radial distribution function g ( r ) plotted as a function of r for all, inter, and intrapairs of monomers, as indicated. (d) Collective structure factor S ( q ) plotted versus the wave factor q . Polymer melts at T = . ε / k B describedby the standard BSM with additional potentials U ( old ) BEND ( θ ) , U ATT ( r ) , and U BEND ( θ ) are shown, as indicated. For long chain simulations, it will not be possible to avoidthis artefact. Also the strong increase of h R ( s ) i of the stan-dard semiflexible polymer model, is an artefact of the modelwhen compared to experiments. This increase in chain stiff-ness is related to the shift of the probability distribution P ( θ ) towards smaller angles as revealed in Figure 4a and whichdirectly connects to the shape of the standard bending poten-tial [22]. In contrast, the new excluded volume and bendingpotential not only leads to a conformational very close matchwith the old one at T = . ε / k B , but it also avoids a signifi-cant temperature shift. Figure 4b demonstrates for N =
50 that h R ( s ) i / s becomes independent of T within the error bars.As a consequence we also do not observe the maximum in h R ( s ) i for N =
500 as a function of temperature (Figure 3d).These observations fit to the T dependence of the distribution P ( θ ) , Figure 4b, which only becomes somewhat sharper butdoes not reveal any shift of the maximum.Finally we report some preliminary results for our newmodel in the glass transition region. As we are not interestedhere in details of the transition itself, we focus on N = n c = Γ = . × − ε / ( k B τ )) ,which, however, allows for a full relaxation of the system upto the region very close to T g , the observed glass transition temperature. T g can be determined from the change of den-sity ρ or volume V as a function of temperature [20]. Theintersection of linear extrapolation of ln V ( T ) between theliquid branch (ln V liquid = a liquid + α liquid T ) and glass branch(ln V glass = a glass + α glass T ) gives a good estimate of T g . Here α liquid and α glass are thermal expansion coefficients for poly-mer melts in the liquid state and the glass state, respectively.Results of ln V plotted versus T are shown in Figure 5a. Theglass transition occurs around T g = . ε / k B . To investigatethe mobility of chains at T > T g , we perform additional NVTMD simulations with a weak coupling Langevin thermostatfor polymer melts at k B T / ε = .
0, 0 .
95, 0 .
90, 0 .
85, 0 . .
75, 0 .
70, and 0 .
65. The initial configuration and volumeof the polymer melt at each temperature T are taken from thelast configuration of the NPT run in the cooling process. Ac-cording to the Rouse model [39], the mean square displace-ment (MSD) of monomers, g ( t ) , is expressed in terms ofthe Rouse rate W = k B T / ( πζσ ) as g ( t ) = σ ( Wt ) / .Here ζ ( ∝ D − ∝ η ) being the monomeric friction coefficientis related to the self-diffusion coefficient D = k B T / ( N ζ ) andthe viscosity η using the Stokes-Einstein relation. Results of g ( t ) taking from the average MSD of inner 12 monomers areshown in Figure 5b. We also include the data at T = . ε / k B (a) (old)BEND ( θ )n c = 2000 < R ( s ) > / s [ σ ] s k B T/ ε =0.500.600.700.800.901.00 (b) BEND ( θ )n c = 2000 < R ( s ) > / s [ σ ] sk B T/ ε = 0.500.600.700.800.901.00 (c) c = 1000, N = 500U (old)BEND ( θ ) < R ( s ) > / s [ σ ] s k B T/ ε =0.500.600.700.800.901.00 (d) c = 1000, N = 500U BEND ( θ ) < R ( s ) > / s [ σ ] s k B T/ ε =0.500.600.700.800.901.00 FIG. 3. Rescaled mean square internal distance, h R ( s ) i / s , plotted as a function of chemical distance s for polymer melts described by thestandard BSM with the original and new bond bending potentials, U ( old ) BEND ( θ ) (a)(c), and U BEND ( θ ) (b)(d), respectively, at P = . ε / σ . Thetheoretical prediction for FRC with [32] h cos θ i = . n c = N = π /4 π /2 3 π /4N = 50U (old)BEND ( θ )n c = 2000 P ( θ ) θ k B T/ ε = 0.500.600.700.800.901.00 (b) π /4 π /2 3 π /4N = 50U BEND ( θ )n c = 2000 P ( θ ) θ k B T/ ε = 0.500.600.700.800.901.00 FIG. 4. Probability distribution of bond angle θ for polymer melts described by the standard BSM with the original and new bond bendingpotential U ( old ) BEND ( θ ) (a), and U BEND ( θ ) (b)), respectively. for the old model for comparison. The Rouse rate W de-pending on the temperature is determined by the best fit ofa straight line with slope 1 / T = . ε / k B , the Rouse rate for the old model ( W = . τ − ) is faster than the new model ( W = . τ − ).From the well-known Vogel-Fulcher-Tammann (VFT) equa-tion [40–42], log η = A + BT − T , where A , B , and T areconstants and T is the absolute temperature, Angell [8, 9] (a) B T g / ε = 0.64 liquidglassn c = 2000, N = 50 l n ( V / σ ) k B T / ε (b) -2 -1 -1
1 10 10 k B T/ ε t t (old) n c = 2000, N = 50 g ( t ) [ σ ] t [ τ ]1.001.000.900.800.750.700.65 (c) c = 2000, N = 50 l og ( / W ) T g / TVFT EquationFragility parameter m FIG. 5. (a) Logarithm of volume of the system, ln V / σ , plotted versus temperature k B T / ε . The two linear lines give the best fit of our dataalong the liquid branch ( a liquid = . α liquid = . k B / ε ) and the glass branch ( a glass = . α glass = . k B / ε ). (b) Time evolution ofmean square displacement of inner monomers, g ( t ) at various chosen temperatures T . The predicted scaling laws by Rouse model are shownby straight lines. (c) Common logarithm of the inverse of the rate constant W estimated from (b), log ( / W ) , plotted versus T g / T . Data for k B T / ε > . m ( T ) = . ( T g / T ) . + . ( / W ) = A + B / ( T − T ) with A = . B = . ε / k B , and T = . ε / k B is shown by a solid curve. has proposed that the fragility parameter m , defined by [43]: m = d ( log η ) / d ( T g / T ) | T = T g . Thus, plotting log ( / W ) versus T g / T in Figure 5c, we obtain the characteristic behav-ior of a polymer approaching the glass transition.In summary, based on the standard BSM, we have in-troduced a new non-bonded short range attractive potential U ATT ( r ) and bond bending potential U BEND ( θ ) for studyingpolymer melts subject to cooling. The functional form ofthese two new interaction potentials also is directly applicableto other standard BSM models with different stiffness [25] justby adjusting the coefficients. By keeping α = . ε , whichresults in a density of 0 . σ − for all longest ( N = a θ = . ε for 0 ≤ k θ / ε ≤ .
0, and b θ = .
32, 1 .
40, 1 .
50, and 1 .
70 for k θ / ε = .
5, 1 . .
5, and 2 .
0, respectively. The new coarse-grained model cap-tures the major features of glass-forming polymers, and pre-serves the Kuhn length as well as internal distances and can also be used to study systems with free surfaces. By con-struction it can directly take advantage of available simulationdata of standard BSM models at T = . ε / k B and can be ap-plied to available large deformed polymer melts [33, 34] andfor understanding the viscoelastic behavior of these polymericsystems.Acknowledgement: We are grateful to B. Dünweg for a crit-ical reading of the manuscript. This work has been supportedby European Research Council under the European Union’sSeventh Framework Programme (FP7/2007-2013)/ERC GrantAgreement No. 340906-MOLPROCOMP. We also gratefullyacknowledge the computing time granted by the John vonNeumann Institute for Computing (NIC) and provided onthe supercomputer JUROPA at Jülich Supercomputing Cen-tre (JSC), and the Max Planck Computing and Data Facility(MPCDF). [1] J. H. Anantrao, J. C. Motichand, and B. S. Narhari, Int. J. Adv.Res. , 671 (2017). [2] R. H. Boyd, R. H. Gee, J. Han, and Y. Jin, J. Chem. Phys. ,
788 (1994).[3] C. A. Angell, Science , 1924 (1995).[4] M. Paluch, J. Gapinski, A. Patkowski, and E. W. Fischer, J.Chem. Phys. , 8048 (2001).[5] L. Berthier and G. Biroli, Rev. Mod. Phys. , 587 (2011).[6] V. B. F. Mathot, Calorimetry and thermal analysis of polymers (Hanser, Munich, 1994).[7] R. Bird, C. Curtis, and R. Armstrong,
Dynamics of PolymerFluids, 2nd ed. (Wiley, New York, 1987).[8] C. A. Angell, J. Phys. Chem. Solids , 863 (1988).[9] C. A. Angell, J. Non-Cryst. Solids , 13 (1991).[10] K. Binder and W. Kob, Glassy Materials and Disordered Solids (World Scientific, Singapore, 2005).[11] J.-L. Barrat, J. Baschnagel, and A. Lyulin, Soft Matter , 3430(2010).[12] F. H. Stillinger and P. G. Debenedetti, Annu. Rev. Condens.Matter Phys. , 263 (2013).[13] M. D. Ediger and J. A. Forrest, Macromolecules , 471 (2014).[14] A. Chaimovich, C. Peter, and K. Kremer, J. Chem. Phys. ,243107 (2015).[15] K. Kremer and G. S. Grest, J. Chem. Phys. , 5057 (1990).[16] B. Dünweg, G. S. Grest, and K. Kremer, Conf. Proc. of theIMA Workshop (Minneapolis MN, May 1996) (Berlin:Springer,1997).[17] A. Kopf, B. Dünweg, and W. Paul, J. Chem. Phys. , 6945(1997).[18] C. Bennemann, W. Paul, K. Binder, and B. Dünweg, Phys. Rev.E , 843 (1998).[19] K. Binder, Comp. Phys. Commun. , 168 (1999).[20] J. Buchholz, W. Paul, F. Varnik, and K. Binder, J. Chem. Phys. , 7364 (2002).[21] K. Binder, J. Baschnagel, and W. Paul, Prog. Polym. Sci. ,115 (2003).[22] B. Schnell, H. Meyer, C. Fond, J. P. Wittmer, andJ. Baschnagel, Eur. Phys. J. E , 97 (2011). [23] S. Frey, F. Weysser, H. Meyer, J. Farago, M. Fuchs, andJ. Baschnagel, Eur. Phys. J. E. , 11 (2015).[24] R. Everaers, S. K. Sukumaran, G. S. Grest, C. Svaneborg,A. Sivasubramanian, and K. Kremer, Science , 823 (2004).[25] C. Svaneborg and R. Everaers, arXiv:1808.03503 (2018).[26] C. Svaneborg, H. A. Karimi-Varzaneh, N. Hojdis, F. Fleck, andR. Everaers, arXiv:1808.03509 (2018).[27] G. S. Grest, J. Chem. Phys. , 141101 (2016).[28] M. Wind, R. Graf, A. Heuer, and H. W. Spiess, Phys. Rev. Lett. , 155702 (2003).[29] L. J. Fetters, D. J. Lohse, and R. H. Colby, in Physical Proper-ties of Polymers Handbook, 2nd , edited by J. E. Mark (Springer,New York, 2007) Chap. 25, pp. 447–454.[30] G. Zhang, L. A. Moreira, T. Stuehn, K. C. Daoulas, and K. Kre-mer, ACS Macro Lett. , 198 (2014).[31] L. A. Moreira, G. Zhang, F. Müller, T. Stuehn, and K. Kremer,Macromol. Theor. Simul. , 419 (2015).[32] H.-P. Hsu and K. Kremer, J. Chem. Phys. , 154907 (2016).[33] H.-P. Hsu and K. Kremer, ACS Macro Lett. , 107 (2018).[34] H.-P. Hsu and K. Kremer, Phys. Rev. Lett. , 167801 (2018).[35] M. Rubinstein and R. H. Colby, Polymer Physics (Oxford Uni-versity Press, Oxford, 2003).[36] G. J. Martyna, D. J. Tobias, and M. L. Klein, J. Chem. Phys. , 4177 (1994).[37] D. Quigley and M. I. J. Probert, J. Chem. Phys. , 11432(2004).[38] J. D. Halverson, T. Brandes, O. Lenz, A. Arnold, S. Bevc,V. Starchenko, K. Kremer, T. Stuehn, and D. Reith, Comput.Phys. Commun. , 1129 (2013).[39] P. R. Rouse, J. Chem. Phys. , 1272 (1953).[40] H. Vogel, Phys. Z , 645 (1921).[41] G. S. Fulcher, J. Am. Ceram. Soc. , 339 (1925).[42] G. Tammann, W. Hesse, and Z. Anorg, Allgem. Chem. ,245 (1926).[43] M. L. F. Nascimento and C. Aparicio, J. Phys. Chem. Solids68