Persistence Length, End-to-End Distance, and Structure of Coarse

Subscriber access provided by UniSA Library

Polymers

Persistence Length, End-to-end Distance and Structure of Coarse-Grained Polymers K. Michael Salerno, and Noam Bernstein J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01229 • Publication Date (Web): 12 Mar 2018 Downloaded from http://pubs.acs.org on March 20, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Persistence Length, End-to-end Distance and Structure of Coarse-Grained Polymers K. Michael Salerno⇤,† and Noam Bernstein‡ †NRC Research Associate, Resident at Center for Computational Materials Science, Naval Research Laboratory, Washington, DC 20375 ‡US Naval Research Lab, 4555 Overlook Ave SW, Washington, DC 20375 E-mail: [email protected]

Abstract Coarse-grained (CG) polymer simulations can access longer times and larger lengths than all-atom (AA) molecular dynamics simulations, however not all CG models correctly reproduce polymer properties on all length scales. Here we coarse grain atomistic position data from polyethylene (PE) and polytetrafluoroethylene (PTFE) melt simulations by combining

backbone carbon atoms in a single CG bead. Resulting CG

variables have correlations along the chain backbone that depend on the coarse-graining scale , and are generally not reproduced by independent bond-length, bond-angle and torsion-angle distributions. By constructing distributions of CG variables equivalent to those from simulated CG potentials we are able to evaluate the bond-orientation correlation for di↵erent CG models at reduced computational cost. CG models and potentials that include only nonbonded, bond-length, and bond-angle interactions computed by Boltzmann inversion correctly reproduce the CG variable distributions, but do not necessarily reproduce the chain sti↵ness, overestimating the persistence length Lp and end-to-end distance hR2 i1/2 with increasing . While CG models that include

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

an independent torsion angle match the bond-orientation correlation and hR2 i1/2 better, only approximate models that include correlations between bond and torsion angles match the true bond-orientation correlation.

Introduction statement to reviewer 1 regarding PE distributions Modeling long polymer molecules presents a computational challenge that requires sophisticated modeling techniques. Understanding and modeling polymer length scales is a continuing theoretical and practical challenge. 1–4 In typical applications polymer chains are long and entangled. Though the end-to-end distance and entanglement length govern stress relaxation relevant to processing, many important properties also depend on the persistence length and local bond and angle structure. 5–7 This wide range of length scales and corresponding time scales makes traditional all-atom (AA) or united-atom molecular-dynamics (MD) simulations computationally demanding. Many coarse-grained (CG) polymer models have been developed to decrease the computational cost of simulating real polymers. Since CG interactions derived systematically from AA simulation data generally represent specific polymers better than comparable phenomenological interaction parameters, 8,9 many models of this type have been developed for a range of CG bead sizes, from a small number of repeat units or atoms per CG bead, 5,6,10–14 to much larger CG beads with 16-20 backbone monomers per CG bead. 15–19 Some previous work maps CG bead position to the center-of-mass for a group of atoms, 8,11 while other work chooses specific atom or bond positions, 6,13,20 with many viable choices in principle. 21 Throughout this paper we refer to a particular choice of CG variables for representing a single-chain structure as a CG representation, distinct from a CG model, which we define below. For example, the choice of all bond lengths between pairs of adjacent beads, all bond angles between adjacent CG bead triplets, and all torsion angles between each set of four adjacent CG beads is sufficient to fully determine the relative positions of all the 2


Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CG beads that represent a chain, but do not specify the positions of all the atoms that underlie the CG bead representation. The equilibrium values of the variables chosen for a particular CG representation are in principle described by a joint-probability distribution of all of them. The potential that reproduces this distribution through Monte-Carlo sampling or MD simulation is proportional to the negative of the distribution’s logarithm, and is a function of all the CG variables. Approximating a full joint distribution with a product of lower-dimensional distributions of small subsets of CG variables is equivalent to expressing the potential as a sum of fewbody terms, which is in practice necessary for its computation. Here and throughout this text we use the term model to indicate a set of approximations to the full joint-probability distribution of CG variables for a particular CG representation. In approximating a representation with a model choices are made about which variables and which correlations among them are important. For example, some models of CG polymers define the energy as a sum of bond-length, bond-angle and torsion-angle terms, independent of position along the chain. Such a potential that is a sum of single-variable functions is equivalent to the assumption that the CG variables are uncorrelated, and that their joint-probability distribution can be approximated by a product of single-variable distributions, with no higher order correlations. For any particular choice of approximate CG variable distributions, Boltzmann inversion can be used to calculate a set of interaction potentials that precisely reproduces the chosen probability distributions. 22 For each CG model (i.e. approximate distribution of CG variables) there is an equivalence between evaluating observables using the CG variable probability distributions and the corresponding potential, and each has di↵erent advantages. Using the probability distribution it is simple to calculate the expectation value for any analytical expression of the CG variables. This calculation can be made without developing potentials through Boltzmann inversion. Developing a set of CG model potentials has a cost, and using them to simulate or sample a system has a high computational cost compared to calculating an expectation value directly

