Prediction of structural and thermomechanical properties of polymers from multiscale simulations

Abstract

We report mesoscale simulations of polymer melts and crosslinked polymer networks by using realistic coarse-grained (CG) models that are developed from atomistic simulations of polymer melts. We apply this multiscale strategy to different polymers in order to predict quantitatively some structural and thermomechanical properties such as the melt density, the end-to-end distance, the entanglement mass and the plateau modulus. The temperature dependence of the CG models is investigated through the calculation of the melt specific volumes at different temperatures and the calculation of the isothermal compressibility gives some insight into the pressure transferability of the CG models. We also show that the CG models can be applied successfully to high molecular weight chains. We test the performance of the CG models by calculating directly the plateau modulus of a crosslinked PIB network from mesoscopic simulations under a tensile stress. We compare the value of the plateau modulus with that calculated from the autocorrelation of the stress tensor during equilibrium simulations.

---

1 Introduction

The improvement of the desired properties and performance of a polymeric material requires investigating the relationship between the structure at the atomic or molecular length and the key physical and chemical macroscopic properties. Although molecular theories developed in polymer science have made paramount progress in relating the structure of polymer chains to their macroscopic properties, such theoretical methods suffer from a limiting predictive power due to the large variety of chemical constituents and molecular architectures that are at the origin of the viscoelastic properties. Molecular modelling^1–3 has become an efficient tool for checking the theoretical predictions of properties, for giving a molecular description of the analysis of experimental work and for designing advanced materials with specific properties. However, many problems at the leading edge of materials science involve collective phenomena that occur over a range of time and length scales that are difficult to capture in atomistic simulations. As a result, the ability to perform computer simulations of materials over length scales that are relevant to experiments represents a grand challenge in computational material science. Different approaches have been developed to model the matter at two or more different length and time scales and participate to the overall strategy of multiscale modelling.^4–16

Indeed, a powerful investigation of the relationship between molecular structure and mechanical properties of polymer melts requires statistical averages on well-equilibrated melts of very long chains. For a polymer chain of length N, the longest relaxation of an entangled polymer^2,17 scales as N^3.4 with a computational cost growing as N.⁴ Atomistic simulations require solving the equations of motion on a timescale of 10⁻¹⁵ s and a length of Angstroms whereas collective relaxations can occur on a scale of micrometers and a time that can exceed 1 ms (see Fig. S1 (ESI†) for a comparison between atomistic and mesoscopic simulations for the calculation of the end-to-end vector autocorrelation function). To simulate materials on larger time scales, there is no alternative but to simplify the model. These simplified/coarse-grained (CG) models allow for longer length and time scales than atomistic models that are unable to reach the complete relaxation of the polymer melts^7,8,18 with a reasonable computational effort. A number of CG models^19–28 have been applied in the past to reproduce universal (scaling laws) properties of polymer melts but the quantitative prediction of properties as a function of the polymer chemistry remains much more challenging. The key problem is how to simplify the model without loosing the essential physics and being able to relate the parameters of the CG model to molecules of specific chemistry. Top-down and bottom-up parameterization schemes represent then two different alternatives to derive CG models. Top-down approaches derive parameters from macroscopic properties (compressibility, diffusion)^29–31 and requires the use of a well-developed experimental database. Bottom-up approaches use the configurations at the atomistic level to develop interaction forces and parameters for mesoscopic model.^{5,9,16,18,32–34} It represents an attractive alternative for designing new polymeric materials from the structure–properties relationship.

Within the bottom-up approach, the Iterative Boltzmann Inversion (IBI) method is an iterative method for potential inversion from distribution functions using the potential of mean force approach.³⁵ CG intramolecular and intermolecular potential functions can be then deduced from appropriate distribution functions of atomistic configurations.^36–39 This method has been successfully applied to a variety of polymer melts formed by vinyl polymer chains,⁴⁰ dendrimers,⁴¹ polystyrene^7,8,32 and polyethylene chains.^37–39 These coarse-grained models have been implemented either in Molecular Dynamics^7,32,40–42 or in Dissipative Particle Dynamics (DPD)^8,37–39 methods. There are some major advantages to use the DPD method: (i) in addition to the conservative force, the dissipative and random forces are short-ranged and pairwise additive so that the hydrodynamic interactions are preserved; (ii) the random and dissipative forces are independent of the conservative forces and are coupled to each other through the fluctuation–dissipation theorem to act as a thermostat. However, the use of IBI method for building CG potentials addresses a number of fundamental issues:^7,40 (i) we have no guarantee that the potentials developed from atomistic configurations to match the distributions functions are able to reproduce all the properties of the original system; (ii) the way of grouping the atoms into the CG element and the degree of coarse-graining are not unique; (iii) the dynamics were often found to be significantly faster than the atomistic simulations.

