Realistic Coarse-Grain Model of - ACS Publications - American

Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

pubs.acs.org/Macromolecules

Realistic Coarse-Grain Model of cis-1,4-Polybutadiene: From Chemistry to Rheology ́ y,† A. Dequidt,*,† M. Couty,‡ and P. Malfreyt*,† K. Kempfer,†,‡ J. Devem †

CNRS, SIGMA Clermont, Institut de Chimie de Clermont-Ferrand, Université Clermont Auvergne, F-63000 Clermont-Ferrand, France ‡ Manufacture Française des Pneumatiques Michelin, 23, Place des Carmes, 63040 Clermont-Ferrand, France

Macromolecules Downloaded from pubs.acs.org by UNIV OF NEW ENGLAND on 03/21/19. For personal use only.

S Supporting Information *

ABSTRACT: We present a practical coarse-grain (CG) force field for dissipative particle dynamics (DPD) simulations of cis1,4-polybutadiene (cPB). The parametrization of the model was realized using the recent Bayesian bottom-up scheme based on trajectory matching. The CG force field enables the correct description of cPB melts in a broad range of temperatures and pressures. We demonstrate its robustness by measuring the thermal expansion coefficient and isothermal compressibility consistent with experimental values. By reducing the computational cost by a factor 21 with respect to classical united-atom force fields while keeping accuracy, we can reach larger time and length scales suitable for the investigation of the rheological and the mechanical properties of such polymeric systems. In particular, the rubbery plateau and the entanglement mass for monodisperse cPB entangled melts are quantitatively predicted using this CG model.

■

INTRODUCTION Molecular dynamics (MD) is very expensive in terms of computational cost. Hence, classical MD simulations usually do not last more than a few tens of nanoseconds. This limitation becomes critical in the case of polymeric systems, which involve large time and length scales. In the tire industry, the molar mass of a single polymer chain can reach values up to 200 kg mol−1; i.e., one molecule contains thousands of repeating units and even more atoms. Worse still, the viscoelastic properties of polymer melts, which are key physical features for the design of new materials, result from a close cooperation between different polymer chains through entanglements (for un-cross-linked rubber) or cross-links (for vulcanized rubber), leading to even larger systems. The investigation of the mechanical and rheological properties of such demanding systems at the atomistic scale seems like an impossible dream. Increasing the number of central processing units (CPU) could partially resolve the box size problem (thanks to spatial decomposition), but under no circumstances will it address the main challenge that is to reach time scales of the order of magnitude of the micro- or millisecond within a reasonable timeline. Indeed, the speed of integration of the equation of motion for a given particle mainly relies on the CPU clock frequency that is constant and about 3−4 GHz since 2005 due to overheating and higher power consumption.1 For some years now,2−5 coarse-grain (CG) models have become very attractive to tackle macromolecular systems. The idea is to reduce the number of degrees of freedom by grouping multiple atoms together into one single coarse grain, also called a bead, to reduce computational cost. The prediction of target © XXXX American Chemical Society

characteristic properties from the mere knowledge of the underlying microstructure of a given material can be achieved by taking a bottom-up route,6,7 in which realistic CG force fields are derived from high-resolution atomistic simulations. This strategy is to put in contrast with top-down pathways that assume we already know some target properties (obtained for instance from experimental measurements), damaging the predictive aspect of the CG model. A great effort has been dedicated to building realistic CG force fields keeping enough resolution with respect to classical atomistic descriptions. Nevertheless, being quantitative remains challenging. For charged biomolecules, the Martini force field developed by Marrink et al.8 has been derived from both bottom-up (for the bonded interactions) and top-down (free energy calculations for the nonbonded interactions) approaches, resulting in a poor reproduction of the polymer conformations. In the field of polymers, the popular iterative Boltzmann inversion (IBI) method9,10 and its derivatives11 produce tabulated CG potentials faithful to the structure, as they are based on it. However, this technique gives no clue about dynamics, and the transferability in terms of pressure or temperature is not direct.7,12 On the other hand, the force matching method fits both the conservative forces and dissipative friction coefficients from the analysis of constrained atomistic simulations.13,14 A new method based on statistical Received: December 29, 2018 Revised: March 8, 2019

A

DOI: 10.1021/acs.macromol.8b02750 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules trajectory matching has been proposed recently,15,16 which can also fit the nonconservative interactions. In this paper, we will demonstrate the capabilities of this new method for polymers. Accordingly, we will propose a multitemperature, multipressure CG force field for dissipative particle dynamics simulations (DPD) of cis-1,4-polybutadiene (cPB). This material is of particular interest in the tire industry and has been studied extensively in the past, but no practical CG force field has been suggested yet. First, the parametrization procedure of the CG model will be detailed and the optimal parameters will be provided. The accuracy of our model will then be quantified against MD reference simulations. Finally, we will investigate the rheological and mechanical behavior based on 10 μs long DPD simulations of realistic un-cross-linked cPB melts. A summary will be given in the last section.

Coarse-Grain Model. Each monomer of the cPB chains was mapped into one coarse grain (CG) all along the different reference MD trajectories. The coordinates of the CG beads are thus defined as the centers of mass of the four UA defining the cPB repeating units. We recently demonstrated that this level of coarse-graining was the most suitable for cPB.24 This process effectively reduces the number of particles in the system by a factor 4 as described in Figure 1. More precisely, the resulting so-

■

FORCE FIELD OPTIMIZATION Atomistic Reference Trajectories. The coarse-grain model has been developed from the bottom up. It is not necessary to simulate highly entangled polymers at the atomistic scale to derive convenient CG potentials that are expected to be transferable to model longer chains. As an input, high-resolution molecular dynamics (MD) simulations of cis-1,4-polybutadiene (cPB) melts at different temperatures and pressures were performed using LAMMPS software.17 Cubic simulation boxes consisting of a periodic amorphous cell containing 90 cPB chains holding each 60 repeating units (their relaxation in MD is obtained after few nanoseconds) were prepared at a density close to its equilibrium value for given target thermodynamic conditions. The united-atom (UA) description proposed by Smith and Paul18 and Tsolou et al.19 was taken to model the polymer. The parameters of the atomistic force field are given in Table S1 of the Supporting Information. The cutoff for LennardJones interactions was set at 14 Å. The equation of motion was integrated using the standard velocity-Verlet algorithm with a time step of 2 fs. All the simulations were performed in the constant-NPT ensemble using Berendsen barostat20 and Langevin thermostat21,22 with 1 and 10 ps for the barostat and thermostat relaxation times, respectively. The bulk modulus parameter for Berendsen algorithm was set at 100 MPa. The equilibration protocol is divided into three stages. The initial configurations are first relaxed for 10 ns at high temperature 500 K in the constant-NVT ensemble. Then, a linear temperature ramp is applied for 2 ns to reach the target temperature. An additional 8 ns run in the constant-NPT ensemble at target pressure forms the last step of the equilibration procedure. Nine temperatures evenly spaced over the interval from 200 to 400 K (200, 225, ..., 400 K) at P = 0.1 MPa were selected to probe the transferability in temperature. Note that all the temperatures are here above the glass transition temperature, which was found at around 160 K in accordance with previous simulations10 and experimental values.23 In the same spirit, eight pressures ranging from −10 to −5, −1, −0.1, 0.1, 1, 5, and 10 MPa at constant temperature T = 300 K were picked to investigate the transferability in pressure. In total, 16 pairs of thermodynamic conditions (P, T) were considered. After equilibration, the 16 high-resolution atomistic reference trajectories ; were recorded every 25 time steps during a 25 ps production stage performed in the constant-NPT ensemble with independent control of the x, y, and z components of the external barostat stress tensor. The resulting trajectories ; comprise 501 configurations equally spaced by the same time interval Δt = 50 fs.

Figure 1. Coarse-graining of a cPB chain. One repeating unit (4 UA) of the polymer is replaced with one coarse particle.