3



from the probability distribution. The main benefit of simulating with a CG potential is that nearly any quantity or expectation value can be calculated from these simulations. Previous studies that have included interactions beyond independent bond-length and bond-angle distributions (i.e. additive bond-length and bond-angle interactions) have not conclusively determined when it is important to include torsion angles or consider the correlations among di↵erent CG variables. From the literature, these features appear important on a case-by-case basis. Studies of polyisoprene using iterative Boltzmann inversion (IBI) 13 and polyethylene using conditional reversible work 23 have found torsion-angles significant while IBI studies of poly(vinylalcohol) 5,14 did not. However there have not been studies to determine how the distribution CG torsion-angles depends on the coarse-graining level. Researchers have long noted how some CG degrees of freedom, for example the bond angle and torsion angle in Boltzmann inverted CG models of PE 24 or PS, 25 can be correlated. Further, studies of poly(vinylalcohol) 5 and trans-1,4-bpolybutadiene 6 indicated that correlations may be even more important in the solid state. In this study we develop and compare CG models of PE and PTFE, as examples of structurally simple and similar polymers with significantly di↵erent sti↵ness. Bulkier, sti↵er PTFE, provides a comparison to our previous study of PE, and a check that our current results are valid and independent of the overall sti↵ness of the polymer. We construct CG representations from underlying all-atom (AA) simulations using two di↵erent procedures at three di↵erent coarse-graining scales. Using di↵erent approximate models based on these CG representations we show that changes in the large-scale polymer structure can result from di↵erent approximations to the CG variable probability distributions. We systematically construct CG nonbonded, bond-length, bond-angle, and torsion-angle interactions using IBI to test the uniform (no torsion-angle potential term) and independent (additive torsion-angle potential term) torsion-angle approximations. Simulations using these CG interactions show that approximations that reproduce the CG bond-length and bond-angle distributions but neglect the distribution of torsion-angles do not in general reproduce the persistence length

4


Page 4 of 34

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Lp , end-to-end distance hR2 i1/2 , and entanglement length Ne . Our results indicate that there is a single coarse-graining scale that balances atomistic torsion sti↵ness against the CG bond length and bond-angle sti↵ness to reproduce both local and nonlocal structure without torsion interactions, and that this scale depends on coarse-graining procedure.

Methods Coarse Graining Representations Two procedures were used to generate CG representations from the AA simulations described below. First, an averaging procedure placed a CG bead at rA , the position of the center of P 1 a a mass of every consecutive backbone carbons rA j = i=1 r(j 1) +i , where ri is a backbone

carbon position. Second, a decimation procedure placed a CG bead at the position of the first of every

= ra(j backbone carbons rD j

1) +1 .

Our results refer almost exclusively to

the averaging coarse-graining procedure, with decimation results used only to confirm the generality of certain results. We refer to CG representations generated from underlying atomistic simulations, as A

for averaging or D

for decimation and refer to atomistic

quantities with a superscript “a”. These symbols all represent data generated from all-atom simulations, not data generated from simulations using CG potentials. Throughout this text undecorated symbols have an implied A or a on all quantities, i.e. expressions given in this form are true within any representation. Three di↵erent scales that represent

were chosen with beads

= 2, 3 or 4 CX2 units of mass 28 (100), 42 (150) or 56 (200) g/mol for PE

(PTFE).

All-atom Simulations The A CG representations described above were created based on underlying AA simulation data. AA simulations of PE were reported in previous work. 8,26 For the reference AA data a T = 500 K melt of PE with 345 chains of length N = 96 CH2 groups was simulated using the 5



Page 6 of 34

OPLS-AA potential for over 100 ns, significantly longer than the 10 ns time to di↵use their end-to-end distance. 8 A second AA melt of 216 chains of length N = 480 was simulated for several hundred ns to compare entanglement length and end-to-end distance. Simulations of PTFE used parameters from the OPLS-AA force field. 27 The reference AA simulation contained 2760 chains of N = 96 CF2 repeat units. The temperature was chosen to be 653 K for comparison with previous experimental work. 28,29 The density of the PTFE melt was chosen to be 1.46 g/cm3 , based on the final density from a 1 ns NPT simulation with target pressure 50 atm and barostat time constant 1 ps. For the density chosen, the average pressure at 653 K was P = 13 atm during the production NVT simulation. For both nonbonded Lennard-Jones and electrostatic interactions a 1.2 nm real-space cuto↵ was used. Beyond this cuto↵ electrostatic and attractive Lennard-Jones interactions were computed using a PPPM solver 30 with an accuracy of 4 ⇥ 10

5

relative to the Lennard-