In this paper, we propose to combine the IBI approach with the DPD method to reproduce thermomechanical and structural properties of different polymer melts (see Fig. 1). In a first step, we develop CG potentials for three polymers (cis-1,4-PB, PDMS and PIB) to calculate the density, end-to-end distance, entanglement mass and plateau modulus of melts. We examine the temperature and pressure transferabilities throughout the calculation of the coefficients of thermal expansion and isothermal compressibility, respectively. We also investigate the transferability to polymer chains of higher molecular weight. In a second step, we apply these potentials to different methodologies for the calculation of the plateau modulus of a crosslinked PIB network. We opt for a crosslinked network in order to avoid too large system-sizes that are expected to show a plateau in polymer melts.^43–46 The plateau modulus is then calculated from non-equilibrium simulations by means of stress–strain curve and by equilibrium simulations using the stress tensor elements.

Fig. 1 cis-1,4-Polybutadiene (cis-1,4-PB), poly(dimethylsiloxane) (PDMS), polyisobutylene (PIB) chemical structures along with the corresponding beads formed by five monomers.

The outline of this paper is organized as follows. In Sec. II, we present the DPD method and the coarse graining procedure. Sec. III presents the results of this work concerning the calculation of the thermomechanical and structural properties of PIB, PB and PDMS melts and of the shear modulus of a crosslinked network of PIB. The main conclusions of this work are summarized in Sec. IV.

2 Developments of the coarse-grained (CG) potentials

2.1 Atomistic molecular dynamics (MD) simulations

The atomistic simulations of the cis-1,4-PB, PDMS and PIB polymers were performed using the all-atom representation of the COMPASS force field⁴⁷ and the Forcite module from MATERIALS STUDIO Accelrys package.⁴⁸ The bulk atomistic configurations consisted of 40 chains of 200 monomers. Constant-NPT simulations were carried out at 298 K and 0.1 MPa. The equations of motion were integrated using the Verlet leapfrog algorithm scheme with a timestep of 2 fs. The cutoff radius for the Lennard–Jones interactions was fixed to 14 Å. A typical MD simulation consists of a first period of 10 ns followed by an additional production phase of 50 ns. The dimensions of the simulation cell were chosen in order to avoid self-entanglements. The validation of these atomistic force field was carried out by comparing the polymer melt density and its temperature dependence as well as the end-to-end distance. We have also checked some methodological aspects of the atomistic simulations (see Fig. S2 (ESI†)) establishing that the description of the polymer melts is as good as possible with an atomistic model. As an example, for the cis-1,4-PB polymer, the atomistic simulations give an end-to-end distance of 0.74 ± 0.03 Ų mol g⁻¹ for an experimental value of 0.76 Ų mol g⁻¹ with a ratio to the mean square end-to-end distance to the mean square radius of gyration equal to 6.1. These atomistic simulations¹⁸ have established a quantitative prediction of the thermal expansion coefficient of 6.5 × 10⁻⁴ K⁻¹ against an experimental value of 6.7 × 10⁻⁴ K⁻¹. Once this atomistic force field was validated, the bottom-up strategy used here consists of designing CG potentials from atomistic molecular dynamics simulations.

2.2 Dissipative particle dynamics (DPD)

The dissipative particle dynamics^49–53 method solves the Newtonian equations for particles subject to conservative, dissipative and random forces. Thus,where f_i is a pairwise-additive force defined as the sum of three contributionswhere f^C_ij, f^R_ij and f^D_ij are the conservative, random and dissipative forces, respectively. The conservative repulsive force f^C_ij derives from a soft interaction potential and is expressed aswhere w_c(r_ij) is the conservative potential, r_ij is the relative displacement of particles i and j and r̂_ij is the corresponding unit vector. r_c is the cutoff radius and is fixed to 4.0 nm. The dissipative and random forces are given by

f^D_ij = −γω^D(r_ij)(r̂_ij·v_ij)r̂_ij (4)

where δt is the time step. v_ij = v_i − v_j is the relative velocity, the terms ω^D(r_ij) and ω^R(r_ij) are dimensionless weighting functions. σ is the amplitude of the noise, θ_ij is a random Gaussian number with zero mean and unit variance. γ and σ are the dissipation strength and noise strength, respectively. Español and Warren⁵⁴ have shown that the system will sample the canonical ensemble and obey the fluctuation-dissipation theorem if the following conditions are satisfied.where k_B is Boltzmann's constant, T is the temperature and ω^D(r_ij) = (1 − r/r_c).²

The conservative potential w_c(r_ij) sums intramolecular and intermolecular interactions. The equations of motion were integrated using a modified version of the velocity-Verlet (DPD-VV) algorithm.⁵⁵ The development of integration schemes in DPD simulations is an area of active research.^8,56–59 We have compared the modified velocity-Verlet scheme⁵⁵ with that developed by Pagonabarraga et al.⁵⁶ using a self-consistent technique to calculate forces and velocities (DPD-SC). With the values of time step used here, we did not detect any significant deviations between the target and the calculated temperatures calculated by the DPD-VV and DPD-SC algorithms.

