J. Phys. Chem. B 2005, 109, 18609-18619
18609
Mapping Atomistic Simulations to Mesoscopic Models: A Systematic Coarse-Graining Procedure for Vinyl Polymer Chains Giuseppe Milano*,† and Florian Mu1 ller-Plathe‡ School of Engineering and Science, International UniVersity Bremen, D-28759, Bremen, Germany ReceiVed: May 6, 2005; In Final Form: July 12, 2005
The paper introduces a systematic procedure to coarse-grain atomistic models of the largest family of synthetic polymers into a mesoscopic model that is able to keep detailed information about chain stereosequences. The mesoscopic model consists of sequences of superatoms centered on methylene carbons of two different types according to the kind of diad (m or r) they belong to. The corresponding force-field contains three different bonds, six angle and three nonbonded terms. Recently developed analytical potentials, based on sums of Gaussians for bond and angle terms of the mesoscale force field have been used. For the nonbonded part, numerical potentials optimized by pressure-corrected iterative Boltzmann inversion have been used. As test case we coarse-grained an atomistic all-atom model of atactic polystyrene. The proposed mesoscale model has been successfully tested against structural and dynamical properties for different chain lengths and opens the possibility of relaxing melts of high molecular weight vinyl polymers.
1. Introduction Vinyl polymers are polymers made from vinyl monomers; that is small molecules containing carbon double bonds (see Scheme 1). They make up the largest family of polymers. Polypropylene (RdCH3), poly(vinyl chloride) (RdCl), and polystyrene (RdC6H5), the most common commercial plastics, are examples of this kind of polymer chain. The production of these polyolefins is a large-scale industry with an estimated annual world production of about 100 million tons. The trend is rising, with a predicted growth of over 50% over the next 10 years. Polymeric materials are used for a legion of applications in a wide array of technological areas, and the extension of computational techniques to larger space and time scales is one of the main goals of current research in polymeric materials simulations. A systematic investigation of the structure-property relations for polymeric materials by computer simulations requires the preparation of equilibrated melts of long, entangled chains. This in principle can be achieved by sufficiently long molecular dynamics or Monte Carlo simulations. The range of molecular weight of atactic polystyrene for practical application depends mostly on the fabrication method and, of course, on the end-use application. Just to have an idea, the most common commercial grades of general-purpose polystyrene are three: easy flow (Mn ) 74.000; Mw) 218.000) and medium flow (Mn) 92.000; Mw 225.000), used for injection molding, and high heat polystyrene (Mn 130.000; Mw ) 300.000) used mostly for extrusion applications. There is market also for very low molecular weight (Mw < 10.000) polystyrenes, but it represents a small specialty market (plasticizers, paper coating, pressure-sensitive tapes). Many simulation studies, addressed to characterize bulk and/ or interface properties of this class of polymers, have been * Corresponding author. E-mail:
[email protected]. † Permanent address: Dipartimento di Chimica, Universita ` di Salerno, Via Salvador Allende, I-84081 Baronissi (SA), Italy. ‡ Present address: Eduard-Zintl-Institut fu ¨ r Anorganische und Physikalische Chemie, Technische Universita¨t Darmstadt, Petersenstrasse 20, D-64287 Darmstadt, Germany.
SCHEME 1
reported. Theodorou and co-workers performed Monte Carlo simulations of glassy atactic polypropylene melts1,2 and surfaces3 exposed to vacuum to elucidate microscopic aspects of structure and molecular organization. Atactic poly(vinyl chloride) conformational behavior in melt has been studied by molecular dynamics simulations.4,5 Many simulation studies about atactic as well as stereoregular polystytrene have been reported. Structure and ring dynamics,6 and, in comparison with twodimensional NMR spectroscopy, the local order of bulk amorphous polystyrene7 have been studied by atomistic simulations. Mattice and co-workers studied the orientation of phenyl rings and methylene bisectors at the free surface of atactic polystyrene.8 Molecular mobility in gels based on atactic polystyrene,9,10 calculation of positronium annihilation spectra of atactic polystyrene melts,11 and diffusion of small penetrants in syndiotactic polystyrene12 have been also reported. High molecular weight polymer chains are difficult to relax. The longest relaxation of an entangled polymer melt of length N scales at least as N3, giving at least N4 in CPU time, which is only feasible for relatively short chain lengths. Although computational power increases 10-fold every five years, the huge number of degrees of freedom limits fully atomistic approaches when it comes to investigating long polymer chains. One way to circumvent this problem is to reduce the degrees of freedom by coarsening the models and keeping only those degrees of freedom that are deemed relevant for the particular range of interest. Simple models for the study of meso- and macroscale phenomena in polymers have been used extensively.13,14 Some examples are the “time coarse graining”,15 lattice approaches such as the bond fluctuation model,16-18 and the high-coordination lattice.19-21 On a larger length scale,
10.1021/jp0523571 CCC: $30.25 © 2005 American Chemical Society Published on Web 09/08/2005
18610 J. Phys. Chem. B, Vol. 109, No. 39, 2005 dissipative particle dynamics (DPD) and smoothed particle dynamics (SPD) are frequently used to tackle hydrodynamic problems.22-25 The main drawback of these methods is connected with their generic nature. Most of them, in fact, do not distinguish between chemically different polymers. Recently, methods to automatically adjust force fields to specific systems on various length scales (atomistic and mesoscopic models) have been introduced, which produce highly accurate physical properties.26-29 The main point of this class of methods is the mapping of atomistic features to mesoscopic models. To this end, one groups several atom together into “super-atoms”. The potentials between them are adjusted to reproduce mainly structural properties of the polymer. The structure of an ensemble of polymer chains in an environment is most conveniently described by distribution of geometrical quantities, which can be intramolecular (distance between two adjacent super-atoms, angles between three subsequent superatoms, principal values of radius of gyration tensor, and so forth) or intermolecular (distance between superatoms belonging to different chains, distances between the centers of mass of different chains or chain fragments, and so on). All these parameters can be used as targets to be reproduced by the coarsegrained model. Following this approach, different polymer properties have been calculated from simulations. Coarse-grained force fields reproduce the radius of gyration (Rg) and the hydrodynamic radius (RH) of long polymer chains in aqueous solution in good agreement with experimental data.30 Viscosity for bisphenol A-polycarbonate (BPA-PC) at different temperatures has been calculated from simulation of coarse-grained models.31 Some of the parametrization work has been devoted to showing that atomistic simulations of polyisoprene (PI)32 can be mapped onto a simple bead-spring model incorporating excluded volume, connectivity, and stiffness.33,34 The general problem of creating a realistic model consists of writing down a Hamiltonian that can then be parametrized to reproduce properties of a system as well as possible. For this class of simulation techniques, an important point is the choice of suitable mapping schemes and the corresponding optimization strategy. In this paper we propose a systematic procedure to coarse grain atomistic models of vinyl polymers into a realistic mesoscopic model able to keep information about polymer stereosequences. As test case we coarse-grained an atomistic all-atom model of atactic polystyrene. To our knowledge, this is the first off-lattice coarse-grained model of a vinyl chain able to keep information about stereochemical composition and stereochemical sequences. As for lattice models, a procedure for the mapping of RIS models into a coarse-grained representation on a high coordination lattice modified so that it can accommodate vinyl polymers has been reported by Mattice and co-workers.21 2. Simulation of Oligomers at the Atomistic Level 2.1. Details of the Simulation. The molecular dynamics package YASP35 was used to run the all-atom simulations under constant pressure (Berendsen manostat with coupling time τp ) 2 ps) and constant temperature (Berendsen thermostat36 with τT ) 0.2 ps at a time step of 2 fs). All bond lengths were kept rigid by the SHAKE procedure. The cutoff for nonbonded interactions was rc ) 1.1 nm with a Verlet neighbor list37 cutoff of 1.2 nm. The polystyrene force field has been used already to describe different polystyrene based materials: polystyrene gels,9,10 amorphous PS calculation of positronium annihilation
Milano and Mu¨ller-Plathe spectra,11 anisotropy of diffusion of helium and CO2 in a nanoporous crystalline phase of syndiotactic PS (s-PS).12 Details about the force field and parameters can be found in refs 9 and 12. 2.2. Sample Preparation. Atactic random polystyrene chains of 10 monomers with five different stereochemistries were generated in vacuum at T ) 500 K. For all considered stereochemistry, the chain microstructure statistic is bernoullian with 46% of meso diads. A set of noninteracting 56 chains has been generated by a hybrid pivot Monte Carlo molecular dynamics approach method (PMC),38 as implemented in the GMQ software developed by Brown.39 In the PMC simulations, the Monte Carlo moves are random changes of randomly picked dihedral angles. This method has been proven to successfully generate meltlike conformations of the chains for a broad range of apolar and polar polymers, showing a good agreement with experimental mean radii of gyration (Rg) and end-to-end distances.40-44 It generates ensembles for a phantom chain, in which every atom interacts with neighboring atoms only if they are within a certain number of bonds along the chain. With this range correctly chosen, the ensemble corresponds to chains in a solvent or melt. The bond-based cutoff in the MD-PMC code, normally used to restrict the range of intramolecular nonbonded interaction, was set at 7. This value is large enough to avoid superposition between phenyl rings of different monomers of the same chain, but small enough to avoid chain collapse; further details about the algorithm and its applications can be found in refs 38 and 40-44. Subsequently, these chains were placed into a cubic periodic simulation cell at a reduced density of about 800 kg/m3. At this preparatory stage, no care was taken to avoid atom overlaps between different chains. The system was first energy-minimized to remove the worst close contacts. Energy minimization, however, still left some polymer entanglements, polymer chains poking through phenyl rings, and so on. All of this was efficiently removed by a 50 ps molecular dynamics run using soft-core potentials. Soft-core potentials are implemented in our simulation program YASP35 in the following way. The shortrange part (0-0.35 nm, e.g.) of the nonbonded potential energy function is replaced by a cubic spline. The spline coefficients are chosen to satisfy four conditions: the spline matches the value (i) and the derivative (ii) of the original potential energy function at the crossover distance; its derivative is zero at an interatomic distance r of zero (iii); most importantly, its value V0 at r ) 0 is finite (iv). This gets rid of the singularity at r ) 0 and allows atoms to pass through each other. We typically start these relaxation simulations with V0 ) 4kBT and increase V0 stepwise to 20 or 50 kBT, before turning to simulations with the full potential. These calculations were initially done at constant volume. After that the system was relaxed with the full nonbonded potential coupling the x, y, and z lengths of the simulation box separately to a pressure bath of 1 atm. This allowed the simulation box to distort slightly from cubic to orthorhombic. After 500 ps, the aspect ratio was frozen and the box volume as a whole was coupled to the isotropic pressure for the remaining simulations. The system was equilibrated for 2 ns, this simulation time was more than enough to reach the equilibrium density. A production run of 18 ns followed the equilibration. Values of density at 500 K, mean square end-to-end distance, and mean square radii of gyration result in agreement with available experimental data.47,48.
Mapping Atomistic Simulations to Mesoscopic Models
J. Phys. Chem. B, Vol. 109, No. 39, 2005 18611
Figure 2. Illustration of two possible mapping schemes for atactic polystyrene. (a) One bead corresponds to a diadic m or r unit. The center of these superatoms, as indicated by filled squares, are the methylene carbons. (b) One bead corresponds to an R or S asymmetric -CHR- group. Figure 1. Polystyrene m and r diads in transplanar conformation (hydrogen atoms on phenyl rings are omitted for clarity). The corresponding modified Fisher projections are also reported.
SCHEME 2
3. Mesoscopic Force Field Development 3.1. Technical Details. For CG simulations, the GMQ_num code,45 a version of the GMQ package39 able to handle numerical potentials was used. The bulk simulations were performed at constant temperature (500 K) and constant pressure (P ) 1 bar) for production runs. The time constants for the loose coupling thermostat and manostat were set to 0.15 and 5 ps. A time step of 15 fs was used. Nonbonded interactions were truncated beyond 15 Å. 3.2. Structural Features of Vinyl Chains. An asymmetric vinyl chain can be stereoregular or stereoirregular, due to the presence on alternate skeletal carbon atoms of asymmetric -CHR- groups. Given a direction of the main chain, it is possible to assign the absolute configuration (R or S) to each asymmetric carbon. The shortest distinguishing piece of a sequence is the diad. If two consecutive configurations have the same absolute configuration (e.g. R-R or S-S), the diad is designated as meso (m diad see Figure 1); if they are different (e.g. R-S or S-R), the diad is designated as racemo (r diad see Figure 1). If along the chain all the asymmetric centers have the same configuration, the chain is entirely made of m diads and is called isotactic. If the configuration is inverted along the chain at each carbon (RSRS.. or SRSR.., i.e., in the planar conformation and in the modified Fisher projection the substituents are on opposite sides), then the chain is entirely made of r diads and is called syndiotactic. If the configurations are randomly distributed the chain is called atactic. 3.3. Mapping. The basic idea of the proposed coarse-graining scheme is to consider the diads (m or r) as superatoms in the mesocale effective force field. As depicted in Figure 2, as the center of a superatom we choose the methylene carbon. The superatom will be of type m or r according to the diad to which the methylene carbon belongs. The force field corresponding to this choice of mesocale model will have two types of superatoms corresponding to m or r diads, three different bond types corresponding to the three undistinguishable triads (mm, mr, rr), six angle types corresponding to the six undistinguishable tetrads (mmm, mmr, mrm, mrr, rmr, rrr), as sketched, with the help of modified Fisher projections, in Scheme 2. Alternative mapping schemes, which also give information about chain stereosequences, could be envisaged. A possible choice could be two different superparticles corresponding to
the chirality of the monomers given by the asymmetric -CHRgroup. In this case, we would have two types of superatoms corresponding to R or S configuration, two different bond types corresponding to the two undistinguishable RR (or SS) and RS (or SR) situations, and three kinds of angles: RRR (or the equivalent SSS), RRS (or the equivalent SSR, SRR, RSS), RSR (or the equivalent SRS). To make the best choice among different mapping schemes, a good criterion is the probability distributions for bond lengths and bond angles evaluated for the superatoms based on these two alternative mapping schemes from the atomistic trajectory. Ideally the mapping should give sharp, single peaked distributions. The distributions for RR (or SS) and mm bond length, coming from the different choices of superatoms, are reported in Figures 3 and 4a, respectively. The mapping based on diads, for bond lengths, gives always single peaked and symmetric distributions for all three bonds. Furthermore, the choice based on diads makes a coarsegrained polymer model easier to connect to experimental data. In fact, polymer microstructure (i.e., the degree of tacticity) is characterized by 13C NMR spectroscopy in terms of diads, triads, tetrads, and pentad fractions. Thus, this approach allows to
18612 J. Phys. Chem. B, Vol. 109, No. 39, 2005
Milano and Mu¨ller-Plathe
Figure 3. Histogram of R-R (or SS) asymmetric -CHR- groups, bond length extracted from atomistic simulation of atactic polystyrene at 500 K according to the mapping scheme of Figure 2b.
connect directly the statistics coming from NMR spectra in a suitable coarse-grained model. A viable course of action for the optimization procedure of the mesoscale force field is to successively adjust the terms contributing to the total force field in the order of their relative strength:
Vstretching f Vbending f Vnonbonded 3.4. Atomistic Target Distributions. A realistic coarsegrained model of atactic polystyrene must produce the same distributions of three different bonds (mm, mr, rr), six angles (mmm, rrr, mmr, etc.), plus intermolecular and intramolecular RDFs extracted from the atomistic simulations. The distributions of the three possible bonds between two diads extracted from atomistic simulations at 500 K, are reported (solid lines) in Figure 4a, 4b, and 4c, respectively. The distributions of the six possible angles corresponding to atomistic simulations are reported (solid lines) in Figures 5a-5f. Is worth noting that, in the case of an angle, the distributions are weighted by a factor sin θ and renormalized:
P(θ) ) fnp(θ)/sin (θ)
(1)
where fn is the normalization factor. The complex shape of structural parameter distributions that one encounters in the simulation of polymer chains, such as those depicted in Figure 5, led to the choice of more flexible tabulated numerical potentials as alternatives to complicated analytical functions. In particular, the so-called iterative Boltzmann inversion has been applied to the problem of coarsegraining polymer models.26 The procedure uses the differences in the potentials of mean force between the distribution functions generated from a guessed potential and the target distribution functions to improve the effective potential successively. The numerical potentials achieve a perfect match of the corresponding RDFs or other structural distribution functions. However, they lack adjustable parameters with physical interpretation. Therefore, as explained in the next subsection, we used an analytical potential approach to handle bond and angle terms of a mesoscale force-field which is general and as flexible as the tabulated potentials.46 3.5. Bonded Potential. Mapping atomistic to mesoscopic models is an inverse problem: search for an interparticle
Figure 4. Histogram of (a) m-m, (b) m-r, (c) r-r bond length extracted from atomistic simulation of atactic polystyrene (empty circles) and from coarse-grained simulations (solid lines) at 500 K.
potential which reproduces a given distribution function used as target function. A multipeaked distribution of a structural parameter, say a bond angle θ, can be approximated by a sum of n Gaussian functions characterized by their centers (θci), total area (Ai), and width (wi): n
P(θ) )
∑ i)1
Ai
wix(π/2)
exp-2(θ-θci) /w i 2
2
(2)
Given a distribution P(θ) of some structural parameter θ such as a bond or an angle, a first approximation of the corresponding potential can be derived doing a simple Boltzmann inversion of P(θ).
Mapping Atomistic Simulations to Mesoscopic Models
J. Phys. Chem. B, Vol. 109, No. 39, 2005 18613
TABLE 1: Parameters of the Bond Potential Represented by Eq 6
TABLE 2: Parameters of the Bond Potential Represented by Eq 4
bond type
na
i
Ai
wi [Å]
lci [Å]
angle type
n
i
Ai
wi [deg]
θci [deg]
m-m m-r r-r
1 1 1
1 1 1
0.015 0.018 0.018
0.09 0.09 0.09
2.46 2.46 2.46
m-m-m
2
m-r-m
2
r-r-m
3
r-r-r
3
r-m-r
3
m-m-r
2
1 2 1 2 1 2 3 1 2 3 1 2 3 1 2
0.861 0.078 0.042 1.283 0.047 0.255 3.404 0.168 0.580 0.385 0.043 0.913 0.167 0.800 0.192
13.7 8.9 11.0 13.0 18.0 14.2 13.0 10.1 11.0 15.2 7.0 13.5 10.9 12.0 13.0
147.7 161.5 142.8 165.0 99.6 144.1 165.4 89.4 142.8 164.1 87.1 136.0 158.8 136.2 155.0
a
n is the number of Gaussians used for each force field term.
The corresponding potential obtained by Boltzmann inversion can be written as
Ai
n
V(θ) ) - kT ln
∑ i)1
exp-2(θ-θci) /w i 2
wix(π/2)
2
(3)
If we define gi(θ ) ) Ai/wix(π/2) exp-2(θ-θci)2/w2i, the potential and the corresponding expression for the force can be written in a more compact form: n
V(θ) ) - kT ln
gi(θ) ∑ i)1
(4)
(θ - θci)
n
gi(θ) ∑ i)1 F(θ) ) - 4kT
wi2 (5)
n
gi(θ) ∑ i)1
iterative Boltzmann inversion.26 The iterative Boltzmann inversion uses the differences in the potentials of mean force between the distribution functions generated from a trial potential and the true distribution functions to improve the effective potential successively. An effective nonbonded potential is derived V(r) from a given tabulated starting potential V0(r), targeted to match the radial distribution function g(r). Simulating a system with V0(r) will yield a corresponding g0(r) that is different from g(r). The potential therefore needs to be improved, which is done by a correction term kBT ln[(g0(r)/g(r)]. This step can be iterated as follows by:
RDFj(r) Vj+1(r) ) Vj(r) + kT ln RDFtarget(r)
For the bond potentials it is possible to derive similar equations for potential (eq 5) and forces (eq 6): n
V(l) ) - kT ln
gi(l) ∑ i)1
F(l) ) - 4kT
l
gi(l) ∑ i)1 n
until
ftarget )
(l - lci)
n
1
(6)
wi2 (7)
gi(l) ∑ i)1
These kind of potentials allow an easy physical interpretation of the parameters. If we consider the simplest case with one Gaussian function, the physical meaning of the potential parameters becomes clear. If n ) 1, then eq 3 becomes 2 2 A exp-2(θ-θc) /w ) V(θ) ) - kT ln wx(π/2) 2kT (θ - θc)2 + const (8) 2 w
Equation 8 is the potential of an harmonic oscillator centered at θc with 4kT/w2 as spring constant. Additional information about the multicentered Gaussianbased potentials (MG-potentials) employed here can be found in ref 46. In Tables 1 and 2, parameter values for bond and angle types are reported. 3.6. Nonbonded Potential and Pressure Optimization. As for the nonbonded part of the potential, we use the pressurecorrected CG numerical potential, which is optimized by
(9)
∫0cutoff ω(r)[RDFj(r) - RDFtarget(r)]2 dr
(10)
falls below an initially specified threshold. We apply ω(r) ) exp(-r) as a weighting function in order to specifically penalize deviations at small distances. Pressure corrections can be introduced in polymer bulk systems by adding a weak linear potential term ∆V to the attractive long-range part of Vj(r), as in most cases the initially optimized structure yields a pressure different from the one at which the atomistic simulations are performed. We then postoptimize the structure of the system according to the above scheme in turn with linear additions until also the pressure matches the atomistic system. The so-called ramp correction is of the form
(
∆V(r) ) V0 1 -
r rcutoff
)
(11)
It vanishes at the cutoff (∆V(rcutoff) ) 0) and has V0 ) ∆V (r ) 0) as the only parameter. Depending on the pressure being above or below the target value, V0 is negative or positive. Taking this as a starting condition, we reoptimize the system against the structure. This combined method was iterated until the structure and the pressure were in agreement with target values. The pressure correction mainly acts on the longer range attraction, whereas the structure is dominated by the short-range repulsion. The calculated density obtained with force field optimized on oligomeric systems generally increases with the chain length. This is due to the fact that end monomers
18614 J. Phys. Chem. B, Vol. 109, No. 39, 2005
Milano and Mu¨ller-Plathe
Figure 5. Histogram of (a) m-m-m, (b) m-m-r, (c) r-m-r, (d) r-r-r, (e) m-r-m, (f) r-r-m angle extracted from atomistic (empty circles) and from coarse-grained simulations (solid lines) at 500 K.
contribute more to the free volume than to the inner monomers, and in a system with a given number of particles, the proportion of end monomers decreases when the chain length increases. The force field is required to be applicable to a broad range of molecular weight, and therefore it should be optimized differently for the end beads that for inner beads. The ramp correction to the nonbonded potentials has been tested on seven different coarse-grained systems, with different chain lengths. In the system with 350 diads, 99% of the nonbonded interactions are of the type mm, mr, or rr; therefore, the influence of all other interactions involving end monomers is relatively small. The mm, mr, and rr potentials were ramped to have the correct pressure of ∼1 bar. Then, to have the correct pressure for low molecular weight as well, the em and er potentials were modified iteratively with one ramp, and the ee (e is the end group diad) potential with a second ramp of the
opposite sign. This procedure finally gives the values of density ranging from 0.94 g/cm3 for the system with N ) 9 to 0.97 g/cm3 for N ) 350. Accurate experimental density measurements of monodisperse polystyrene samples (Mn ) 51.000; Mw/Mn < 1.06) above the Tg in the range 373-595 K have been done.47 In this range, volume changes at small intervals have been measured and the data were found to fit a second-order polynomial. From the coefficients reported in ref 47, an experimental density of 0.95 g/cm3 at 500 K is obtained. In Figure 6, the intermolecular pair correlation functions coming from atomistic simulations (solid line) of a-PS are compared with those obtained with the coarse-grained models (empty circles). In all cases the distributions obtained at atomistic level and by the coarse-grained model are in good agreement. The pressure corrected potentials for mm, mr, and rr interactions are reported in Figure 7.
Mapping Atomistic Simulations to Mesoscopic Models
J. Phys. Chem. B, Vol. 109, No. 39, 2005 18615
Figure 7. Coarse-grained nonbonded numerical potentials for (a) m-m, (b) m-r, (c) r-r interaction.
Figure 8. Values of gyration radii 〈S2〉1/2 obtained by simulations of polystyrene coarse-grained models (9), as function of molecular weight, in comparison with available experimental data (∆ SANS data from melts, O light scattering data from θ solutions).
Figure 6. Interchain (a) m-m, (b) m-r, (c) r-r pair correlation function from atomistic simulations (empty circles) and obtained with the coarsegrained model (solid lines) at 500 K.
4. Simulations of Longer Coarse-Grained Chains Once we have, at a coarse-grained level, a model able to reproduce the target distributions, we are able to access larger time and length scales by MD simulations. In this section, structural and dynamical quantities calculated from coarsegrained simulation trajectories will be compared with available experimental data. 4.1. Structural Parameters. Values of gyration radii 〈S2〉1/2 obtained by coarse-grained simulations of polystyrene melts at T ) 500 K with different chain lengths, as function of molecular weight, in comparison with available experimental data are reported in Figure 8. Unfortunately, small-angle neutron scattering (SANS) data for polystyrene melts are very limited, and there are no data for melts with chains shorter than Mw ) 10000.48 However, it is usually expected that chain statistics in
θ solution are identical to the chain statistics in the melt. This assumption has been confirmed by many SANS measurements and simulations.49-52 From the figure is clear that a good agreement between coarse-grained simulations and experimental measurements in the entire range of molecular weight has been found. From the calculated values of the mean square end-to-end distance R2, it is possible to extract the value of Flory’s characteristic ratio Cn for different chain lengths:
Cn ) 〈R2〉/nl2
(12)
In Figure 9 the values of Cn for different chain lengths are reported. Starting from values close to 200, the Cn starts to saturate at C∞. The value of C∞ ∼ 8 we found is in good agreement with experimental value of 8.3 at 500 K.53 As reported in a recent paper,54 the values of Cn can be reported as function of 1/n (Figure 9b). This allows to check if the data, as expected from the theory for large values of n, adhere to a straight line. In the present case, the curve reported in Figure 9b shows a positive curvature. Reproduction of structural properties of such as end-to-end distance and gyration radius of a-PS in the entire range of molecular weight has been achieved without inclusion of any effective torsional potential. The reason can be understood considering the distributions of torsional angles extracted from atomistic reference simulations.
18616 J. Phys. Chem. B, Vol. 109, No. 39, 2005
Figure 9. Flory’s characteristic ratio Cn (defined by eq 12) calculated from coarse-grained simulations as function of chain length. The straight line at ∼8 corresponds to the experimental value of C∞.
Figure 10. Correlation of DiDi+1 vectors, separating the ith and the (i+1)th superatom of the chain, as function of k (k is the number of bonds separating the two vectors to be correlated).
Milano and Mu¨ller-Plathe coarse-grained simulations. Furthermore, in the framework of hindered rotation model, we evaluated the Cφ ratio ((1 + 〈cosφ〉)/(1 - 〈cosφ〉)), where 〈cos φ〉 ) (∫2π 0 cosφ p(φ) dφ)/ (∫2π 0 p(φ) dφ) and p(φ) are our torsion distributions. This gives an idea on the chain dimensions effects due to different distributions p(φ). The values obtained for 〈cos φ〉 and -0.090 for the distribution obtained from the reference atomistic simulations averaging on all the possible quadruplets is -0.12 (Cφ ) 0.79), and the value for the distribution obtained by coarse-grained simulations without inclusion of torsional terms is -0.09 (Cφ ) 0.83). In both cases Cφ is ∼ 0.8. In the case of a partially full stereoregular polymer, the situation would be different because it is dominated by only one of the possible quadruplets. The value of 〈cos φ〉 calculated from the distribution of the stereoregular rrrr sequence extracted from atomistic simulation is 0.21 (with Cφ) 1.5). This has to be compared with -0.09 (Cφ ) 0.83) obtained with coarsegrained simulations without torsions. In this case, a nonnegligible effect of the torsion term leading to an enlargement of the chain dimension is expected. With respect to this point, simulations of full syndiotactic chains (chains containing only r superatoms) have been performed. For the longest chain (N ) 350), the calculated values of characteristic ratio for a syndiotactic chain, without inclusion of torsional potential, is ∼5. In disagreement with experimental findings, the Cn becomes even lower than the value of ∼8 obtained for the corresponding atactic chain. Preliminary simulations of a syndiotactic chain have been performed, including a numerical torsional potential for the rrrr sequence obtained by iterative Boltzmann inversion of the distribution reported in Figure 11. In this case, in good agreement with experimental findings,56 the calculated characteristic ratio for N ) 350 diads is ∼9.57 However, although we have this agreement, we want to stress that, in principle, the model used in our preliminary simulations has been derived from stereoregular sequences in an atactic polymer. A correct coarse-grained model for pure syndiotactic polystyrene will probably be different. A simple unified description of all ideal chains is provided by the equivalent freely jointed chain. The equivalent chain has the same mean-square end-to-end distance 〈R2〉 and the same maximum end-to-end distance as the actual polymer, but has N freely jointed effective bonds of length b. The effective bond length is called Kuhn length. The Kuhn length can be calculated from the value of C∞ and Rmax in the following way:
b)
Figure 11. Histogram of dihedral angle distribution (O) averaged on all possible quadruplets, (0) r-r-r-r extracted from atomistic simulations and (4) from coarse grained simulations without torsional potential at 500 K. According to the proposed mapping scheme, the r-r-r-r distribution extracted from atomistic simulations has been obtained averaging over both R-S-R-S-R-S and S-R-S-R-S-R sequences. This way of averaging leads to symmetric torsions in the mesoscale model of the syndiotactic polymer.
In Figure 11 the distribution of torsional angles averaged on all possible quadruplets (empty circles) for stereoregular rrrr pentads (empty squares) and for the coarse-grain simulation are plotted. The distribution obtained by averaging an all possible quadruplets is flatter than the one obtained considering stereoregular sequences and close to the distribution obtained from
C∞nl2 Rmax
(13)
where Rmax is the largest possible end-to-end distance and is determined by the product of the number of skeletal bonds n and their projected lengths:
Rmax ) nl cos
θ 2
(14)
For N ) 200, the calculated value of b is 21 Å. The corresponding persistence length Lp (half of b) is 10.5 Å. The persistence length can be also evaluated directly a simulation trajectory from the correlation of DD vectors:
CDD(k) ) 〈(DiDi+1‚Di+kDi+k+1)/|DiDi+1|Di+kDi+k+1|〉 (15) where DiDi+1 is the vector separating the ith and the (i+1)th diad of the chain. The averaging is performed over all positions i along the chain, over all chains, and over all records of the
Mapping Atomistic Simulations to Mesoscopic Models
J. Phys. Chem. B, Vol. 109, No. 39, 2005 18617
trajectory file. The variable k is the number of bonds separating the two vectors to be correlated. In Figure 10, the CDD(k) autocorrelation function of a chain of 200 superatoms is plotted. The autocorrelation function CDD(k) then can be fitted using the following exponential form in order to calculate the persistence length:
CDD(k) ≈ C0 exp(- k/Lp)
(16)
The best least-squares fit of the data of Figure 10 to eq 12 gives a value of Lp of 3.86 bonds. Multiplying Lp by the average bond length (2.56 Å), we obtain a persistence length of 9.88 Å with a corresponding Kuhn length b of 19.8 Å. Both values obtained from Cn and from the best fit of the curve of Figure 10 give a value in good agreement with experimental one of ∼10 Å.50 4.2. Dynamical Properties. 4.3. Chain Relaxation. To characterize the efficiency of the coarse-grained approach in the global relaxation of the chains, the time autocorrelation functions for the end-to-end vector, R, and the square end-to-end distance, R2 ) R‚R, were calculated:
C1(t) ) 〈R(0)‚R(t)〉/〈R2〉
(17)
C2(t) ) 〈R2(0)‚R2(t) - 〈R2〉2〉/(〈R4〉 - 〈R2〉2)
(18)
For the coarse-grained model with N ) 9, C1(t) and C2(t) decayed from one to zero in ∼50 and ∼25 ps, respectively, whereas the value for the atomistic system never went lower than 0.87 and 0.55 in 20000 ps. This means that the coarsegrained model is much more efficient to relax the end-to-end distance than the atomistic model. To obtain a relaxation time characterizing the sampling of conformational space, C1(t) and C2(t) were fitted to the stretched exponential (Kohlrausch-Williams-Watts) functional form:
C(t) ) exp(- (t/R)β)
(19)
The relaxation time equals the time integral of the stretched exponential, which is analytic and can be expressed with the Euler Γ function:
τ)
∫0∞ exp(- (Rt ) ) dt ) RβΓ(β1) β
(20)
Values of relaxation times for different chain lengths range from an order of magnitude of picoseconds for the shortest chains (for N ) 9 τ1 ∼ 10 and τ2 ∼ 2 ps, for N ) 30 τ1 ∼ 200 and τ2 ∼ 40 ps) to an order of magnitude of nanoseconds for the longest ones (for N ) 200 τ1 ∼ 17 and τ2 ∼ 3 ns, for N ) 350 τ1 ∼ 100 and τ2 ∼ 50 ns). The values of τ1 and τ2 for the reference atomistic system with N ) 9 cannot be properly evaluated because the corresponding C1(t) and C2(t) correlation functions do not decrease significantly during the simulation and the fitting is not reliable. To evaluate the reliability of brute force MD equilibration of polystyrene melts with the proposed coarse-grained model, the behavior of chain relaxation times with chain length can be considered. In particular, the values we obtained with different chain lengths from N ) 9 to N ) 350 have been fitted by the an exponential function in the form of eq 21
τ ≈ c1eN/ν - c0
(21)
Fixing a maximum simulation time (tmax in eq 22), from the values of c1, ν, and c0 we can estimate that the maximum chain
length of a melt can be relaxed by brute force MD
Nmax ≈ ν ln
(
)
tmax + c0 c1
(22)
From the best fit of eq 21, for τ1 we got c1 ) 2.2 ns, ν ) 92, and c0 ) 2.7 ns; for τ2 we got c1 ) 0.1 ns, ν ) 55, and c0 ) 0.11 ns. A reasonable choice for tmax is ∼10.000 ns.58 In both cases (for τ1 and τ2), we estimated a maximum chain length for brute force MD equilibrations of melts of N∼500 diads, corresponding to a molecular weight Mw ∼ 80.000. The application of the procedure recently proposed by Everaers and Kremer,61 a combination of MD relaxation, double bridging, and slow push-off, can lead to an efficient and controlled preparation of equilibrated melt up to N ) 7000 coarse-grained beads. For our coarse-grained model of polystyrene, this corresponds to a Mw of ∼700.000. This value is well beyond the highest molecular weights used for extrusion applications (Mn 130.000; Mw ) 300.000). 4.4. Diffusion Coefficients. Tracer diffusion coefficients of the chains can be calculated from the mean-squared displacements ∆r2 of the centers of mass via the Einstein relation.
D≈
1d 〈∆r2〉 6 dt
(24)
The calculated diffusion coefficients at 500 K range from 4.8×10-4 cm2/s for the shortest chain length (N ) 9) to 4.6×10-11 cm2/s for the longest one (N ) 350). The techniques to map the structural static properties of polymers which have been described in the paper lead inherently to larger time scales as the fastest inherent degrees of freedom are now motions of superatoms of the size of diads. If quantitative information about dynamical properties is desired, one has to find a correct mapping of the time scales of the atomistic simulation to the mesoscale. An obvious candidate for calibrating the time scale is the chain diffusion coefficient obtained by reference atomistic simulations. The time scale of the coarse-grained model can be scaled to bring the diffusion coefficient to the same value. Following this approach, the time scaling factor we obtain is ∼200. Considering the time coarse-graining for our mesoscale models, the calculated diffusion coefficients will be 2.4 × 10-6 cm2/s for N ) 9 and 0.23 × 10-12 cm2/s for N ) 350 this. The value obtained by from viscosity data on atactic polystyrene decamers (i.e. N ) 9 diads) using the Stokes-Einstein relation is 5 × 10-7 cm2/s at 500 K.62 As for larger molecular weights, the self-diffusion coefficients of monodisperse polystyrenes, which vary from 33.000 to 943.000, have been evaluated at 447 K by marker displacement measurements.63 The authors report a fit of the data D ) 0.007 M-2(0.1 cm2/s. From this equation the estimated value of D for N ) 350 diads is between 1.89 × 10-12 and 1.54 × 10-11 cm2/s. Figure 12 shows the simulation results for the dependence of the diffusion coefficients on the chain length. The data are plotted as DN as function of chain length N. For short chain lengths, ranging from N ) 9 to N ) 70, the calculated DN align on a straight line. In agreement with the Rouse ,model the diffusion coefficients scale as N-1. From the picture is clear that when the chain becomes longer (from N > 100), the slope changes to negative value of approximately -1. This is in agreement with a scaling behavior of the diffusion coefficient of N-2. The behavior of the diffusion coefficient found in the MD simulations agrees with the prediction of the reptation theory. The number of points in the plot is not enough and the
18618 J. Phys. Chem. B, Vol. 109, No. 39, 2005
Milano and Mu¨ller-Plathe Foundation for a fellowship. F.M.-P. thanks the Fonds der Chemischen Industrie for financial support. This work was partially supported by the MURST of Italy (grants PRIN 2002 and FISR 1999), Universita` di Salerno (Fondo Medie Apparecchiature 2003). We thank CINECA for allowing us CPU time (Progetti di Supercalcolo convenzione CINECA/INSTM). We thank Prof. Gaetano Guerra and Dr. Guglielmo Monaco of University of Salerno for useful discussions. References and Notes
Figure 12. Dependence of the diffusion coefficients on the chain length. The data are plotted as DN as a function of chain length N. For short chain lengths, the calculated DN align on a straight line. In agreement with the Rouse model, the diffusion coefficients scales as N-1. When the chain becomes longer, the slope changes to negative value of approximately -1. This, in agreement with the prediction of the reptation theory, shows a scaling behavior of the diffusion coefficient of N-2.
entanglement length can be evaluated only roughly around Ne ) 100 diads, which corresponds to Me ∼ 11.000, while the experimental estimations give Me ∼ 13.000, corresponding to Ne ) 120 diads. Still, these results are very encouraging, considering that the force-field has been optimized only against structural static properties such as atomistic RDFs and geometric distributions. 5. Conclusions A systematic procedure has been introduced to coarse grain atomistic models of vinyl polymers, the largest family of synthetic polymers, into a mesoscopic model that is able to keep detailed information about polymer chain stereosequences. The mesoscopic model consists of sequences of superatoms centered on methylene carbons of two different types according to the kind of diad (m or r) to which they belong. No extra superatoms for the side groups need be introduced to carry information about a given stereosequence. As a test case, we coarse-grained an atomistic all-atom model of atactic polystyrene. The change to a coarse-grained scale led to an effective speed-up of 2000 for the computational efficiency of the relaxation of the chains. The proposed mesoscale model has been successfully tested against structural and dynamical properties of the melt for different chain lengths and opens the possibility of relaxing melts of high molecular weight vinyl polymers. For the considered polymer, atactic polystyrene, melts up to a Mw of ∼700.000 can be prepared. This value is well beyond the highest molecular weights used for extrusion applications (Mn 130.000; Mw ) 300.000). A further relevant use of the proposed coarse-grained model is the generation of well-equilibrated atomistic amorphous phases of vinyl polymers of high molecular weight. According to the proposed mapping scheme, the diad distribution can be taken in to account in a possible back-mapping strategy. In the cases where the chirality of one chain end group has been fixed at random as R or S, the diad distributions can be translated in terms of absolute configurations of the CHR groups and the approach indicated by Tscho¨p et al.64 can be readily adapted to the present mapping scheme. Acknowledgment. We are indebted to Prof. David Brown and Dr. Severine Queyroy of Universite´ de Savoie for providing the GMQ_num code. G.M. expresses gratitude to the Humboldt
(1) Theodorou, D. N.; Suter, U. W. Macromolecules 1985, 18, 1467. (2) Dodd, L. R.; Boone, T. D.; Theodorou, D. N. Mol. Phys. 1993, 78, 961. (3) Mansfield, K. F.; Theodorou, D. N. Macromolecules 1990, 23, 4430. (4) Smith, G. D.; Jaffe, R. L.; Yoon, D. Y. Macromolecules 1993, 26, 298. (5) Tanaka, G.; Mattice, W. L. Macromolecules 1995, 28, 1049. (6) Rapold, R. F.; Suter, U. W.; Theodorou, D. N. Macromol. Theory Simul. 1994, 3, 19. (7) Robyr, P.; Tomaselli, M.; Grob-Pisano, C.; Meier, B. H.; Ernst, R. R.; Suter, U. W. Macromolecules 1995, 28, 5320. (8) Clancy, T. C.; Jang, J. H.; Dhinojwala, A.; Mattice, W. L. J. Phys. Chem. B 2001, 105, 11493. (9) Mu¨ller-Plathe, F. Macromolecules 1996, 29, 4782. (10) Witt, R.; Sturz, L.; Dolle, A.; Mu¨ller-Plathe, F. J. Phys. Chem. A 2000, 104, 5716. (11) Schmitz, H.; Mu¨ller-Plathe, F. J. Chem. Phys. 2000, 112, 1040. (12) Milano, G.; Guerra, G.; Mu¨ller-Plathe, F. Chem. Mater. 2002, 2977. (13) De Gennes, P.-G. Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca, New York, 1979. (14) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press: New York, 1986. (15) Forrest, B.; Suter, U. J. Chem. Phys. 1995, 102, 7256. (16) Carmesin, I.; Kremer, K. Macromolecules 1988, 21, 2819. (17) Paul, W.; Binder, K.; Kremer, K.; Heermann, D. Macromolecules 1991, 24, 6332. (18) Baschnagel, J.; Binder, K.; Paul, W.; Laso, M.; Suter, U.; Batoulis, I.; Jilge, W.; Bu¨rger, T. J. Chem. Phys. 1991, 95, 6014. (19) Rapold, R.; Mattice, W. J. Chem. Soc., Faraday Trans. 1995, 91, 2435. (20) Cho, J.; Mattice, W. Macromolecules 1997, 30, 637. (21) Halilogˇlu, T.; Mattice, W. J. Chem. Phys. 1998, 16, 6989. (22) Hoogerbrugge, P.; Koelman, J. Europhys. Lett. 1992, 19, 155. (23) Schlijper, A.; Hoogerbrugge, P.; Manke, C. J. Rheol. 1995, 39, 567. (24) Groot, R.; Warren, P. J. Chem. Phys. 1997, 107, 4423. (25) Espan˜ol, P.; Serrano, M.; Zuniga, I. J. Mod. Phys. C 1997, 8, 899. (26) Reith, D.; Pu¨tz, M.; Mu¨ller-Plathe, F. J. Comput. Chem. 2003, 24, 1624. (27) Meyer, H.; Biermann, O.; Faller, R.; Reith, D.; Mu¨ller-Plathe, F. J. Chem. Phys. 2000, 113, 6264. (28) de Vries, A. H.; Mark, A. E.; Marrink, S. J. J. Phys. Chem. B 2004, 108, 2454. (29) For a recent review see: Mu¨ller-Plathe, F. ChemPhysChem 2002, 3, 754. (30) Reith, D.; Meyer, H.; Mu¨ller-Plathe, F. Macromolecules 2001, 34, 2335. (31) Tscho¨p, W.; Kremer, K.; Batoulis, J.; Bu¨rger, T.; Han, O. Acta Polym. 1998, 49, 61. (32) Faller, R.; Mu¨ller-Plathe, F.; Doxastakis, M.; Theodorou D. Macromolecules 2001, 34, 1436. (33) Faller, R.; Kolb, A.; Mu¨ller-Plathe, F. Phys. Chem. Chem. Phys. 1999, 1, 2071. (34) Faller, R.; Mu¨ller-Plathe, F.; Heuer, F Macromolecules 2000, 33, 6602. (35) Mu¨ller-Plathe, F. Comput. Phys. Commun. 1993, 78, 77. (36) Berendsen, H. J. C.; Postma, J.; van Gunsteren, W.; Di Nola, A.; Haak, J. J. Chem. Phys. 1984, 81, 3684. (37) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (38) Neyertz, S.; Brown, D. J. Chem. Phys. 2001, 115, 708. (39) Brown, D. The gmq User Manual Version 3, http://www.univsavoie.fr/labos/lmops/brown/gmq.html. (40) Brown, D.; Clarke, J. H. R.; Okuda, M.; Yamazaki, T. J. Chem. Phys. 1994, 100, 1684. (41) Brown, D.; Clarke, J. H. R.; Okuda, M.; Yamazaki, T. J. Chem. Phys. 1994, 100, 6011. (42) Brown, D.; Clarke, J. H. R.; Okuda, M.; Yamazaki, T. J. Chem. Phys. 1996, 104, 2078.
Mapping Atomistic Simulations to Mesoscopic Models (43) Neyertz, S.; Brown, D.; Clarke, J. H. R. J. Chem. Phys. 1996, 105, 2076. (44) Neyertz, S.; Brown, D. J. Chem. Phys. 1995, 102, 9725. (45) Queyroy S. Simulations moleculaires dynamiques de surfaces de polymeres amorphes: cas de la cellulose; Ph.D. Thesis; Universite de Savoie, 2004. (46) Milano, G.; Goudeau, S.; Mu¨ller-Plathe, F. J. Polym. Sci., Part B: Polym. Phys. 2005, 43, 871. (47) Hocker, H.; Blake, G. J.; Flory, P. J. Faraday Soc. Trans. 1971, 67, 2251. (48) Ding, Y.; Kisliuk, A.; Sokolov, A. P. Macromolecules 2004, 37, 161. (49) Wignall, G. D.; Ballard, D. G. H.; Schelten, J. Eur. Polym. J. 1974, 10, 861. (50) Cotton, J. P.; Decker, D.; Benoit, H.; Farnoux, B.; Higgins, J.; Jannink, G.; Ober, R.; Picot, C.; des Cloizeaux, J. Macromolecules 1974, 7, 863. (51) Ballard, D.; Rayner, S.; Shelten, J. Polymer 1976, 17, 349. (52) Yoon, D. Y.; Flory, P. J. Macromolecules 1976, 9, 294. (53) At 500 K a value of C∞ of 8.3 can be obtained from a value of d ln C∞/d T -0.90 × 103 calculated by Flory and Yoon in ref 53. (54) Mattice, W. L.; Helfer, C. A.; Sokolov, A. P. Macromolecules 2004, 37, 4711. (55) Yoon, D. Y.; Sundararajian, P. R.; Flory, P. J. Macromolecules 1975, 8, 776.
J. Phys. Chem. B, Vol. 109, No. 39, 2005 18619 (56) Sto¨lken, S.; Ewen, B.; Kobayashi, M.; Nakaoki, T, J. Polym. Sci.: Part B Polym. Phys. 1994, 32, 881. (57) As far as we know, there are no experimental measurements of d ln C∞/d T for syndiotactic polystyrene. Taking the value -0.90 × 10-3 K-1, we used for atactic polystyrene and the value of 10.6 measured at 300 K by Ewen et al. (ref 56), we obtain at 500 K for syndiotactic polystyrene a value of C∞ of 8.8. (58) To have an idea of the simulation costs in terms of CPU time, the simulation of a system of about 600 superatoms (64 chains of nine monomers) for 1 ns using the generalized parallel domain decomposion version of gmq code (ddgmq see refs 60 and 61) on eight processors (SGi 3800, MIPS R14000, 500 MHz) requires 3 min. (59) Brown, D.; Maigret, B. Speedup 1999, 12, 33. (60) Brown, D.; Minoux, H.; Maigret, B. Comput. Phys. Comm. 1997, 103, 170. (61) Auhl, R.; Everaers, R.; Grest, G. S.; Kremer, K.; Plimpton, S. J. J. Chem. Phys. 2003, 119, 12718. (62) Ru¨llman, M., private communication. (63) Green, P. F.; Palmstrøm, C. J.; Mayer J. W.; Kramer, E. J. Macromolecules 1985, 18, 501. (64) Tscho¨p, W.; Kremer, K.; Han, O.; Batoulis, J.; Bu¨rger, T. Acta Polym. 1998, 49, 75.