called reference coarse-grained trajectories contain 5400 beads instead of 21600 UA before coarse-graining. The grain velocities vi and the total force f i acting on each of them were computed at each time step from the analysis of the coarse-grained reference trajectories using the Verlet integration scheme and Newton’s equation of motion as described in previous work.15 The other coarse degree of freedom is the volume evolution of the simulation boxes in the reference trajectories, which allows for the development of CG potentials able to reproduce the MD density.15,16 Coarse-grain DPD potentials were derived using the recent Bayesian bottom-up approach.15,16,24 The conservative nonbonded force was parametrized by a linear combination of 4 piecewise Bernstein polynomials.25 The force exerted by particle j on i is expressed by f jC→ i

l 3 o o o o o∑ xiPi(r )eij for r < rc =o m i=0 o o o o o o 0 for r ⩾ rc n

(1)

where xi are the weighting factors of the Bernstein basis polynomials Pi, r is the distance between i and j, eij is the unit vector in the direction ri − rj, and rc = 20 Å is the conservative cutoff. For completeness, full details about these Pi bases functions are given in the Supporting Information (see Figure S1 and eqs S1−S5). This nonbonded interaction is used between pairs of beads belonging to different molecules or to the same molecule but that are not directly connected through a mesoscopic bond. In addition, this interaction is half-screened in the case of two grains connected by two mesoscopic bonds. In other words, the weighting coefficients for pairwise conservative energy and force between 1−2 (bond), 1−3 (angle), and 1−4 (dihedral) partners were set at 0, 0.5, and 1, respectively. Intramolecular bonding and bending potentials were also included to account for chain stiffness. Bonds and angles were modeled by quartic and harmonic potentials, respectively. Four parameters are required to model the bond while only two are needed to model the angle formed by three consecutive monomers. B


Article

Macromolecules The dissipative force field uses the classical DPD component, leading to the following additional interaction between particles i and j

f jD→ i