2.3 Coarse-graining procedure

From 10 independent atomistic configurations, we built bead–bead intermolecular and intramolecular potential functions calculated by coarse-graining atomistic configurations. This procedure allows a direct mapping between the atomistic and CG systems. The degree of coarse-graining was fixed to 5, indicating that each bead corresponds to 5 monomers. This degree of coarse-graining is then much smaller than the entanglement mass and allows the use of a relatively larger time step. By using a value of 5, no bond crossings are detected during the dynamics of the polymer chains. As a consequence, we expect a good reproduction of the entanglement effects. The time step was fixed to 50 fs and the temperature was kept to 300 K. Three different distributions functions are developed simultaneously: the radial distribution function g_bond(r_i,i+1) between two consecutive beads in the polymer chain; the radial distribution function g_bend(θ_i,i+1,i+2) between three consecutive beads in the polymer chain; the radial distribution function g_nb(r_i,j) between beads i and j of different polymer chains and beads i and j of the same molecule separated by more than one bond. From these radial distribution functions (RDF), we use the potential of mean force approach³⁵ to derive the corresponding CG potential functions as

w⁰_bond(r_i,i+1) = −k_BT ln g_bond(r_i,i+1) (7)

w⁰_bend(θ_i,i+1,i+2) = −k_BT ln g_bend(θ_i,i+1,i+2) (8)

w⁰_nb(r_i,j) = −k_BT ln g_nb(r_i,j) (9)

where w⁰_bond(r_i,i+1), w⁰_bend(θ_i,i+1,i+2), w⁰_nb(r_i,j) are the CG potentials corresponding to the bonding, bending and nonbonded interactions developed from the CG atomistic configurations. To improve the statistics, we have adopted the strategy³⁷ consisting of using different ways of grouping the five monomers in the atomistic configurations. The resulting potential functions have been tabulated: the force are then obtained by a cubic spline interpolation.⁶⁰

From these w⁰_bond(r_i,i+1), w⁰_bend(θ_i,i+1,i+2), w⁰_nb(r_i,j) coarse-grained potentials, we develop constant-DPD simulations of polymer melts formed by 40 polymer chains of 40 beads. A first DPD simulation of 10 000 steps was performed with the initial w⁰_bond(r_i,i+1), w⁰_bend(θ_i,i+1,i+2) and w⁰_nb(r_i,j) potentials. From these mesoscopic simulations, we re-calculate the corresponding radial distribution functions g_bondⁿ(r_i,i+1), g_bendⁿ(θ_i,i+1,i+2) and g_nbⁿ(r_i,j), respectively. We observe that the target RDFs in eqn (7), (8) and (9) are not accurately reproduced. By calculating the ratio between these intermediary RDFs and the target RDFs according to the Iterative Boltzmann Inversion process, we obtain new potential functions as

where gⁿ(r) and wⁿ(r) refer to the bonding, bending and non-bonded distribution and potential functions at the step n, respectively. g(r) is the target RDF. This scheme is reiterated until convergence is obtained.

An iteration consists of a constant-NVT DPD simulation of 10 000 steps and 5 iterations were required to obtain a good convergence. Additional constant-NPT DPD simulations were performed to modify these intermolecular potentials by using an attractive linear function¹⁸ in order to make the pressure in line with a pressure of 0.1 MPa. Fig. 2a shows that the intramolecular target RDFs obtained from atomistic simulations are very well reproduced by the CG models after the iterative process for the cis-1,4-PB, PIB and PDMS polymers. The resulting CG potential functions used in the DPD simulations are given in part b of Fig. 2. We observe some significant changes in the distribution functions between the different polymers indicating that the presence of the methyl groups and their resulting excluded-volume in the PDMS and PIB chains are well-reproduced by oscillations in the distribution functions in spite of a coarse-graining degree of five. Fig. 2c and e show the bonding and bending distribution functions and the corresponding CG potentials are represented in Fig. 2d and f, respectively. It is worth noting that the CG potential for the bond stretching resembles a simple harmonic potential that has been widely applied to model the interactions between consecutive beads^36,50,51,61 at the mesoscale. The deviation from this simple harmonic potential are polymer-microstructures dependent. From Fig. 2f, we observe that the CG bending potentials tend to disfavor the small angles in line with the bending potentials used for flexible and semi-flexible polymer chains.⁴⁶

Fig. 2 (a) Nonbonded radial distribution functions calculated from MD and DPD simulations and (b) nonbonded CG potential functions used in the DPD simulations obtained from the CG procedure. (c) Bonding and (e) bending CG distribution functions built from MD and DPD simulations. Corresponding (d) bonding and (f) bending CG potential functions obtained for a degree of coarse-graining of 5 monomers.