Jones or electrostatic force on two charges separated by 1 ˚ A. For the PTFE simulations the RESPA multi-time step integrator was used to integrate bond, angle, torsion and nonbonded interactions with a 1 fs time step and reciprocal-space interactions with a 4 fs time step. 31 Snapshots of the simulation configurations were saved every 100 ps over the final 80 ns of a 100 ns simulation. Over the initial 20 ns of simulation hR2 i rose by about 10% from its unequilibrated value, while over the subsequent 80 ns quantities measured in this paper as well as hR2 i changed by about 1% with no clear drift. The 100 ps sampling interval was measured to be approximately the decorrelation time of the CG bond-length, bond-angle and torsion-angle variables. The simulation duration was longer than the approximately 55 ns time for the PTFE chains to di↵use their end-to-end distance, indicating adequate sampling of independent chain conformations. Separate N=480 PTFE simulations were simulated for comparing entanglement length. Because it is prohibitive to simulate as long as the Rouse time for this chain length, simulations were run until the chain end-end distance and radius of gyration showed no drift over a 30 ns window. All simulations were performed using LAMMPS. 32

6


Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Coarse Grained Potentials We applied the A coarse-graining procedures to the AA simulation position data, and from each A representation we developed bond-length, bond-angle, torsion-angle and nonbonded interaction potentials using Boltzmann inversion. 33 For the first set of potentials we consider each CG variable to be independent, and the torsion-angle potential to be strictly zero, which is equivalent to a model that assumes a uniform torsion-angle distribution. We refer to models and potentials using this approximation as uniform-torsion models, and label results from these models and simulations using these potentials as BA(A ). Assuming a uniform torsion-angle distribution is a common approximation in many CG polymers, because of a known instability in including an independent torsion-angle interaction when bond-angles have a finite probability to reach 180 . In previous work, 8,26 bond-length and bond-angle interactions for PE were computed through a single Boltzmann inversion step while nonbonded potentials were calculated using an IBI scheme. Here an iterative scheme is also applied to the bond and angle interactions for PTFE, with 2 (3, 4) iterations for bond and 4 (3, 1) iterations for bond-angle interactions for the BA(A2) (BA(A3), BA(A4)) models. The IBI of the PTFE nonbonded interaction required 60 (46, 45) iterations to converge the BA(A2) (BA(A3), BA(A4)) radial distribution function (RDF) to the reference A RDF. For potentials based only on the A representation, converging to the reference A RDF leads to incorrect pressures. We corrected the BA(A2) (BA(A3), BA(A4)) nonbonded potential by applying 4 (5, 6) pressure-correction iterations 34–36 to reduce the average pressure from P = 70 (160, 380) atm to P = 12 (15, 11) atm. These final pressures matched the pressure in a long-time AA NVT simulation. We also develop models and potentials that assume that the torsion-angle is independent from other CG variables, but not uniformly distributed. We refer to models and potentials using this approximation as independent-torsion models, and label results from these models and simulations using these potentials, which include additive bond-length, bond-angle, and torsion-angle terms, as BAT (A ). The nonbonded and bond-length interactions are identical 7



to those in the BA(A ) models above. Due to the numerical instability described above, the form used for the torsion-angle potential is not strictly independent, but is actually a function of the adjacent bond-angles ✓ (as described in more detail below), thus an iterative scheme that alternates the bond-angle and torsion-angle updating was employed. Optimization of the potentials required 4 (3, 1) bond-angle iterations and 0 (1, 0) torsion-angle iterations for the BAT (A2) (BAT (A3), BAT (A4)) models. Using the BA and BAT (A ) potentials, polymer chains with lengths 24-480 beads were simulated. These chain lengths allowed us to compute properties like hR2 i and the primitive path length hL2pp i1/2 for longer chains. 37,38 Below we also include results for CG models labeled BA0 (A4) with angle interactions phenomenologically weakened to match the AA Lp . For PE we also include a model BA00 (A4) with angle interactions weakened to match hR2 i.

Cuto↵ Torsion Potential Torsion angle interactions are typical of all-atom force fields that are used to model extended molecules. Unfortunately, in many coarse-grained models the bond-angle interaction has a minimum near the ✓ = 180 , configuration. This feature leads to a finite probability for the bond-angle to be arbitrarily close to 180 and for the resulting torsion-angle interaction to produce rapidly varying forces over a very short distance. This near-discontinuity has been a problem for researchers working in CG modeling, 39 and several solutions have been developed to regularize torsion-angle interactions so that they may be used in CG models. The solution of coupled bond-angle/torsion-angle interactions have long been used in all-atom models to better reproduce certain qualitative features. 40 Our goal is to create CG models with independent potentials of the bond-angle and torsion-angle that are well-behaved for bond-angles near 180 . Essentially, this requires cutting o↵ the torsion-angle interaction for a range of bond-angles by multiplying the potential Ud (') by a cuto↵ function. We have tested two di↵erent cuto↵ forms, one that is constant

8


Page 8 of 34

Page 9 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in the bond-angle ✓ except for a small range of ✓ near 180 where it decreases quadratically:

f (✓) =1 ✓ f (✓) = 1

(✓

⇡ + ✓) ( ✓)2

✓

Persistence Length, End-to-End Distance, and Structure of Coarse

Recommend Documents