2 l o o r zyz jij o o o−γ jjj1 − zzz vij ·eij ·eij for r < rd rd { =o m k o o o o o o0 for r ⩾ rd n

model (Multi-T) by combining the references trajectories in the temperature range between 200 and 400 K at constant pressure 0.1 MPa. The resulting potentials (bonded, angular, and nonbonded) roughly superimpose with the single-temperature ones obtained at 300 K. Contrary to their work, this functional form was not sufficient for polymers as we will discuss in the next section. To ensure transferability in temperature, we propose an extension of the method in which we set the conservative parameters xCMT to vary as a linear function of temperature; that is, xCMT = a + RTb, where a and b are two vectors of the same length than xCMT and R is the universal gas constant. The resulting force field is called Multi-LinT. The equations for the optimal a and b are given in the Appendix. We chose to keep γ̅MT independent of the temperature. Indeed, for the Single-T optimizations, we found that its value is almost constant, comprised between 21.57 and 23.54 kg mol−1 ns−1. For completeness, the conservative and dissipative optimal parameters obtained using the individual Single-T models and the Multi-T and Multi-LinT ones are shown in Figure S2. We also developed a multipressure (Multi-P) CG model by combining the eight coarse-grained reference trajectories in the pressure range between −10 MPa (traction) and 10 MPa (compression) at constant temperature T = 300 K. Finally, we propose a linear multitemperature, multipressure (Multi-LinT Multi-P) that combines the 17 reference trajectories used to build the Multi-LinT and the Multi-P CG models. Optimal Model Parameters. The optimal set of conservative parameters a̅ and b̅ for the Multi-LinT Multi-P CG model are summarized in Table 1. From them, one can compute the optimal parameters x̅CMT = a̅ + RTb̅ of the model at any temperature. The first four parameters {x̅0, x̅1, x̅2, x̅3} correspond to the weighting factors of the nonbonded force as expressed in eq 1. The next four parameters {x̅4, x̅5, x̅6, x̅7} model the quartic bond between two chemically linked monomers. The last two parameters {x̅8, x̅9} model the bending potential between three consecutive monomers. As mentioned above, bonds were modeled by a quartic potential that can be written as a linear combination of ri power law basis functions, directly optimized by the Bayesian method, such as

(2)

where γ is the friction coefficient, vij is the relative velocity vi − vj, and rd = 20 Å is the dissipative cutoff. The optimal γ was obtained from analysis of the coarse-grained reference trajectories ; using the Bayesian method. We also used the standard DPD random force to account for lost degrees of freedom eliminated by coarse-graining. The random interaction between particles i and j is

f jR→ i

l i o r yz 1 o j o o eij for r < rd o σ jjjj1 − zzzzθij r Δt =m d{ k o o o o o for r ⩾ rd o n0

(3)

where θij ∼ 5(0, 1) is a random Gaussian variate, σ is the amplitude of this interaction, and Δt = 50 fs is the DPD time step. The fluctuation−dissipation theorem imposes that σ = (2kBTγ)1/2 as shown by Español and Warren.26 f Dj→i and f Rj→i act between all pairs of particles including bond neighbors. Together, they form the DPD thermostat. Parametrization of the CG Model. The CG force field was optimized using statistical trajectory matching (Bayesian method). This method requires a reference trajectory from MD, recorded at time steps Δt. It optimizes the force field parameters so as to maximize the probability that the model will produce the reference trajectory. For DPD equations of motion, this probability can be computed analytically from eq 3. The optimal parameters can then be found explicitly or through iterative optimization. The details of the method can be found in refs 15 and 16. Single-temperature CG potentials (Single-T) at P = 0.1 MPa and single-pressure CG potentials (Single-P) at T = 300 K were optimized using each coarse-grained reference trajectory ; independently of the others. The resulting Single-T nonbonded potentials are plotted in Figure 2. As in the case of pentane,16 the higher the temperature, the more shifted the potential well to small distances. Following the path proposed by Solano Canchaya et al.,16 we developed a multitemperature CG

4

E bond(r ) =

∑ xi̅ + 3r i + E0 i=1

(4a)

with E0 an irrelevant energy constant. This expression can be rewritten in a more conventional form as 4

E bond(r ) =

∑ kbi(r − r0)i i=2

(4b)

with r0 the equilibrium bond distance and kbi the stiffness of the bond. r0 is the positive root of a third-order polynomial in r as expressed in eq 5a, which is the one with the physical meaning. 0 = 4x 7̅ r 3 + 3x6̅ r 2 + 2x5̅ r + x4̅

(5a)

kb4, kb3, and kb2 are calculated from the following system of equations: kb4 = x 7̅ kb3 = x6̅ + 4kb4r0 Figure 2. Single-temperature (Single-T) pairwise nonbonded conservative potential for cPB at 200, 250, 300, 350, and 400 K.

kb2 = x5̅ + 3kb3r0 − 6kb4r0 2 C

(5b) DOI: 10.1021/acs.macromol.8b02750 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Table 1. Coarse-Grain Linear Multitemperature, Multipressure (Multi-LinT Multi-P) Force Field Parameters for cis-1,4Polybutadiene (CPB)a

a

Their respective meaning is defined in the text. The effective parameters are obtained by computing x̅CMT = a̅ + RTb̅ for a chosen temperature in the range between 200 and 400 K. The units of RT are kJ mol−1. One bead represents one CPB repeating unit. Its mass is mb = 54.0904 g mol−1, and its charge is zero. The 1−2, 1−3, and 1−4 weighting factors of the nonbonded interaction are set to 0, 0.5, and 1, respectively. a̅ and b̅ are provided with more precision in the Supporting Information (cispb.ff).

intramolecular and intermolecular structures were investigated by computing the nonbonded radial distribution function (RDF) gnb(r) between two beads that are not directly connected together through a mesoscopic bond. These three distributions were computed over a 50 ns DPD production run and compared to their target functions (MDCG) obtained beforehand from the analysis of the coarse-grained molecular dynamics reference trajectories. Figure 3 shows these functions obtained at 300 K and 0.1 MPa. Even though the Bayesian method does not strictly focus on reproducing the local structure around beads, as opposed to more standard routines such as iterative Boltzmann inversion (IBI), we found a very good agreement between our model and the coarse-grained atomistic description. The average bond length lb obtained from DPD 3.87 Å is only slightly lower than the MDCG target value of 3.96 Å. Moreover, for both MD and DPD, the local structure of the chain is only barely modified by temperature. Using the CG model, the average angle between three beads is 126° at 400 K. This quantity decreases down to 117° at 200 K, traducing minor folding of the polymer chains with lower temperature. The RDFs are also only little impacted over the 200−400 K temperature range (see Figure S3) although the effect on density is huge (drop by 14% from 200 to 400 K, Figure 6a). This last result suggests that the RDF alone is not a good descriptor to build transferable potentials. At larger length scales, it is crucial to check that the CG polymer chains hold (statistically) their conformation described by the mean-square end-to-end vector ⟨Ree2⟩ and the meansquare radius of gyration ⟨Rg2⟩ (⟨ ···⟩ denotes the ensemble average over all chains and all time steps). ⟨Ree2⟩ is directly related to the packing length p defined such as27

Angles were modeled by a harmonic potential E bend(θ ) = ka(θ − θ0)2

(6)

where θ0 and ka are the equilibrium angle and angle stiffness, respectively. The two bending parameters are thus ka = x̅8 and θ0 = x̅9/ x̅8. Concerning the dissipative component, the optimal friction coefficient γ̅MT = 22.40 kg mol−1 ns−1 was deduced from eq 23 (see the Appendix). For training purposes, ready-to-use LAMMPS input files are provided in the Supporting Information. They comprise an initial polymer melt configuration (data.lmp), tabulated pairwise force and potential suitable at T = 300 K written in LAMMPS format (pair.table), and the input script to run LAMMPS (input.lmp). Speed-Up. Let us first discuss about computational costs. We compared the speed of execution required to produce 1 ns of simulation between DPD and MD using the optimal CG and united-atom descriptions, respectively. The simulation box consists in 90 monodisperse polymer chains of length 60 monomers each. As mentioned in the previous section, the number of particles is reduced by a factor 4 in the CG model, but the cutoff is increased from 14 to 20 Å. The real gain comes from the choice of the DPD time step for the integration of the equation of motion, which was allowed to be 25 times larger than the one used in MD (without any numerical instabilities). In the end, the DPD models are around 21 times faster than MD. We also noted that the relative performance between DPD and MD is independent of the number of processors used to make this benchmarking (from 1 to 729 CPUs).

■

VALIDATION OF THE CG MODEL Structural Properties. We now focus on the ability of the Multi-LinT Multi-P force field to conveniently reproduce the structural properties of the cPB. The intramolecular skeleton of the polymer chain, modeled by bonds and angles between two and three consecutives beads, has been investigated by means of the 1−2 bonding distribution gbond(r) and the bending distribution gbend(θ). The long-range

p=

M ⟨R ee ⟩ρ 5A 2

(7)

with M the polymer chain molar mass, ρ the equilibrium density, and 5A the Avogadro number. This quantity controls28 the viscoelastic properties of polymer melts as we will discuss in the next section. As a result, a good estimation of p is required to D


Article

Macromolecules

Figure 4, with Mij = mbNij the mass of the subchain between beads i and j. This makes possible the comparison with

Figure 4. Normalized mean-square internal distance as a function of the number of beads Nij inside the subchain between beads i and j for four polymer melts of different monodisperse molecular weights. The horizontal dotted line shows the experimental value of the normalized mean-square end-to-end distance ⟨Ree2⟩/M.28 The measurements were performed at T = 300 K and P = 0.1 MPa. For comparison, the dashed line represents the theoretical prediction for freely rotating chains (FRC)31 as expressed in eq 9.

experiments. The last point of the curve actually corresponds to ⟨Ree2⟩/M and is a bit overestimated with respect to the experimental value of 0.758 Å2 mol g−1. We highlight that this quantity is very sensitive since the end-to-end vector describes a broad Gaussian distribution as given in the Supporting Information in Figure S4 and in eq S6. A deviation of only 4% between simulation and experiment is thus first rate. On the other hand, from the first point of the curve, i.e., Nij = 2, we compute the average bond length lb between two neighboring repeating units (this can be achieved as well by looking at the 1− 2 bonding distribution gbond(r), Figure 3a). We then deduce the Kuhn length l 2 expressed as32

Figure 3. Bonding distribution (a), bending distribution (b), and nonbonded radial distribution function (c) at T = 300 K and P = 0.1 MPa. Comparison between coarse-grained MD (MDCG) and DPD. The DPD dynamics were performed using the generalized linear multitemperature, multipressure CG potential (Multi-LinT Multi-P).

expect realistic and quantitative rheological behavior. Using the Multi-LinT Multi-P model, we obtained under atmospheric conditions ⟨Ree2⟩/M = 0.787 ± 0.008 Å2 mol g−1 and p = 2.30 Å, while the analysis of the coarse-grained atomistic trajectory yields ⟨Ree2⟩/M = 0.724 ± 0.015 Å2 mol g−1 and p = 2.50 Å. The cPB polymer chains become therefore a bit more rigid in the coarse-grained model with respect to MD, but the matching is still outstanding since the mean square end-to-end vector is not directly taken as a target quantity when building the CG force field. Experimentally, p equals 2.44 Å for cPB28 (96% cis content, unprecised molecular weight and polydispersity), which is in good agreement with both MD and DPD simulations. Next, we verified that the ratio ⟨Ree2⟩/⟨Rg2⟩ approximately equals 6 for both the atomistic (6.14 ± 0.11 at 300 K and 0.1 MPa) and the CG model (6.13 ± 0.02 in the same conditions). We also checked that the probability distribution function of the end-toend distance is Gaussian29 (see Figure S4). These two verifications ensure the ideal behavior of the polymer chains as expected by Flory’s theory30 in the case of polymer melts. As stated above, the Multi-LinT Multi-P force field was optimized for chains of molecular weight of 3.2 kg mol−1 (60 beads per chain), which are long enough to exhibit self-similarity since the normalized mean-square internal distance ⟨Rij2⟩/Mij converges to a plateau value for large subchains, as described in

l2 =

⟨Ree 2⟩ (N − 1)lb

(8)

where N = 60 is the number of beads per chain and (N − 1)lb is the fully extended chain length (also called contour length). The Kuhn length represents the bond length of the equivalent freely jointed ideal chain and is characteristic of the microstructure. We report about 11% error between the atomistic (l 2 = 10.05 Å ) and the CG model (l 2 = 11.19 Å ). This deviation originates from both a 9% overestimation of the mean-square end-to-end vector and a 2% underestimation of the bond length (lMDCG = 3.96 Å, lCG b b = 3.87 Å). To go further, we now consider six simulation boxes consisting in periodic amorphous cells of monodisperse cPB melts with molecular weights M of 3.2, 6.5, 9.7, 16.2, 32.4, and 48.6 kg mol−1. These molecular weights correspond to 60, 120, 180, 300, 600, and 900 beads per chain, respectively. To have equivalent computing times, each box contains in total 5400 beads, and the number of chains is chosen accordingly. The six systems were prepared by putting cPB phantom chains together into a periodic box at a density of 0.92 g mol−1, i.e., at a density near its equilibrium using the CG model (see Figure 6). The initial conformations of the chains were controlled statistically E


Article

Macromolecules

polymer chains are nearly ideal and because no dihedral potential was included in the CG force field (although a subtle dihedral structure can be observed in the reference trajectories), the conformations of the subchains can be compared to the freely rotating chain (FRC) model,31 which predicts the evolution of the internal polymer chain structure as ÄÅ ÉÑ ÅÅ 2c(1 − ( −c)Bij ) ÑÑÑÑ 1−c 2 2Å Å ÑÑ + ⟨R ij ⟩ = Bij lb ÅÅÅ ÅÅ 1 + c Bij (1 + c)2 ÑÑÑÑ (9) ÅÇ Ö

using the coarse-grained bonding and bending distributions (independent of the chain length) as well as the target end-toend vector (Figure 4). This way, the initial structure of the melt is a good representation for cPB. In particular, the relaxation using the CG force field is facilitated. To avoid atom loss and bond breaking, excluded volume is introduced slowly between beads. We adapted the push-off method proposed by Auhl et al. in 200333 later extended for DPD.34 Instead of using nondiverging excluded volume potentials to separate overlapping beads, we directly took our more realistic Multi-LinT Multi-P model, in which the conservative nonbonded interaction was switched on progressively by linearly increasing the corresponding weighting factors as given in eq 1. The pushoff was realized over 100 ns in the constant-NVT ensemble. The systems are then relaxed for 100 ns in the constant-NPT ensemble before launching a 10 μs production run (2 × 108 steps) in the constant-NVT ensemble at the equilibrium volume (average box size over the last nanosecond of the equilibration). Figure 5 shows the end-to-end vector autocorrelation function

with Bij = Nij − 1 the number of bonds inside the subchain between beads i and j and c = ⟨cos θ⟩. No fitting parameters are needed because the two parameters of the model (average bond length and angle between two consecutive bonds) are directly estimated from the DPD trajectory (lb = 3.87 Å, ⟨cos(θ)⟩ = −0.504). Figure 4 shows perfect matching of our measurements with the theoretical prediction whatever the size of the subchain considered. These last findings confirm both the good equilibration (at all length scales) of the polymer melts and the transferability of the CG model in chain length. Transferability in Temperature and Pressure. Starting from a well-equilibrated coarse-grained atomistic configuration comprising 90 cPB chains of length 60 repeating units, DPD simulations in the constant-NPT ensemble were performed at different temperatures and pressures. Figure 6a shows the temperature dependence of the equilibrium density over the 200−400 K temperature range under isobaric conditions P = 0.1 MPa. In the same spirit, Figure 6b displays the equilibrium density ranging from −10 MPa (expansion) to 10 MPa (compression) under isothermal conditions T = 300 K. Above the glass transition, the density follows linear temperature and linear pressure dependencies that yield the thermal expansion coefficient αP and the isothermal compressibility χT expressed as

αP = Figure 5. Decorrelation of the end-to-end vector at T = 300 K in function of the molecular weight of the polymer chain. The full lines correspond to the fit using Kohlrausch−Williams−Watts stretched exponentials.35

1 ∂Vs Vs ∂T

χT = −

(10)

P

1 ∂Vs Vs ∂P

T

(11)

−1

with Vs = ρ the specific volume. First, the united-atom force field used in MD allows a quantitative prediction of the thermal response as αP = (6.87 ± 0.07) × 10−4 K−1 in accordance with the experimental value of (6.63−7.35) × 10−4 K−1 36−38 and previous MD simulations (αP = (7.1 ± 0.1) × 10−4 K−1 39 at T = 310 K). Under atmospheric conditions, the equilibrium density ρ = 0.92 g cm−3 obtained in MD is only 2% higher than the expected experimental value of 0.90 g cm−3,28,32 in accordance with earlier simulations using the same force field.18,19,39−41 The pressure behavior is more surprising. In the region between −10 and 10 MPa, χT = (7.00 ± 1.07) × 10−4 MPa−1 equals within the error bar the experimental value of 7.25 × 10−4 MPa−1 measured by DiBenedetto.36 Anderson et al.42 reported a lower experimental value of 3.6 × 10−4 MPa−1, which emphasizes the difficulty of measuring low compressibilities. Locally, between −1 and 1 MPa, the compressibility of the system increases 3-fold with χT reaching (26.86 ± 3.02) × 10−4 MPa−1. By construction, the CG force field has to be faithful to MD. As mentioned in the previous section, the Multi-T model is not able to reproduce the correct thermal response as αP = (1.22 ± 0.01) × 10−4 K−1 (see Figure 6a). On the contrary, the most advanced Multi-LinT model reproduces with excellent agreement the atomistic densities. Only a small discrepancy appears

obtained for these systems. The shortest chains lose their initial orientation after 100 ns, whereas the longest would require few hundreds of microseconds to fully relax. This result emphasizes the importance of building well the initial polymer melt. Furthermore, the characteristic shape of ⟨Rij2⟩/Mij as drawn in Figure 4 (steady increase with the number of beads) is a proof of both good relaxation of the polymer chains and sufficient sampling of the phase space of the system, even for high molecular weights up to 16.2 kg mol−1 (300 beads per chain). The local structure of the polymer chains is independent of their length, and the ratio ⟨Rij2⟩/Mij converges to the same plateau value. The statistics falls as Nij increases because less subchains can be considered for computation, which explains the slight fluctuations of the plateau when the subchain grows. For the two highest molecular weights (not shown), the statistics is even lower for Nij ⩾ 200, leading to greater uncertainty of ⟨Rij2⟩/Mij because the corresponding simulation boxes contain only few polymer chains. This uncertainty could be reduced by increasing the number of chains (created such as their overall conformations follow a Gaussian distribution) without any loss of computational time if increasing the number of CPUs in consequence, but it is not the purpose of this work. Because the F


Article

Macromolecules

Table 2. Thermal Expansion Coefficient and Isothermal Compressibility experiment MD DPDc

αP [× 10−4 K−1]

χT [× 10−4 MPa−1]

a

7.25b 7.00 ± 107 19.96 ± 19

6.63 6.87 ± 7 7.15 ± 12

Reference 37. bReference 36. cMulti-LinT Multi-P force field (this work). a

the plateau modulus Ge, through the following empirical universal laws: Me = B−2 p3 ρ 5A Ge =

(12)

4 ρ 5AkBT 4 = B2 kBTp−3 5 Me 5

(13)

with B = 0.0565 a universal proportionality constant fitted over 17 polymers, kB the Boltzmann constant, and T the temperature. We prefer the former definitions for Me and Ge rather than the one proposed later by the same authors32,43 as discussed by Larson et al.44 Using these relations, we predict Me = 2.1 kg mol−1 and Ge = 0.87 MPa for cPB with our CG model. In terms of average number of beads between two entanglement strands, we predict Ne = Me/mb = 39. The stress relaxation modulus G(t) is calculated from the autocorrelation of the stress tensor σ(t)45 using the Green− Kubo formula46,47 28

Figure 6. Transferability in temperature and pressure for various atomistic and coarse-grained force fields. The blue squares correspond to the MD equilibrium densities. The Single-T and Single-P yellow diamonds represent each one optimal set of CG potentials for given thermodynamic (P, T) conditions. (a) Temperature dependence of the density at P = 0.1 MPa. The Single-T200 (pink), Single-T300 (gray), and Single-T400 (purple) diamonds are the Single-T models optimized at 200, 300, and 400 K, respectively. These CG models were tried in the full temperature range between 200 and 400 K. The equilibrium densities obtained using the Multi-T (green squares) and Multi-LinT (mauve circles) are also plotted. The results using the generalized Multi-LinT Multi-P model are not displayed as they superimpose with the Multi-LinT ones. (b) Pressure dependence of the density at T = 300 K.

G (t ) =

for the lowest temperature 200 K where the CG system is 1% less compact. Concerning the mechanical response, the Multi-P potentials match nicely with MD in the pressure range between −1 and 1 MPa even though the optimization procedure took the full pressure interval between −10 and 10 MPa into account. This is quite startling since the equilibrium densities using the Single-P models coincide properly with their respective MD targets as drawn in Figure 6b. Nevertheless, the CG system remains relatively incompressible as χT = (19.86 ± 0.19) × 10−4 MPa−1 (in the Appendix, we show that Poisson’s ratio equals ν = 0.499). The generalized Multi-LinT Multi-P force field, which combines the two Multi-LinT and Multi-P models, preserves their respective characteristics without compromising accuracy. In other words, the Multi-LinT Multi-P model yields the same thermodynamical properties than its two parents models, ensuring the transferability in both temperature and pressure. The results are summarized in Table 2.

os os V ⟨σ (t + t ′):σ (t ′)⟩ 10kBT

(14)

Figure 7. Stress relaxation modulus in function of the molecular weight of the cPB chains. The horizontal dotted line corresponds to the experimental rubbery plateau modulus Gexp e = 0.76 MPa as proposed in ref 28. The dashed line represents a slope of −0.5 (Rouse regime). os

where σ (t ) is the traceless symmetric part of σ(t). Figure 7 shows G(t) as a function of the molecular weight of the cPB polymer chains. The shape of the curves can be interpreted in terms of polymer physics.29 As expected, G(t) decreases continuously from high values of the order of practically the GPa (solidlike behavior) until flow is observed (below the MPa). At early times, G(t) does not depend on the molecular

■

RHEOLOGY AND MECHANICS By analogy to elastomer networks, we know from Fetters et al.28 that the packing length p drives important viscoelastic properties of polymer melts, namely, the entanglement molar mass Me and G


Article

Macromolecules

reptation model proposed by de Gennes in 1971.51 It ends with an exponential decay (as for unentangled polymers) when the reptation time τrep ∝ M3 is reached. We do observe the beginning of the rubbery plateau in our simulations with long chains (Figure 7). This implies that the intermolecular interactions are strong enough to prevent bond crossing at entanglements under quiescent conditions. Although thebox size used in the simulations for the longest chain systems does not ensure the absence of self-interactions, it is larger than twice the average distance between entanglements and then only allow interactions between different subchains on a given chain. To confirm our observations and facilitate the determination of all these characteristic times, the long time dynamics of the cPB melts has been investigated. The mean-square displacement (MSD) of the center of mass of the polymer chains is presented in Figure 8a. Whatever the time scale considered, the dynamics of the polymer chain decreases with its molecular weight. For the shortest chains, we clearly distinguish the diffusion regime (scaling as t), starting from about 100 ns. For the longest chains, we observe a slope of 0.5 in the region where G(t) exhibits a plateau value (Figure 7). We now move to smaller length scales by looking at the MSD of the monomers. For the purpose of this analysis, we only considered in Figure 8b the monomers belonging to the middle of the polymer chains. We verified that the monomers at the end of a chain move faintly faster (in the subdiffusive regime) than the middle ones, as reported earlier. Interestingly, all the monomers follow exactly the same dynamics until they reach τe = 15 ± 4 ns. By carrying this value in Figure 7, we get a direct measurement of the plateau modulus Ge = 1.35 ± 0.30 MPa. We then deduce Me = 1.45 ± 0.41 kg mol−1 using the empirical relation 13 and Ne = 27 ± 8. The rubbery plateau modulus seems to be overestimated with respect to either the prediction using eq 13 or the experimental value of 0.76 MPa.28 The entanglement mass is therefore underestimated. However, other experimental studies have estimated Ge between 1.10 and 1.45 MPa depending on the molecular weight,52−54 depicting the difficulty of getting an accurate measure of Ge (log scale). We are thus within the experimental error. As discussed in the previous section, the cPB chains may be a little too stiff, which could also explain the deviation from experiment. Another source of uncertainties comes from the stress relaxation modulus itself because such data are usually noisy. We refer the reader to the work of Likhtman et al.48 for a more detailed discussion about this topic. The smooth G(t) curves displayed in Figure 7 were indeed obtained by filtering the raw G(t) (2 × 108 points each separated by 50 fs) with a sliding average on a logarithmic scale. We can still extract more information from the analysis of the MSD. The monomers are not aware they belong to a polymer chain of molar mass M before some time reaching either τR (unentangled) or τe (entangled). The limit between unentangled and entangled chains appears when τR = τe. At very short time scales, below about 300 fs (only few DPD time steps), the ballistic regime (scaling as t2) is observed, which illustrates the integration of the equation of motion of the beads. For t ⩽ τe, the motion of the monomers in the so-called subdiffusive regime is Rouse-like (scaling of the MSD as t0.5) between 0.1 and 1 ns, in accordance with the shape of G(t) (Figure 7) within the same temporal window. If τR > τe (entangled chains), the theory claims that the MSD of the beads scales as t0.25 for τe < t < τR and as t0.5 for τR < t < τrep. These power laws are all observed in Figure 8b. In particular, the MSD of the beads belonging to the longest polymer chains clearly exhibit the t0.25 power law characteristic

weight of the chain since all the curves collapse onto one master curve. This result also demonstrates the good preparation of the systems. Contrary to the work of Likhtman et al.,48 we do not report any oscillation at very low time scales arising from bond relaxations. We believe that these oscillations are screened because of the use of a large DPD time step. A little bump is observed at around few hundreds of femtoseconds. This event corresponds to the transition between the ballistic and the first relaxation regime of the beads (see Figure 8b). G(t) scales as

Figure 8. Mean-square center-of-mass displacement (a) and meansquare middle chain monomer displacement (b) as a function of the molecular weight of the chain. The two plots are displayed using the same log scale for easy comparison. We also precise here that the dynamics of the beads measured in DPD is about 2 times faster (in the subdiffusive regime) than the reference obtained from the coarsegrained reference trajectories (not shown).

t−0.5 between 0.1 and 1 ns, which corresponds to a Rouse-like behavior.49 Then, the decrease of G(t) slows down as the molar mass of the chain increases. This slowdown is characteristic of entangled polymer dynamics and becomes even more visible if plotting t0.5G(t) (see Figure S5a). For unentangled polymers, G(t) decays exponentially beyond the Rouse time of the polymer chain τR ∝ M2 (longest Rouse mode). For entangled polymers, G(t) reaches the so-called rubbery plateau beyond the Rouse time of an entanglement strand τe ∝ Me2 where Me is the entanglement mass. Because the rubbery plateau decreases slightly with time due to tube length fluctuations,50 the plateau modulus Ge is formally defined as the value of the stress relaxation modulus at τe.29 This regime can be described by the H


Article

Macromolecules

because the MSD is averaged over all beads, allowing for greater statistics even on a long time range. Furthermore, when reaching the diffusion regime (only clearly visible for the shortest chains), the motion is driven by the collective motion of the beads inside a molecule, i.e., by the center of the polymer chain as shown in Figure 8a. In this regime, the MSD of the beads (regardless of their position inside the polymer chain) superimpose with the one of the center of mass of the chain. Accordingly, the combination of Figures 8a and 8b allows a more accurate estimation of the relaxation time of the whole chain (τR or τrep depending on the entanglement rate) that delimits the start of the diffusion regime. The self-diffusion coefficient of the polymer chains + is then computed from the MSD with the relation

of reptation and thus confirms our observations. We highly recommend plotting t−0.5 MSD(t) to help defining τe, τR, and τrep (see Figure S5b). By combining Rouse and reptation theories, one can demonstrate that ij M yz τR = jjj zzz j Me z τe k {

2

(15)

ij M yz = jjj zzz j Me z τe (16) k { The characteristic ratios τR/τe and τrep/τe are plotted in Figure 9a as a function of the molecular weight of the polymer chain. τrep

3

+=

1 d lim ⟨(r(t ) − r(0))2 ⟩ 6 t →+∞ dt

(17)

Figure 9b reports + as a function of M at T = 300 K. We should point out here that + is overestimated for the chain with molecular weights beyond 16.2 kg mol−1 because the transition between the reptation and the diffusion regimes occurs after 10 μs of simulation. Interestingly, the self-diffusion coefficient is accessible experimentally by the means of neutron scattering,55 recoil spectrometry,56 scanning infrared microscopy57, and NMR58 techniques. The data obtained for polybutadiene at T = 413 K, compiled by Lodge in 1999,59 are represented in Figure 9b (molar mass dispersity Đ < 1.06). Assuming the validity of Einstein relation, + ∝ T , which allows us to rescale + by a factor 413/300 = 1.38 for comparison with experiment. We observe an excellent matching between experiment and simulation. The diffusion of the polymer in the simulation is only a bit faster, and the scaling of + as M−2.3, which is experimentally observed, is roughly verified. The entanglement mass was also estimated through primitive path analysis (PPA).60 A homemade implementation of the Zcode61 was run on 9000 sample configurations from the 10 μs trajectories for each chain length. The average entanglement mass ⟨Me⟩ was computed from the lengths of the primitive paths Lpp. Two estimators were used, which respectively slightly overand underestimate ⟨Me⟩:62 ij ⟨Lpp2⟩ yz zz ⟨Me⟩+ = (N − 1)mbjjjj − 1 zz j ⟨R ee 2⟩ z k {

−1

Figure 9. (a) Characteristic ratios τR/τe and τrep/τe as a function of the molar mass M of the polymer chains. The horizontal dotted line corresponds to the limit of the simulation (10 μs) normalized by τe. The other dotted line represents a fit of eq 15, which is the expected evolution of τR/τe = (M/Me)2 according to Rouse and reptation theories. τrep is not easily accessible but should follow the expected behavior τrep/τe = (M/Me)3. (b) Diffusion coefficient of cPB obtained from DPD simulations at T = 300 K (rescaled at 413 K assuming + ∝ T ) and measured experimentally at 413 K as a function of the molecular weight of the polymer. The dotted and dashed lines represent slopes of −2 (Rouse) and −2.3 (experimental observations55−59).

⟨Me⟩− = (N − 1)mb

(18)

⟨R ee 2⟩ ⟨Lpp⟩2

(19)

The mean of both estimates was calculated and found to not depend significantly on the chain length, yielding ⟨Me⟩ = 1.82 ± 0.18 kg mol−1. This value is intermediate between estimations from the MSD and from the G(t) analyses. All estimates of ⟨Me⟩ are summarized in Table 3.

■

CONCLUSIONS In conclusion, we propose a realistic coarse-grain (CG) force field to model cis-1,4-polybutadiene (cPB) polymer melts. The parametrization of this new force field was conducted by adopting the Bayesian bottom-up strategy based on trajectory matching, that integrates the development of CG potentials from a set of high-resolution atomistic reference trajectories at different thermodynamic states. The transferability in temperature between 200 and 400 K and pressure between −10 MPa

Because τrep is too close to the simulation time (or higher), its precise determination was not accessible (except for the shortest chains). The simulation confirms the validity of relation 15, from which we deduce Me = 2.44 ± 0.75 kg mol−1 (Ne = 45 ± 14). Me is in excellent agreement with the experimental value of 2.3 kg mol−1 (Ne = 43).28 From a technical point of view, this method based on the dynamics of the beads is also the most accurate I


Article

Macromolecules

target property, we obtained a quantitative value of the plateau modulus in accordance with experiment. The mean-square displacement of the polymer chains and the beads corroborate our observations. Their extensive analysis showed transitions from Rouse to reptation and diffusion regimes, from which we could extract a precise measurement of the entanglement mass. In summary, starting from the unique knowledge of the microstructure of the cPB, we made the bridge between chemistry and rheology. From a mechanical point of view, since the compressibility of the rubber is well reproduced, we expect this CG model to be useful if probed under mechanical sollicitations. At this level, we are now able to provide relevant mesoscopic and macroscopic quantities such as the rubbery plateau modulus or the entanglement mass for a given material of interest. In the future, we expect to give some new quantitative insights into structure/property relationship of polymer-based materials. Finally, this work opens the door for future developments of quantitative CG force fields specific to the chemical nature of the species considered.

Table 3. Average Entanglement Mass Estimated from Our Simulations Using Different Methodsa method

⟨Me⟩ [kg mol−1]

G(t) (plateau) MSD (Rouse) PPA

1.45 ± 41 2.44 ± 75 1.82 ± 18

a

The value estimated via G(t) depends on the definition of the plateau. The weighted average of all the estimations is ⟨Me⟩ = 1.79 ± 0.16 kg mol−1 while the experiments28 give 2.3 kg mol−1.

(traction) and 10 MPa (compression) is ensured. The corresponding thermal expansion coefficient and isothermal compressibility were found to be consistent with both molecular dynamics (MD) and experiment. The structural properties reproduced by the coarse-grained model are representative of cPB. Even though the Bayesian optimization is not directly based on the reproduction of the local structure, in contrast with more commonly used methods such as iterative Boltzmann inversion (IBI), the matching of the bonding and bending distributions in the DPD simulations with respect to their coarse-grained reference counterparts (MDCG) were found to be very good, especially since their modeling is simply ensured by quartic and harmonic potentials, respectively. The nonbonded radial distribution function, which quantifies the longrange intramolecular and intermolecular structuring, is also well reproduced in the full temperature range, while keeping the correct pressure behavior. In addition, we showed that this new CG model can conveniently reproduce the conformations of the cPB polymer chains. The distribution of the end-to-end vector, which follows the expected Gaussian behavior for ideal chains in a melt, superimposes with the one obtained from MD simulations. Finally, we checked the transferability of the CG model in chain length through the mean-square internal distance, which also matches well with theory using the freely rotating chain (FRC) model. In this CG model, one coarse grain represents one cPB monomer. By dividing the number of particles by 4 with respect to the united-atom description, we reached an overall speed-up of 21, while keeping both the thermodynamic and the structural properties of cPB. The dynamics of the dissipative model was found to be only slightly faster than the coarse-grained atomistic reference. In this study, the dissipative cutoff rd was chosen equal to the conservative cutoff rc = 20 Å, but there is no general rule for it, and these two nonlinear parameters could have been chosen independently. In addition, in a recent paper,24 we introduced the concept of reduced mass to conveniently reproduce the short time dynamics. This mass correction relies on the simulation algorithm (equipartition of energy, velocityVerlet) and depends on the coarse time step one wants to use. Further refinements of the CG model could thus be achieved by choosing a reduced mass of the bead of 82.235 g mol−1 instead of 54.0904 g mol−1 (with a DPD time step of Δt = 50 fs). Nevertheless, the current model already fits quantitatively to most of the physical properties. We took advantage of the DPD efficiency to perform 10 μs long DPD runs to investigate the rheological properties of cPB. Using our CG model, we were able to quantify the viscoelastic properties of monodisperse cPB polymer melts for various molecular weights. The transition from unentangled to entangled polymer chains has been clearly identified. The rubbery plateau modulus is accessible after 3 × 106 DPD time steps. Although the CG force field was not parametrized on this

■

APPENDIX

Optimization of Parameters with a Linear Dependence on Temperature

If we impose that the model parameters depend linearly on temperature, xCMT = a + RTb, the Bayesian minimization procedure becomes ∂3(x) = ∂a

∑

∂3(x) = ∂b

C ref ∑ ∑ A t*Bt+(A t xMT ̅ − Ft ) = 0

;

;

1 T;

C ref ∑ A t*Bt+(A t xMT ̅ − Ft ) = 0

(20)

t∈;

(21)

t∈;

where 3(x) is the negative log-likelihood of the parameters x (to be minimized) and T; is the temperature of trajectory ; . The matrices At, Fref t , Vt, and Bt are defined in ref 15. The superscripts ∗ and + mean respectively matrix transposition and matrix pseudoinverse. Using x̅CMT = a̅ + RTb̅, eqs 20 and 21 can be rewritten in a matrix format, yielding the optimal a̅ and b̅ ij y 1 jj∑ ∑ A t*Bt+A t R ∑ ∑ A t*Bt+A t zzzz j zz ij a ̅ yz jjj ; T; t ∈ ; ; t∈; zz jj zz = jj zz j b ̅ z jj z + + k { jj∑ ∑ A t*Bt A t R ∑ T; ∑ A t*Bt A t zzz jjj zz ; t∈; k ; t∈; { ij y 1 ref + jj∑ ∑ A t*Bt Ft zzzz jj zz jj ; T; t ∈ ; zz jj zz jj jj ∑ ∑ A *B+F ref zzz t t t jjj zzz (22) k ; t∈; { +

The optimal friction parameter γ̅MT is then obtained by solving ∂3(x) = 0. For multitemperature reference trajectories, γ̅MT is ∂γ given by 1 nt

∑

∑

1 T;

2 0 = γMT ̅

−

1 nt

;

;

1 T;

2NR ∑ V t*Bt Vt + γMT ̅ Δt

t∈;

∑ (Ftref )*Bt+(Ftref

− A t a ̅ + RT; A t b ̅ )

t∈;

(23) J


Article

Macromolecules with N the number of degrees of freedom and nt = ∑; ∑t ∈ ; 1 the number of configurations. Equation 23 shows that γ̅MT is just the positive root of a second-order scalar polynomial in γ̅MT.

ORCID

Bulk Modulus and Poisson’s Ratio

Notes

The mechanical properties discussed in the main text can be formulated in terms of equivalent quantities (bulk modulus, Poisson’s ratio, Young’s modulus). The bulk modulus K, which quantifies the resistance to isotropic compression, is simply the reciprocal of the isothermal compressibility χT.

The authors declare no competing financial interest.

K=

1 χT

A. Dequidt: 0000-0003-1206-1911 P. Malfreyt: 0000-0002-3710-5418

■

ACKNOWLEDGMENTS This work was partially funded by the Investissements d’Avenir program “Développement de l’Economie Numérique” through the SMICE project.

■

(24)

Unsurprisingly (same discrepancy for χT), the bulk modulus K = 501 ± 5 MPa is 3 times lower than the experimental value of 1.38 GPa.36 Poisson’s ratio ν, which measures volume changes due to stretching or compression of a material, can be deduced from K and Ge prior to running any further simulation. ν=

3K − 2Ge 2(3K + Ge)

(25)

One can notice that as long as K ≫ Ge, ν equals 0.5, traducing the actual incompressibility of the material. This behavior is expected for rubbers, since their volume remains constant during uniaxial traction tests. The application of eq 25 yields ν = 0.499. From a mechanical point of view, the CG model preserves the incompressibility of the cPB. Young’s modulus E, which quantifies the stiffness of a material in the linear elasticity regime (small deformations) under uniaxial deformation, is then obtained as well from E = 2(1 + ν)Ge = 3(1 − 2ν)K

(26)

The application of eq 26 yields E = 4.05 ± 0.90 MPa.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.8b02750.

■

REFERENCES

(1) Sutter, H. The Free Lunch Is Over: A Fundamental Turn Toward Concurrency in Software. Dr. Dobb’s Journal 2005, 10. (2) Goujon, F.; Malfreyt, P.; Tildesley, D. J. Interactions Between Polymer Brushes and a Polymer Solution: Mesoscale Modelling of the Structural and Frictional Properties. Soft Matter 2010, 6, 3472−3481. (3) Ghoufi, A.; Emile, J.; Malfreyt, P. Recent Advances in Many Body Dissipative Particles Dynamics Simulations of Liquid-Vapor Interfaces. Eur. Phys. J. E: Soft Matter Biol. Phys. 2013, 36, 10. (4) Li, Y.; Abberton, B. C.; Kröger, M.; Liu, W. K. Challenges in Multiscale Modeling of Polymer Dynamics. Polymers 2013, 5, 751− 832. (5) Español, P.; Warren, P. B. Perspective: Dissipative Particle Dynamics. J. Chem. Phys. 2017, 146, 150901. (6) Padding, J. T.; Briels, W. J. Systematic Coarse-Graining of the Dynamics of Entangled Polymer Melts: The Road from Chemistry to Rheology. J. Phys.: Condens. Matter 2011, 23, 233101. (7) Maurel, G.; Goujon, F.; Schnell, B.; Malfreyt, P. Prediction of Structural and Thermomechanical Properties of Polymers from Multiscale Simulations. RSC Adv. 2015, 5, 14065−14073. (8) Marrink, S. J.; de Vries, A. H.; Mark, A. E. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 108, 750−760. (9) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24, 1624−1636. (10) Maurel, G.; Schnell, B.; Goujon, F.; Couty, M.; Malfreyt, P. Multiscale Modeling Approach toward the Prediction of Viscoelastic Properties of Polymers. J. Chem. Theory Comput. 2012, 8, 4570−4579. (11) de Oliveira, T. E.; Netz, P. A.; Kremer, K.; Junghans, C.; Mukherji, D. C-IBI: Targeting Cumulative Coordination Within An Iterative Protocol to Derive Coarse-Grained Models of (MultiComponent) Complex Fluids. J. Chem. Phys. 2016, 144, 174106. (12) Ohkuma, T.; Kremer, K. Comparison of Two Coarse-Grained Models of Cis-Polyisoprene With and Without Pressure Correction. Polymer 2017, 130, 88−101. (13) Hijón, C.; Español, P.; Vanden-Eijnden, E.; Delgado-Buscalioni, R. Mori-Zwanzig Formalism as a Practical Computational Tool. Faraday Discuss. 2010, 144, 301−322. (14) Lemarchand, C. A.; Couty, M.; Rousseau, B. Coarse-Grained Simulations of Cis- and Trans-Polybutadiene: A Bottom-up Approach. J. Chem. Phys. 2017, 146, 074904. (15) Dequidt, A.; Canchaya, J. G. S. Bayesian Parametrization of Coarse-Grain Dissipative Dynamics Models. J. Chem. Phys. 2015, 143, 084122. (16) Solano Canchaya, J. G.; Dequidt, A.; Goujon, F.; Malfreyt, P. Development of DPD Coarse-Grained Models: From Bulk to Interfacial Properties. J. Chem. Phys. 2016, 145, 054107. (17) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (18) Smith, G. D.; Paul, W. United Atom Force Field for Molecular Dynamics Simulations of 1,4-Polybutadiene Based on Quantum Chemistry Calculations on Model Molecules. J. Phys. Chem. A 1998, 102, 1200−1208.

United-atom force field parameters for cPB; Bernstein bases polynomials Pi used to build the pairwise conservative nonbonded force; optimal parameters versus temperature for the Single-T, Multi-T, and Multi-LinT models; nonbonded radial distribution function for various temperatures at constant pressure P = 0.1 MPa obtained in DPD and MDCG; distribution of the end-toend vector as a function of the molecular weight of the polymer chain for monodisperse cPB melts; stress relaxation modulus and mean-square displacement of the monomers in the middle of the polymer chains normalized by the Rouse behavior (PDF) Tabulated CG force and potential written in LAMMPS format (pair.table) for T = 300 K; initial polymer melt configuration (data.lmp); input script to run LAMMPS (input.lmp); text file (cispb.ff) containing a̅ and b̅ vectors at high precision (ZIP)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. K


Article

Macromolecules

(43) Fetters, L. J.; Lohse, D. J.; Milner, S. T.; Graessley, W. W. Packing Length Influence in Linear Polymer Melts on the Entanglement, Critical, and Reptation Molecular Weights. Macromolecules 1999, 32, 6847−6851. (44) Larson, R. G.; Sridhar, T.; Leal, L. G.; McKinley, G. H.; Likhtman, A. E.; McLeish, T. C. B. Definitions of Entanglement Spacing and Time Constants in The Tube Model. J. Rheol. 2003, 47, 809−818. (45) Thompson, A. P.; Plimpton, S. J.; Mattson, W. General Formulation of Pressure and Stress Tensor for Arbitrary Many-Body Interaction Potentials Under Periodic Boundary Conditions. J. Chem. Phys. 2009, 131, 154107. (46) Daivis, P. J.; Evans, D. J. Comparison of Constant Pressure and Constant Volume Nonequilibrium Simulations of Sheared Model Decane. J. Chem. Phys. 1994, 100, 541−547. (47) Lee, W. B.; Kremer, K. Entangled Polymer Melts: Relation between Plateau Modulus and Stress Autocorrelation Function. Macromolecules 2009, 42, 6270−6276. (48) Likhtman, A. E.; Sukumaran, S. K.; Ramirez, J. Linear Viscoelasticity from Molecular Dynamics Simulation of Entangled Polymers. Macromolecules 2007, 40, 6748−6757. (49) Rouse, P. E. A. Theory of the Linear Viscoelastic Properties of Dilute Solutions of Coiling Polymers. J. Chem. Phys. 1953, 21, 1272− 1280. (50) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Clarendon Press: 1988. (51) de Gennes, P. G. Reptation of a Polymer Chain in the Presence of Fixed Obstacles. J. Chem. Phys. 1971, 55, 572−579. (52) Raju, V. R.; Menezes, E. V.; Marin, G.; Graessley, W. W.; Fetters, L. J. Concentration and Molecular Weight Dependence of Viscoelastic Properties in Linear and Star Polymers. Macromolecules 1981, 14, 1668−1676. (53) Carella, J. M.; Graessley, W. W.; Fetters, L. J. Effects of Chain Microstructure on the Viscoelastic Properties of Linear Polymer Melts: Polybutadienes and Hydrogenated Polybutadienes. Macromolecules 1984, 17, 2775−2786. (54) Baurngaertel, M.; De Rosa, M. E.; Machado, J.; Masse, M.; Winter, H. H. The Relaxation Time Spectrum of Nearly Monodisperse Polybutadiene Melts. Rheol. Acta 1992, 31, 75−82. (55) Bartels, C. R.; Crist, B.; Graessley, W. W. Self-Diffusion Coefficient in Melts of Linear Polymers: Chain Length and Temperature Dependence for Hydrogenated Polybutadiene. Macromolecules 1984, 17, 2702−2708. (56) Crist, B.; Green, P. F.; Jones, R. A. L.; Kramer, E. J. Self-Diffusion of Hydrogenated Polybutadiene by Forward Recoil Spectroscopy. Macromolecules 1989, 22, 2857−2858. (57) von Seggern, J.; Klotz, S.; Cantow, H. J. Reptation and Constraint Release in Linear Polymer Melts: An Experimental Study. Macromolecules 1991, 24, 3300−3303. (58) Pearson, D. S.; Fetters, L. J.; Graessley, W. W.; Ver Strate, G.; von Meerwall, E. Viscosity and Self-Diffusion Coefficient of Hydrogenated Polybutadiene. Macromolecules 1994, 27, 711−719. (59) Lodge, T. P. Reconciliation of the Molecular Weight Dependence of Diffusion and Viscosity in Entangled Polymers. Phys. Rev. Lett. 1999, 83, 3218−3221. (60) Everaers, R.; Sukumaran, S. K.; Grest, G. S.; Syaneborg, C.; Sivasubramanian, A.; Kremer, K. Rheology and Microscopic Topology of Entangled Polymeric Liquids. Science 2004, 303, 823−826. (61) Kröger, M. Shortest Multiple Disconnected Path for the Analysis of Entanglements in Two- and Three-Dimensional Polymeric Systems. Comput. Phys. Commun. 2005, 168, 209−232. (62) Hoy, R. S.; Foteinopoulou, K.; Kröger, M. Topological Analysis of Polymeric Melts: Chain-Length Effects and Fast-Converging Estimators for Entanglement Length. Phys. Rev. E 2009, 80, 031803.

(19) Tsolou, G.; Mavrantzas, V. G.; Theodorou, D. N. Detailed Atomistic Molecular Dynamics Simulation of cis-1,4-Poly(butadiene). Macromolecules 2005, 38, 1478−1492. (20) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to An External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (21) Schneider, T.; Stoll, E. Molecular-Dynamics Study of a ThreeDimensional One-Component Model for Distortive Phase Transitions. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 17, 1302−1322. (22) Dünweg, B.; Paul, W. Brownian Dynamics Simulations Without Gaussian Random Numbers. Int. J. Mod. Phys. C 1991, 02, 817−827. (23) Makhiyanov, N.; Temnikova, E. V. Glass-Transition Temperature and Microstructure of Polybutadienes. Polym. Sci., Ser. A 2010, 52, 1292−1300. (24) Kempfer, K.; Devémy, J.; Dequidt, A.; Couty, M.; Malfreyt, P. Development of Coarse Grain Models for Polymers by Trajectory Matching. ACS Omega 2019, under revision. (25) Natanson, I. P. Constructive Function Theory I: Uniform Approximation; Ungar Pub Co.: 1964; Vol. I. (26) Español, P.; Warren, P. Statistical Mechanics of Dissipative Particle Dynamics. Europhys. Lett. 1995, 30, 191. (27) Witten, T. A.; Milner, S. T.; Wang, Z.-G. In Multiphase Macromolecular Systems; Culbertson, B. M., Ed.; Springer: New York, 1989. (28) Fetters, L. J.; Lohse, D. J.; Richter, D.; Witten, T. A.; Zirkel, A. Connection between Polymer Molecular Weight, Density, Chain Dimensions, and Melt Viscoelastic Properties. Macromolecules 1994, 27, 4639−4647. (29) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press: New York, 2003. (30) Flory, P. J. The Configuration of Real Polymer Chains. J. Chem. Phys. 1949, 17, 303−310. (31) Flory, P. J.; Volkenstein, M. Statistical Mechanics of Chain Molecules. Biopolymers 1969, 8, 699−700. (32) Fetters, L. J.; Lohse, D. J.; Colby, R. H. In Physical Properties of Polymers Handbook; Mark, J. E., Ed.; Springer: New York, 2007; pp 447−454. (33) Auhl, R.; Everaers, R.; Grest, G. S.; Kremer, K.; Plimpton, S. J. Equilibration of Long Chain Polymer Melts in Computer Simulations. J. Chem. Phys. 2003, 119, 12718−12728. (34) Sliozberg, Y. R.; Andzelm, J. W. Fast Protocol for Equilibration of Entangled and Branched Polymer Chains. Chem. Phys. Lett. 2012, 523, 139−143. (35) Milano, G.; Müller-Plathe, F. Mapping Atomistic Simulations to Mesoscopic Models: A Systematic Coarse-Graining Procedure for Vinyl Polymer Chains. J. Phys. Chem. B 2005, 109, 18609−18619. (36) DiBenedetto, A. T. Molecular Properties of Amorphous hHgh Polymers. I. A Cell Theory for Amorphous High Polymers. J. Polym. Sci., Part A: Gen. Pap. 1963, 1, 3459−3476. (37) Paul, D. R.; DiBenedetto, A. T. Diffusion in Amorphous Polymers. J. Polym. Sci., Part C: Polym. Symp. 1965, 10, 17−44. (38) Barlow, J. W. Measurement of the PVT Behavior of cis-1,4Polybutadiene. Polym. Eng. Sci. 1978, 18, 238−245. (39) Tsolou, G.; Harmandaris, V. A.; Mavrantzas, V. G. Temperature and Pressure Effects on Local Structure and Chain Packing in cis-1,4Polybutadiene from Detailed Molecular Dynamics Simulations. Macromol. Theory Simul. 2006, 15, 381−393. (40) Sharma, P.; Roy, S.; Karimi-Varzaneh, H. A. Validation of Force Fields of Rubber through Glass-Transition Temperature Calculation by Microsecond Atomic-Scale Molecular Dynamics Simulation. J. Phys. Chem. B 2016, 120, 1367−1379. (41) Smith, G. D.; Paul, W.; Monkenbusch, M.; Willner, L.; Richter, D.; Qiu, X. H.; Ediger, M. D. Molecular Dynamics of a 1,4Polybutadiene Melt. Comparison of Experiment and Simulation. Macromolecules 1999, 32, 8857−8865. (42) Anderson, J. E.; Davis, D. D.; Slichter, W. P. Pressure Dependence of Molecular Motion in Some Elastomers. Macromolecules 1969, 2, 166−169. L


Realistic Coarse-Grain Model of - ACS Publications - American

Recommend Documents