Since the dissipative forces are present in the simulations, another issues concern the choice of the friction coefficient^7,38,39 that has been shown to significantly impact on the dynamical properties of the polymer chains. It is well-known that the dynamical properties of polymer melts can be significantly affected by the choice of these adjustable parameters^8,38 r_c and γ. Based on the work³⁸ of Lahmar and Rousseau who have studied the influence of r_c and γ on the global and local dynamics of a polymer melt, we checked that γ = 300 kg mol⁻¹ ns⁻¹ was a good choice for the observation of long-time dynamical processes. Indeed, by using this value of γ, we have shown in a previous work¹⁸ that the diffusion coefficient of the cis-1,4-PB is slightly faster than that calculated from atomistic simulations but it remains on the same order of magnitude.

3 Discussions

3.1 Prediction of the thermomechanical and structural properties

It is now well-established that the combination of the CG models with DPD is able to equilibrate polymer melts containing high molecular mass chains^18,32 with a significant computational reduction (see Fig. S1 and S3 (ESI†)). Fig. 3a shows the density of the different polymer melts as a function of time. These trajectories resulting from stable constant-NPT simulations allow to extract time average densities that are given in Table 1 for each polymer. Interestingly, the CG potentials allow a very good reproduction of the melt density with a maximum deviation of 1% from experiments.⁶² Additionally, the change in the density due to the chemical nature of the monomers is also very well predicted by these mesoscale simulations. Fig. 3b shows the end-to-end distance of the polymer chains versus time. The different trajectories show that the CG models are able to distinguish the different polymers on this structural property even if the chemical details of the polymer are no longer explicitly described. The comparison with the experimental end-to-end distance shows a quantitative prediction with a maximum deviation of 7% for PIB (see Table 1). The ratio of the end-to-end distance to the radius of gyration, given for each polymer in Table 1 is in line with the expected value of 6 for ideal gaussian chains.⁶³ The bottom-up approach developed in the constant-NPT ensemble becomes predictive for the density and reproduces also accurately the end-to-end distance.

Fig. 3 (a) Simulated and experimental polymer melt densities as a function of time for cis-1,4-PB, PIB and PDMS polymers; (b) mean square end-to-end distances as a function of time. The experimental data, shown in dotted lines in (a) and (b) are taken from ref. 62.

Table 1 Entanglement molecular mass (M_e), melt density(ρ), plateau modulus (G⁰_N), end-to-end distance (R_ee/M), ratio to the mean square end-to-end distance to the mean square radius of gyration, thermal expansion coefficient α_p calculated from CG DPD simulations. The experimental values for each polymer are taken from ref. 62. The subscripts give the accuracy of the last decimal(s), i.e., 0.78₈ means 0.78 ± 0.08

	Me (g mol^−1)		ρ (g cm^−3)		(G^0N) (MPa)		Ree^2/M (Å^2 mol g^−1)		Ree^2/RG^2	αp (10^−4 K^−1)
	CG	Exp.	CG	Exp.	CG	Exp.	CG	Exp.		CG	Exp.
cis-1,4-PB	2500	2347	0.9109	0.900	0.731	0.76	0.788	0.76	6.1	6.51	6.7
PIB	6238	5686	0.9255	0.918	0.301	0.32	0.615	0.57	6.2	5.21	5.5
PDMS	10588	9613	0.9735	0.970	0.181	0.20	0.423	0.42	6.2	5.61	9.0

---

We now focus on the transferability of the CG potentials to the polymer chain length, temperature and pressure. Whereas the CG potentials have been developed from simulations of polymer chains of 200 monomers (40 beads), we use these potentials to simulate larger polymer chains of 800 monomers corresponds to a length of 160 beads. The polymer melt density and the end-to-end distances are given in Table S1 (ESI†). The deviations from experiments are not larger than those calculated with smaller chain lengths. This is an interesting result that shows the possibility of developing the CG potentials from shorter polymer chains. We now turn attention to the temperature transferability of the CG potentials by calculating the specific volume of the polymer melt in the [200–500] temperature range above the glass transition temperature. The simulated specific volumes result from constant-NPT simulations and their evolution with the temperature are represented in Fig. 4 at different temperatures. The linear temperature-dependence of the specific volume on the temperature allows to determine the thermal expansion coefficient defined by eqn (11). The value of the thermal expansion, given in Table 1, is calculated from the slope of the different lines at 300 K. We also report for comparison the experimental¹⁷ and simulated values of α_p at 300 K by dotted and solid lines in Fig. 4, respectively. We observe a good agreement between experiments and simulations for cis-1,4-PB and PIB polymers with a deviation of less of 5%. The deviation is more significant with PDMS leading to a deviation of 35%. We can conclude however that the order of magnitude for α_p is reproduced with the CG models developed in this work. A better agreement between experimental and calculated values of α_p for PDMS is possible but would require some refinements of the CG potential. The temperature dependence of the polymer melt density is then well-reproduced by the CG models in the temperature range above the glass transition temperature.

Fig. 4 Specific volume (1/ρ) of the three polymer melts as a function of temperature calculated using the CG models where ρ is the average polymer melt density. The solid and dotted lines represent to the simulated and experimental values of α_p at 300 K.

To complete the analysis, we investigate the pressure dependence of these CG potentials through the calculation of the isothermal compressibility κ_T from the volume fluctuations in the constant-NPT ensemble as

The CG models give values of κ_T at 300 K equal to 1090 × 10⁻¹⁵, 130 × 10⁻¹⁵, 120 × 10⁻¹⁵ bar⁻¹ for cis-1,4-PB, PIB and PDMS polymers, respectively. The corresponding experimental values¹⁷ for cis-1,4-PB, PIB and PDMS polymers are 7.2 × 10⁻¹⁵, 4.8 × 10⁻¹⁵, 11 × 10⁻¹⁵ bar⁻¹, respectively. This comparison establishes the inability of the CG potentials to reproduce any pressure dependence leading to deviations from experiments of two orders of magnitude as already observed for other CG models.^32,64 The procedure applied here allows to reproduce some thermodynamic properties at a target pressure but excludes any transferability to other pressures due to the deficiency of reproducing the incompressibility of the polymer material at equilibrium.

The explanation of this poor reproduction of the incompressibility comes from various origins: (1) the degree of coarse-graining, (2) the correction of the pressure during the development of the CG potentials. First, concerning the degree of coarse-graining, we develop CG potentials for a level of coarse-graining of 4 leading to harder potentials. The results are presented in Table S1 (ESI†). Fig. S4 (ESI†) shows that the trajectory of the simulated polymer density shows weaker fluctuations of the density than with λ = 5. This clearly establishes a relationship between the softness of the potential and the compressibility property. Second, the modification of the CG potentials to reproduce the pressure of the atomistic configurations introduces a pressure dependence in these CG potentials.^32,64

We propose here to analyze the CG configurations of the different polymer melts by applying the primitive path analysis (PPA⁶⁵) methodology, thoroughly described in ref. 18. Fig. 5 shows the distributions of the molecular weight (M_e) between entanglements. The average entanglement mass calculated from this distribution is reported in Table 1 along with the experimental corresponding property. We establish here a very good performance of the prediction with a maximum deviation of 10% from experiments. Additionally, the differences between the entanglement molecular weight of each polymer are also very well reproduced: the ratios M^PDMS_e/M^PIB_e and M^PDMS_e/M^cis-1,4-PB_e are equal to 1.7 and 4.2 for CG simulations against 1.7 and 4.1 for experiments, respectively. For PIB, an average entanglement mass of 6238 g mol⁻¹ corresponds to an entanglement length of 9 beads (45 monomers) and to the presence of approximately 17 entanglement segments per chain formed by 160 beads. As a result, we model here polymer chain lengths on the order of 20 entanglement lengths N_e whereas the longest chains that can be simulated by atomistic models within a reasonable computational cost represent roughly 2N_e. Since the plateau modulus G⁰_N is related to M_e through the relationship⁶² , it is then possible to estimate the plateau modulus G⁰_N of linear polymer melts. However, this plateau modulus could be calculated directly from CG simulations without making any assumption in any theoretical model but would require very long chains. The values of G⁰_N are reported in Table 1 and the agreement with experiments is quantitative with a maximum deviation of 10% for PDMS and 4% for cis-1,4-PB.

Fig. 5 Distributions of the entanglement molecular weight for cis-1,4-PB, PIB and PDMS polymers. The dotted curves represent the fits obtained from a Poisson distribution.

3.2 Application of an uniaxial deformation

We now turn attention to the transferability of the CG models on the direct calculation of the modulus G⁰_N of a crosslinked network of PIB. A crosslinked network allows to observe a plateau regime within a reasonable computational effort whereas the same calculation for polymer melts would require very large chains of polymers. We consider then a randomly tetrafunctional crosslinked PIB network where the crosslinking is ensured by an insertion of a few amounts of isoprene units which allows chemical bonding between chains through vulcanization. We have checked that the CG potential is not affected with a degree of coarse-graining of 5. The average molecular weight (M_c = 6379 g mol⁻¹) between crosslink units is approximately the entanglement molecular weight of the PIB melt leading to no trapped entanglement between crosslinks. In a first step, we aim to calculate the stress–strain curve of the network by applying a tension in the z-direction and by maintaining a constant pressure in the x and y-directions. More precisely, uniaxial tensile is applied by stretching the z-dimension of the simulation box under a negative normal pressure P_zz. We also apply a positive pressure P_xx = P_yy = 0.1 MPa in the x and y-directions. We used the anisotropic Berendsen barostat⁶⁶ for the control of the pressure. The magnitude of the elongation along the z-axis is measured throughorwhere L_z and L_z,0 are the box dimensions along the z-axis at times t and t = 0, respectively. This property is reported in Fig. 6a at three different values of stress tensor defined as

σ_zz = P^isotropic_zz − P_zz (15)

where P^isotropic_zz = 0.1 MPa was the value of the normal component in the system when the pressure was isotropically maintained. As P_zz varies from −0.2 to 0, σ_zz decreases from 0.3 to 0.1 MPa.

Fig. 6 Stress–strain curve of uniaxial deformation as a function of (λ − λ⁻²) over (a) a large range of deformations and (b) a limited range of deformations obtained from CG simulations, affine and phantom models. (c) Stress–strain curve as a function of ε in order to focus on the region of weak deformations and on the linear behavior (black dotted line) of the fitting curve (red curve).

Fig. 6a shows the stress–strain curve in a wide range of deformations from the elastic regime to the high elongation regime. We also plot the stress–strain curves deduced from the classical theories, the phantom network model⁶⁷ and the affine deformation model^68,69 of the rubber elasticity. We observe that the curve deviates from the phantom behavior for an elongation ratio of ε_z = 0.6 leading to an hardening of the network in the regime of strong deformations. We also observe that the stress–strain curve of the crosslinked PIB network exhibits a linear portion in the range of weak deformations and matches very well the phantom model. The agreement between the simulated and theoretical curves of the phantom model is not surprising since our network corresponds to a weakly constrained system due to the absence of additional topological constraints (no entanglement) between two crosslinks whereas the affine model describes a strongly constrained network. For strong deformations (Fig. 6a), the stress–strain curve deviates from both the affine and phantom theories as expected for the hardening of the network. In the linear regime (see Fig. S6 (ESI†)), the stress–strain curve allows to determine the Young's modulus E for uniaxial deformation where E is defined as the slope of where (ε = λ − 1) → 0. We find that E = 0.57 MPa. Since G = E/3 in the elastic regime for incompressible systems, i.e., with a Poisson's ratio of 0.5 (see Fig. 6b and S5 (ESI)†), we obtain a numerical value of the plateau modulus G⁰_N of 0.19 MPa in line with that calculated from MPa.

It should be noticed that the CG models are able here to reproduce the Poisson's coefficient characterizing incompressible systems whereas these same models were unable to reproduce the incompressibility of polymer melts (see Section 3.1). It means that adding topological constraints improves the description of the compressibility behavior. It seems that decreasing the interpenetration between the beads leads to a better reproduction of the compressibility in line with the impact of harder potentials on the density fluctuations of pure polymer melts. In the case of a crosslinked network of polymer, we can conclude that the CG models developed here perform very well in the reproduction of the Young's modulus and plateau modulus through the calculation of the stress–strain curve.

Another alternative for a direct calculation of the plateau modulus from equilibrium simulations is to consider the shear relaxation modulus G(t) from the autocorrelation of the stress tensor elements⁴⁴ by running very long equilibrium simulations. The shear relaxation modulus can be calculated from the autocorrelation of the stress tensor elements⁴⁴ as

with N_αβ = σ_αα − σ_ββ where α, β denote the x, y and z components. The components of the stress tensor are defined bywhere f^C,α_ij is the α component of the conservative force defined in eqn (3) and v^α_i the α component of the velocity of bead i.

Fig. 7 reports the calculated stress relaxation function of the crosslinked PIB network. As expected, the G(t) curve exhibits a plateau region. Very interestingly, the value of the plateau modulus fits very well with the values predicted by the nonequilibrium CG simulations and the phantom theory. We show here in the case a crosslinked network that the two different methodologies converge to the same value of the plateau modulus establishing that the combination of the CG models with the DPD method as an efficient tool for determining this property.

Fig. 7 Time evolution of the shear modulus G calculated from the stress relaxation function.

4 Conclusions

One of the challenge of the CG models in reducing the number of degrees of freedom is to incorporate the chemical nature of the molecules into the CG element. This is an essential step toward the development of realistic effective potentials. Additionally, the bottom-up approach addresses fundamental questions concerning the thermodynamic and structural consistencies of the CG models with the underlying atomistic configurations. This consistency is required for a quantitative prediction of thermomechanical, structural and viscoelastic properties of polymeric materials that are often only accessible at the mesoscopic scale.

We have developed CG models of different polymers (cis-1,4-PB, PIB and PDMS) by using the iterative Boltzmann inversion method. Nonbonded and intramolecular (bonding and bending) potentials have been built from distribution functions of the atomistic configurations. We took the route of developing the CG models in the constant-NPT ensemble in order to reproduce the pressure of the atomistic simulations. These models have been incorporated in the DPD method. We did not aim to focus on the dynamical properties of the polymer melt since we have shown recently that the dynamics of our CG models¹⁸ is slightly faster than that of the atomistic models.

Here, we have established that the (IBI + DPD) combined approach can be applied to different polymers and different properties. We have shown that the CG models are able to reproduce accurately the polymer melt density, the end-to-end distance, the entanglement mass of different polymer melts. From theoretical models, it is then possible to estimate the plateau modulus with a reasonable computational effort without running simulations of very long chains. The temperature and pressure transferabilities of these potentials have been discussed. It results a relatively good reproduction of the thermal expansion coefficient whereas the calculation of the isothermal compressibility shows significant deviations of two orders of magnitude from experiments. We have proposed some assumptions that could explain this behavior requiring however further investigations.

In order to perform a direct calculation of the plateau modulus with these CG models within a reasonable system-size, we have opted for a crosslinked network of PIB. We have shown that the two ways of calculating the plateau modulus either by applying a deformation either by calculating the autocorrelation of the stress tensor elements lead to the value expected by the classical theories of linear elasticity. This work underlines the performance of the bottom-up approach for quantitative predictions of thermomechanical properties of polymer with a reasonable computational cost.

References

1. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1989.
2. Monte Carlo and Molecular Dynamics Simulations in Polymer Science, ed. K. Binder, Oxford University Press, New York, 1995.
3. Simulation Methods for Polymers, ed. M. Kotelyanskii and D. N. Theodorou, Marcel Dekker, Inc., New York, 2005.
4. F. Müller-Plathe, ChemPhysChem, 2002, 3, 754–769.
5. V. A. Harmandis, N. P. Adhikari, N. F. A. van der Vegt and K. Kremer, Macromolecules, 2006, 39, 6708–6719.
6. C. Peter and K. Kremer, Soft Matter, 2009, 5, 4357–4366.
7. H. J. Qian, P. Carbone, X. Chen, H. A. Karimi-Varzaneh, C. C. Liew and F. Müller-Plathe, Macromolecules, 2008, 41, 9919–9929.
8. H. J. Qian, C. C. Liew and F. Müller-Plathe, Phys. Chem. Chem. Phys., 2009, 11, 1962–1969.
9. V. A. Harmandis and K. Kremer, Macromolecules, 2009, 42, 791–802.
10. D. Fritz, K. Koschke, V. A. Harmandaris, N. F. A. van der Vegt and K. Kremer, Phys. Chem. Chem. Phys., 2011, 13, 10412–10420.
11. Y. Li, S. Tang, B. C. Abberton, M. Kröger, C. Burkhart, B. Jiang, G. J. Papakonstantopoulos, M. Poldneff and W. K. Liu, Polymer, 2012, 53, 5935–5952.
12. Y. Li, B. C. Abberton, M. Kröger and W. K. Liu, Polymers, 2013, 5, 751–832.
13. R. Potestio, S. Fritsch, P. Espanol, R. Delgado-Buscalioni, K. Kremer and R. E. D. Donadio, Phys. Rev. Lett., 2013, 110, 10831.
14. K. Johnston and V. Harmandaris, Soft Matter, 2013, 9, 6696–6710.
15. K. Johnston and V. Harmandaris, Macromolecules, 2013, 46, 5741–5750.
16. S. Trément, B. Schnell, L. Petitjean, M. Couty and B. Rousseau, J. Chem. Phys., 2014, 140, 134113.
17. J. E. Mark, K. Ngai, W. Graessley, L. Mandelkern, E. Samulski, J. Koening and G. Wignall, Physical Properties of Polymers, Cambridge University Press, Cambridge, 3rd edn, 2004.
18. G. Maurel, B. Schnell, F. Goujon, M. Couty and P. Malfreyt, J. Chem. Theory Comput., 2012, 8, 4570–4579.
19. G. S. Grest and K. Kremer, Phys. Rev. A, 1986, 33, 3628–3631.
20. A. Baumgärter and K. Binder, J. Chem. Phys., 1981, 75, 2994–3005.
21. Computational Modelling of Polymers, ed. K. Binder and M. Dekker, New-York, 1992.
22. G. S. Grest, Curr. Opin. Colloid Interface Sci., 1997, 2, 271–277.
23. K. Kremer and K. Binder, Comput. Rep., 1988, 7, 259.
24. K. Kremer and G. S. Grest, J. Chem. Phys., 1990, 92, 5057–5086.
25. G. S. Grest, Adv. Polym. Sci., 1999, 138, 149.
26. T. Kreer, M. H. Muser, K. Binder and J. Klein, Langmuir, 2001, 17, 7804–7813.
27. J. T. Padding and W. Briels, J. Chem. Phys., 2001, 115, 2846–2859.
28. J. T. Padding and W. J. Briels, J. Chem. Phys., 2002, 117, 925–943.
29. R. D. Groot and K. L. Rabone, Biophys. J., 2001, 81, 725–736.
30. A. Ghoufi and P. Malfreyt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 051601.
31. A. Ghoufi and P. Malfreyt, J. Chem. Theory Comput., 2012, 8, 787–791.
32. T. Spyriouni, C. Tzoumanekas, D. Theodorou, F. Müller-Plathe and G. Milano, Macromolecules, 2007, 40, 3876–3885.
33. T. Mulder, V. A. Harmandaris, A. V. Lyulin, N. F. A. van der Vegt, K. Kremer and M. A. J. Michels, Macromolecules, 2009, 42, 384–391.
34. L. Lu, J. F. Dama and G. A. Voth, J. Chem. Phys., 2013, 139, 121906.
35. J. G. Kirkwood, J. Chem. Phys., 1935, 3, 300–313.
36. D. Reith, M. Pütz and F. Müller-Plathe, J. Comput. Chem., 2003, 24, 1624–1636.
37. X. Guerrault, B. Rousseau and J. Farago, J. Chem. Phys., 2004, 121, 6538–6546.
38. F. Lahmar and B. Rousseau, Polymer, 2007, 48, 3584.
39. F. Lahmar, C. Tzoumanekas, D. N. Theodorou and B. Rousseau, Macromolecules, 2009, 42, 7485–7494.
40. G. Milano and F. Müller-Plathe, J. Phys. Chem. B, 2005, 109, 18608–18619.
41. P. Carbone, F. Negri and F. Müller-Plathe, Macromolecules, 2007, 40, 7044–7055.
42. P. Carbone, H. Varzaneh and F. Müller-Plathe, J. Chem. Phys., 2008, 128, 064904.
43. O. Byutner and G. S. Smith, Macromolecules, 2002, 35, 3769–3771.
44. A. E. Likhtman, S. K. Sukumaran and J. Ramirez, Macromolecules, 2007, 40, 6748–6757.
45. W. B. Lee and K. Kremer, Macromolecules, 2009, 42, 6270–6276.
46. A. E. Likhtman and S. K. Sukumaran, Macromolecules, 2010, 43, 3980–3983.
47. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364.
48. Accelrys Software Inc., Material Studio Modeling Environment, Release 6.2, Accelrys Software Inc., San Diego, 2007.
49. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett., 1992, 19, 155–160.
50. P. Malfreyt and D. J. Tildesley, Langmuir, 2000, 16, 4732–4740.
51. F. Goujon, P. Malfreyt and D. J. Tildesley, Macromolecules, 2009, 42, 4310–4318.
52. C. Ibergay, P. Malfreyt and D. J. Tildesley, J. Chem. Theory Comput., 2009, 5, 3245–3259.
53. C. Ibergay, P. Malfreyt and D. J. Tildesley, J. Phys. Chem. B, 2010, 114, 7274–7285.
54. P. Español and P. B. Warren, Europhys. Lett., 1995, 30, 191.
55. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
56. I. Pagonabarraga, M. H. J. Hagen and D. Frenkel, Europhys. Lett., 1998, 42, 377–382.
57. G. Besold, I. Vattulainen, M. Karttunen and J. M. Polson, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62, R7611–R7614.
58. I. Vattulainen, M. Karttunen, G. Besold and J. M. Polson, J. Chem. Phys., 2002, 116, 3967–3979.
59. T. Soddemann, B. Dünweg and K. Kremer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 046702.
60. W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical Recipes: The art of scientific computing, Cambridge University Press, NY, New-York, 2007.
61. F. Goujon, A. Ghoufi, P. Malfreyt and D. J. Tildesley, Soft Matter, 2012, 8, 4635–4644.
62. L. J. Fetters, D. J. Lohse, D. Richter, T. A. Witten and A. Zirkel, Macromolecules, 1994, 27, 4639–4647.
63. M. Rubinstein and R. H. Colby, Polymer Physics, Oxford University Press, Oxford, 2003.
64. H. Wang, C. Junghans and K. Kremer, Eur. Phys. J. E, 2009, 28, 221–229.
65. S. K. Sukumaran, G. S. Grest, K. Kremer and R. Everaers, J. Polym. Sci., Polym. Phys. Ed., 2005, 43, 917–933.
66. H. J. C. Berendsen, J. P. M. Postma, W. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
67. H. M. James and E. Guth, J. Chem. Phys., 1947, 15, 669–684.
68. W. Kuhn, J. Polym. Sci., Polym. Phys. Ed., 1946, 1, 380–388.
69. L. R. G. Treolar, The Physics of Rubber Elasticity, Clarendon Press, Oxford, 1975.

---

Footnote

† Electronic supplementary information (ESI) available: The atomistic and CG models are compared on their ability of simulating well-equilibrated melts of long polymer chains. The results with a smaller degree of coarse-graining and different polymer chain lengths are given for comparison. The methodology corresponding to the simulation of the elongation of the crosslinked polymer network is described. See DOI: 10.1039/c4ra16417b

---