An analytical coarse-graining method which preserves the free energy, structural correlations, and thermodynamic state of polymer melts from the atomistic to the mesoscale

An analytical coarse-graining method which preserves the free energy, structural correlations, and thermodynamic state of polymer melts from the atomistic to the mesoscale

J. McCarty    A. J. Clark    J. Copperman    M.G. Guenza111Author to whom correspondence should be addressed. Electronic mail: mguenza@uoregon.edu Department of Chemistry and Institute of Theoretical Science, University of Oregon, Eugene, Oregon 97403
July 3, 2019
Abstract

Structural and thermodynamic consistency of coarse-graining models across multiple length scales is essential for the predictive role of multi-scale modeling and molecular dynamic simulations that use mesoscale descriptions. Our approach is a coarse-grained model based on integral equation theory, which can represent polymer chains at variable levels of chemical details. The model is analytical and depends on molecular and thermodynamic parameters of the system under study, as well as on the direct correlation function in the limit, . A numerical solution to the PRISM integral equations is used to determine , by adjusting the value of the effective hard sphere diameter, , to agree with the predicted equation of state. This single quantity parameterizes the coarse-grained potential, which is used to perform mesoscale simulations that are directly compared with atomistic-level simulations of the same system. We test our coarse-graining formalism by comparing structural correlations, isothermal compressibility, equation of state, Helmholtz and Gibbs free energies, and potential energy and entropy using both united atom and coarse-grained descriptions. We find quantitative agreement between the analytical formalism for the thermodynamic properties, and the results of Molecular Dynamics simulations, independent of the chosen level of representation. In the mesoscale description, the potential energy of the soft-particle interaction becomes a free energy in the coarse-grained coordinates which preserves the excess free energy from an ideal gas across all levels of description. The structural consistency between the united-atom and mesoscale descriptions means the relative entropy between descriptions has been minimized without any variational optimization parameters. The approach is general and applicable to any polymeric system in different thermodynamic conditions.

I Introduction

Molecular dynamics (MD) simulations have become instrumental in developing and testing theories in polymer physics since they can provide direct access and physical insight into the time evolution of complex fluids. However, explicit atom MD simulations are computationally costly which limits their range of applicability to small length and time-scales. This limitation has stimulated the development of numerous coarse-graining methods, which are highly efficient because they represent the system at a lower resolution, and thus greatly reduce the number of degrees of freedom.Hall (); Dama (); Martini (); Arun (); Kremer (); depablo () Nonetheless, a major limitation in coarse-graining, which has delayed the widespread use of such models in engineering and material science is the lack of thermodynamic consistency between various representations. For many numerical coarse-graining schemes, there is no guarantee that the resulting behavior of the coarse system will be consistent with what would have been observed by using a more detailed, more-expensive, model.JOHNSON (); LOUIS (); FALLER2008 () Much of the difficulty in developing consistent coarse-graining approaches stems from the fact that it is not always clear how various many-body effects are incorporated into simple two-body interactions, and how errors in the numerical optimization of a pair potential are propagated to thermodynamic properties.Anthony2 ()

In a multiscale procedure, one seeks to link simulations at different levels of description into an “unified” description of the liquid at all length scales. This requires understanding, in a well-defined and clear way, how to represent a molecule into coarser and coarser units successively, and then to reintroduce the finer structure of the liquid “a posteriori” in a reverse mapping procedure.Kremer () This procedure needs to be done with coarse-graining models that ensure the consistency of the structure and the thermodynamics, independent of the degree of coarseness of the molecular model adopted. Conserving both the thermodynamic and dynamical properties while representing the same system with two different degrees of structural refinement, i.e. atomistic and effective unit models, is not trivially done, with the consequence that even the most popular coarse-graining techniques generally produce systems that have the correct structure but too high compressibility and in most cases incorrect free energy. This deficiency can have important consequences when coarse-grained models are used in mesoscale molecular dynamics simulations to study large systems approaching their phase transitions, or the thermodynamic conditions that lead to the possible emergence of kinetically trapped metastable phases, that have geometries desirable for their engineering applications.

In a series of papersYAPRL (); YA2005 (); SAM2006 (); SAMB (); Anthony (); Anthony2 () we have developed a coarse-graining method which is based on statistical mechanical principles, i.e. the Ornstein-Zernike equation and liquid state theory. The structure of the integral equation theory warrants the proper calculation of the effective pair interactions between coarse-grained units, which emerge from the propagation through the liquid of the many-body atomistic interactions. Recently, we demonstrated in detail how our coarse-graining procedure is thermodynamically consistent with the atomistic description, within polymer integral equation theory, for a soft-sphere coarse-grained representation where we used the PRISM thread model to parameterize the atomistic model.McCarty1 (); PRISM (). We also presented a detailed, formal study of the effective potential, when a polymer is represented as a chain of soft blobs with variable size. Using the analytical solution for the effective pair potential acting between coarse-grained units, we have demonstrated that the thermodynamic inconsistency frequently observed in numerical coarse-graining schemes arises from a failure to properly consider the long range repulsive tail and related attractive well of the effective potentialAnthonyPRL (). It is our contention that with a reasonable molecular model along with the correct parameters, both structural and thermodynamic properties can be simultaneously preserved in coarse-graining, as it is in our approach, and that this can be achieved without recourse to any numerical re-optimization scheme.

The purpose of this paper is to implement our variable-level coarse-graining model to derive thermodynamic and structural properties for the variable levels of our coarse-grained description, to investigate consistency and inconsistency of the thermodynamic properties, and to demonstrate the versatility of the approach to model realistic polymer chains. The approach taken here is slightly different from our previous work. In our previous work, we obtained an analytical expression of the effective coarse-grained potential, that depends on the monomer parameter, . Essentially, the parameterization of the potential reduces to finding the single monomer quantity, . In this work, we propose to obtain this quantity from the solution to the numerical PRISM equations for a homopolymer semi-flexible chain. This model requires an effective hard sphere diameter, which is chosen such that the pressure calculated from the predicted equation of state agrees with the observed pressures in united atom simulations. Once this parameter is known, the potential is uniquely specified.

Using the calculated potential we perform simulations of the coarse-grained systems. We then compare thermodynamic quantities and structural quantities of interest from these coarse-grained simulations with united atom simulations. Figure 1 shows a schematic of the method. There are two sets of simulations: coarse-grained and united atom. The united atom simulations are used only as a test of the method, not too provide information to the coarse-graining approach but to compare thermodynamic properties with the coarse-grained simulations of the same system. Tests are done for a variety of coarse-graining levels, and for systems in different thermodynamic conditions and variable chain length.

Figure 1: Schematic outline of the approach used in this paper. There are two different types of simulations: united atom and coarse-grained. The coarse-grained potentials are obtained by solving the numerical PRISM equation for the parameter , which serves as an input into the analytical theory. The potential is then used in coarse-grained simulations with variable numbers of effective sites. Structural and thermodynamic quantities are then compared directly to united atom simulations.

Here we consider linear chains of polyethylene melts at various chain lengths and various densities where we represent the chains at different level of chemical detail. We show that our coarse-graining procedure is general and applicable to a wide number of interesting systems. Furthermore, the study presented here shows that our approach is useful in multi-scale simulations of dense polymer systems with specific chemical structure because it is system specific and reproduces the correct equation of state across various levels of coarse-graining.

The paper is divided as follows. We begin Section II with a brief discussion of the background theory that summarizes the procedure to compute the effective potential, in terms of the parameter, . In Section III, we introduce the numerical solution of the PRISM approach to calculate the direct correlation function parameter, , which is an input parameter to the coarse-grained representation. Then, we summarize the procedure that we follow to perform molecular dynamics simulations of the coarse-grained liquid, and the united atom molecular dynamic simulations that are used to test the theory. In the following sections, the paper presents calculations of thermodynamic properties, specifically structural correlation functions, isothermal compressibility (Section IV), equation of state, Helmholtz and Gibbs free energies, potential energy and entropy (Section V). For each level of coarse-graining both the analytical formalism for the thermodynamic properties, and the numerical solution of the thermodynamics using mesoscale simulations of the coarse-grained system, as well as the calculation of the thermodynamic properties by united atom molecular dynamic simulations are presented. The purpose of this cross testing using multiple approaches is to demonstrate, without any ambiguity, the consistency (or inconsistency) of the thermodynamic properties for the variable-level coarse-graining model, independent of the chosen level of coarse-graining. The study confirms the agreement of the structure and thermodynamics of the coarse-grained representation, with the united atom description of the same liquid through traditional MD simulations, and validates the analytical description of the thermodynamic properties presented in this paper.

Figure 2: Coarse-graining at multiple block-levels. The far right is a snapshot of a typical configuration from an United Atom Molecular Dynamics simulation of polyethylene. Moving from right to left the same chain is represented by 5 blocks, 3 blocks, to a single soft sphere. Below is a snapshot of a typical mesoscale simulations, where each sphere represents a single polyethylene chain.

Ii Background Information: Variable-level Coarse-Grained Description of Polymer Melts

Our course grained description represents each polymer chain in a melt as a single soft sphereYAPRL (); YA2005 (); SAM2006 () or as a collection of soft connected blobsSAMB (); Anthony (); Anthony2 (). While our previous papers provided a theoretical derivation of the variable-level coarse-grained representation, here the numerical procedure is only summarized for completeness. In a series of snapshots Figure 2 shows how a single polymer chain is represented first by several spheres, and then by a single soft sphere, allowing for the simulation of many more polymers than would otherwise be possible to simulate.

We take the polymer to be comprised of equivalent sites of monomers with a chain density , and an effective segment length with being the polymer radius of gyration and the blob radius of gyration. The number of underlying monomers per blob is given as . Coarse-grained or fictitious interacting sites are taken to be located at the center of mass of the polymer chain for the soft sphere model or at the centers of mass of several monomers along the same chain for the connected blob model. The relation between center of mass fictitious sites and real monomer sites, originally proposed by Krakoviack et alKrakoviack () is derived by solving a generalized matrix Ornstein Zernike equation, and was extended by Clark and Guenza Anthony () for the multi-blob model. We now summarize the model for the case of one, three, and multiple effective sites (see Figure 2).

ii.1 Soft Sphere Potential

For the case where the number of blobs is one, meaning each chain is represented as a single interaction site, one recovers the relation,

(1)

where is the total correlation function between center of mass, fictitious sites. is the total correlation function between monomers, and and are the intramolecular distributions between monomers or between monomers about the center of mass respectively. The monomer total correlation function, , is given in Fourier space by the PRISM integral equation for a homopolymer fluidPRISM1 (); PRISM2 (),

(2)

where is the number density of monomer sites, and is the monomer direct correlation function. As an approximation, valid at large r, we take the direct correlation function to be a constant in Fourier space, . Substitution of Equation 2 into Equation 1 gives,

(3)

For the soft sphere model, the direct correlation function between coarse-grained sites, , is identical to that for a simple fluid,

(4)

The effective soft sphere potential is given by the hypernetted chain (HNC) closureHansenMcDonald (),

(5)

The monomer intramolecular distribution is given by the Debye formula,

(6)

whereas the monomer-cm distribution is given by

(7)

The analytical solution of for polymer melts, i.e. high density mean-field solution, has been presented in a recent publication, together with the study of the scaling behavior.Anthony2 () In this paper we calculate the potential numerically. The procedure formulated in this manner is easy to implement. Using the definitions of the intramolecular and intermolecular distribution functions, is simply calculated and its subsequent numerical Fourier gives the effective potential from Equation 5. This potential is used as an input to Molecular Dynamics simulations of soft spheres. Note that the formalism only depends on molecular ( and ) and thermodynamic parameters (, , and ), as well as on the system compressibility through the direct correlation function at , .

ii.2 Tri-Block Potential

For chains represented by three coarse grained sites per chain (tri-block model) the two terminal coarse grained sites are distinct from the central coarse grained site giving rise to three distinct block combinations of block-block total correlation functions. Using the notation from our previous workAnthony () the intermolecular total correlation functions between block centers, with block indices and , are given as a function of the new omega distributions for the tri-block structure as

(8)
(9)

and

(10)

The direct correlation functions are then expressed in term of the intramolecular distributions by inserting Equations 8-10 into the matrix Ornstein-Zernike equation of the coarse-grained fluid. The resulting direct correlation functions are lengthy but presented here in their full form,

(11)
(12)
(13)

with

(14)

In order to to obtain the tri-block potential, we use the results for the tri-block intramolecular distributions calculated by Clark and GuenzaAnthony2 (). The monomer-monomer distribution is given by the generalized Debye formula,

(15)

For the block-monomer and block-block distributions, distinct block sites are represented by Greek indices and and the block separation is given by . For the tri-block model and the distributions of monomers about blocks on the same site, , is given as

(16)

whereas the block-monomer distribution of monomers on block, , with block is

(17)

Lastly, the block-block distribution is given as

(18)

Once again Equations 15-18 are used to derive expressions for and , using Equations 8-10 and Equations 11-13. The potential is then calculated numerically after Fourier transform using the HNC closure relation, yielding three potentials between block types, , , and , depending if the block is inside the polymer or at its ends.Anthony2 ()

The formalism just presented to describe position-specific correlation functions for tri-block renormalization can be extended to treat five-block renormalization.Anthony2 () The distribution functions are formally complex, but the procedure of deriving them is straightforward.

ii.3 Block Averaged Potential

Although the expressions for site specific block properties are complicated, when the number of blocks is large, the formalism can be drastically simplified. Specifically we have observed that when there are more than five coarse grained sites per chain, end effects become negligible and it is possible to use a block-averaged description which simplifies the expressions. A detailed analysis of the properties of the block-averaged potential is presented in a recent paper together with the analytical solution of the expression for the potential and the analysis of the scaling exponents.Anthony2 ().

We use here a numerical procedure to calculate the characteristic distribution functions and the potentials, which enter the molecular dynamic simulations. The numerical procedure begins the same way by calculating the intramolecular distributions in the block averaged limit. To maintain consistency with the expressions in our previous workAnthony2 (), we introduce the normalized intramolecular distribution , for the block-block (bb), block-monomer (bm), and monomer-monomer (mm) distributions. The normalized block-monomer and the block-block distributions in the chain averaged limit are

(19)

and

(20)

The monomer distribution is again given by equation 15 normalized by .

The next step is to use the result from solving the matrix Ornstein-Zernike equation for ,

(21)

and, from the Ornstein-Zernike relation

(22)

in the direct correlation function, calculated as

(23)

While the analytical solution of the potential has been presented earlier,Anthony2 () in this paper we solve the potential numerically. First, the direct and total distribution functions in real space are calculated from numerical Fourier transform of Equation 22 and Equation 23. Then, the interaction potential is calculated numerically evaluating the HNC potential, as above,

(24)

Iii Methods

iii.1 Numerical PRISM Solution for the Calculation of the Direct Correlation Function

The only non-trivial parameter in our formalism is the direct correlation function in the large scale limit.McCarty1 () For a realistic representation of the polymer chain it is useful to work with the numerical solution of the PRISM equations, which are solved to obtain the monomer direct correlation parameter, , which is then used as an input in the expressions of the coarse-grained potentials. This is the only non-trivial parameter that enters our approach, i.e. the only parameter that does not, trivially, describe physical or molecular quantities. Other parameters, besides the direct correlation , are the thermodynamic properties of temperature, , and density, , as well as the structural properties of , the number of monomers, and the effective segment size, . These parameters allow the method to be readily applied to a variety of polymers in variable experimental conditions.

An analytical definition of the parameter was obtained using the Gaussian thread model, which relies on the description of the polymer chain as infinitely long and infinitely thin, and is the model used to represent polymers in field theory.McCarty1 (). While the PRISM thread model represents an idealized limiting case, it is not expected to give quantitative predictions for real chains of finite length and thickness. In this paper we use a self-consistent numerical solution of the PRISM equation.

The monomer total correlation function, , is given in Fourier space by the PRISM integral equation for a homopolymer fluidPRISM1 (); PRISM2 (),

(25)

where is the number density of monomer sites, and is the monomer direct correlation function. We adopt the semiflexible chain modelHonnellCurro () as depicted in Figure 3. In its simplest form, the model requires three parameters, the monomer hard sphere diameter, , the bond length, , and the bond angle, . For polyethylene chains we use the values of and to agree with the stiffness , which has been shown to be reasonable for linear polyethylene models. To estimate the intramolecular distribution, we then use the Koyama distribution which can be calculated using the method from reference HonnellCurro ().

Figure 3: Semiflexible chain model used as the molecular model for numerical PRISM equations. This model allows for a realistic estimate of , which is input into the coarse-grained model.

Having the intramolecular distribution, we solve the PRISM equation, Equation 2, with the reference molecular mean spherical approximation (R-MMSA) closure given byYethiraj1 (); Yethiraj2 ()

where the asterisks denote convolution integrals, is the direct correlation function of a reference purely hard sphere chain, and is the attractive part of the potential given as a Lennard Jones tail,

Inside the hard core, impenetrability is ensured by enforcing the condition

(26)

Once the closure is specified by Equations III.1 - 26, the PRISM equation is solved numerically using the standard Picard iteration with fast Fourier transform or with the KINSOL nonlinear algebraic solver available through the SUNDIALS suite of programs.sundials ()

Once a solution of is achieved, the value of is determined as

(27)

In performing numerical PRISM calculations we adjust the hard-sphere diameter, , to agree with the predicted equation of state given by Equation 40, so that the pressure agrees with United Atom simulation data and determined from Equation 27. We also check that this value of yields a good representation of the monomer structure as given by the radial distribution function, . Figure 4 shows calculated from PRISM theory and from United Atom simulations for and a monomer density if sites/.

Figure 4: The monomer radial distribution function for PE with at a monomer density sites/ calculated from United Atom simulations (circles) and from PRISM calculation (solid line) using the SFC model with which gives the best agreement between the pressures

.

Figure 5 shows a plot of for different densities for a chain of , , and . Most noticeably, the dependence of on density is stronger than the linear dependence predicted by the idealized thread model. The right side of Figure 5 shows the as a function of chain length for sites/. The value of approaches an asymptotic value for large which corresponds to a leveling off of the chain-length dependence of the pressure. This limiting value of can be used to perform coarse-grained simulations of large chains where all-atom simulations are prohibitive.

It is interesting to compare the effective pair potential obtained from the PRISM calculation using the semi flexible chain model (SFC) with that obtained from the thread model. Figure 6 shows the effective pair potential for the intermediate density, sites/, where the molecular model input into the theory is calculated from PRISM theory for both the semi flexible chain model and the thread model. The SFC model results in a stronger and longer-ranged repulsive core; however the attractive well is small and of comparable magnitude. The inset shows the virial force which is used to calculate the pressure equation. The repulsive contribution is much larger for the case of the semi flexible chain model. This is due to the finite-sized hard sphere core in the semi flexible chain model, which results in a stronger repulsive force. Because the radial distribution function is fairly insensitive to the shape of the potential, the structure is nearly identical between the thread model and the semi flexible chain model; however, because the pressure depends on the virial force, and it is extremely sensitive to the potential, even small deviation in the shape of the force can lead to different compressibility, and thus the semi flexible chain model should be used since the value of chosen is such as to match the pressure from the united-atom simulations.

Figure 5: Left: Calculated value of from numerical PRISM using a simple semi flexible chain chain model as a function of the monomer site density. Three different chain lengths are shown; (circles), (squares) and (triangles). The line is a fit to a quadratic polynomial to serve as a guide to the eye. Right: Calculated values of from numerical PRISM as a function of polymer chain length at fixed density, sites/ (circles). The solid line is a fit to the form .
Figure 6: The effective pair potential for PE100 at an intermediate density of sites/ The solid line is the potential calculated using numerical PRISM with the SFC chain model. The dashed line is the same potential using the thread model value for for comparison. The inset shows the virial force.

iii.2 Mesoscale Molecular Dynamics Simulations

In our coarse-grained models each polymer chain is represented as a collection of soft spheres, where each sphere includes a large section of the polymer chain, which is larger than the polymer persistence length. As a limiting case, we represent the whole chain with a single, soft interaction site, i.e. the soft-sphere model. For this level of coarse-graining, the relevant information is on length-scales larger than a single polymer chain, and the model affords the largest computational gain.

In general, we call simulations performed at these various multi-block levels of description Mesoscale (MS) simulations, since they represent the system at a level of description in between an all-atom and a continuum level. Having a description that systematically represents the chain at variable levels of detail, enables us to connect information in a multi-scale procedure.Jaymultiscale () Mesoscale simulations of point particles allow for the simulation of many more particles than the united-atom simulations. The increased particle number makes these simulations suitable to study bulk, large-scale properties such as shearing or phase transitions.

The soft particles interact via the coarse-grained potentials described in the previous section. For the multi-block model, bonded sites are given a bond potential derived from the direct Boltzmann inversion of the probability distribution of the effective bond length,

(28)

where

(29)

The angle potential is similarly computed

(30)

where the angular probability distribution for a random walk chain is given byLasos ()

(31)

with for long chains. All intramolecular sites separated more than two atoms apart interact via the pair potential of Equation 24.

Initial configurations of particles were generated on a cubic lattice, and classical molecular dynamics simulations were performed using the LAMMPS MD codeLAMMPS (). Due to the soft nature of the pair potential, mesoscale simulations rapidly equilibrate which is a further advantage to UA simulations of the same system. Simulations were run in the NVT ensemble to reproduce the conditions of the United Atom simulations.

iii.3 United Atoms Molecular Dynamics Simulations

The above sections completely describe the coarse-grained model. It is important to note that there is no need to parameterize the model against structural correlations from any more detailed simulation (united atom or atomistic). However, we want to demonstrate that the model can provide quantitative agreement with more detailed simulations. Therefore, in this paper we make direct comparison to united atom simulations, since the parameters of this model for polyethylene are well known. In the united representation, each group is represented as a single interaction site.UAmodel () This model has been used extensively to represent hydrocarbons of different chemical specificity.TRAPPE (); Jaramillo (); Jaramillo1 () We ran two sets of simulations: the first set at with a variety of densities and chain lengths, specifically , , , , , each with the set of densities , , , , , and sites/. The second group of simulations was at at a fixed density of sites/ and a variety of chain lengths, specifically , , , , , , , , and .

The UA simulations use the TRAPPE forcefield parameters,TRAPPE () and are performed using the LAMMPS molecular dynamics codeLAMMPS (). Chains were randomly generated, and overlapping chains in the initial configuration were slowly pushed apart by a soft repulsive potential for 1 ns of simulation time. Then, the full non-bonded potential was switched on with a small time step, and the system was run for an additional 1ns while ramping up the timestep to 1.25 fs. Subsequently, chains were allowed to equilibrate for 10 ns before a final production run of 11 ns was used to collect trajectories and calculate static properties and thermodynamic averages. Although these trajectories are shorter than the longest relaxation times of long chains, here we are not concerned with dynamics; therefore these trajectories provide ample time for static properties to be collected. Simulations were run in parallel using the SDSC Trestles clusters available through the XSEDE program. All simulations were run in the NVT ensemble using Nose Hoover thermostat.

Lennard-Jones (kcal/mol) ()
0.0912 3.95
Bond parameters () (kcal/mol )
1.54 900
Angle parameters (degrees) (kcal/mol rad)
114.0 123.75
Dihedral parameters (kcal/mol) (kcal/mol) (kcal/mol)
1.4110 -0.27084 3.143
Table 1: Polyolefin Melt Simulation Parameters

Non-bonded interactions between sites separated by more than 4 sites away are governed by the Lennard-Jones potential,

(32)

with the parameters specified in Table 1. Adjacent intramolecular sites interact through a harmonic potential

(33)

where is the bond length. Angle interactions between triplets of intramolecular sites are governed by a harmonic potential of the form

(34)

Finally, dihedral interactions between quadruplets of atoms is governed by the OPLS dihedral potentialOPLS ()

(35)

Iv Structural Correlation Functions and Isothermal Compressibility

The particles in the coarse grained model correspond to fictitious center of mass sites of the underlying real chain, either at the single soft sphere level, or at the mutiblock level. In other words, many atoms in the underlying detailed description are “mapped” onto a single site in the coarse grained representation. Thus, the structure of the coarse grained model reproduces the structure of the detailed model only on the level of the mapping, and the structures are equivalent if the distribution functions are equivalent at that level of description and larger. Many atomistic configurations could have the same center of mass configuration, and a coarse grained procedure is unable to distinguish between these “finer” levels of detail. This leads to a reduction in the dimensionality of the configuration space which is also associated with a smoothing of the probability distributions which have a consequence for the entropy of the system, which will be discussed below.

Here we are interested in comparing the structural distributions between coarse-grained and united atom simulations in the configuration space of the coarse-grained model. Figure 7 shows a depiction of a random-walk chain mapped onto a tri-block chain with three interaction sites located at the center of mass of a group of monomers. The bond vector linking coarse-grained sites is given as R and the angle, , defines the angle between two bond vectors, and . As a representative case, we map a polyethylene chain with monomers onto a 3-block model, each block with 75 underlying monomers. The top panel of Figure 8 displays the intramolecular block-block correlation function in reciprocal space, for the central block in a chain represented as a tri-block chain. The distribution is calculated with the analytical expression of the tri-block model, and directly compared to the united atom simulation, and the mesoscale simulation of the coarse-grained, tri-block representation.

The bottom panel of Figure 8 shows the intermolecular correlation function, , between block centers from the analytical expression, and both the coarse-grained and united atom simulations. The agreement is quantitative, with the only parameter that enters the coarse-grained potential, , calculated from numerical PRISM. Typically the ratio between this distribution computed in both the coarse-grained and united atom simulation would be used in an iterative numerical scheme to optimize the pair potential or the related force. Our approach does not require any optimization since the two distributions are equivalent. This is because the coarse-grained potential reproduces the underlying free energy surface along the coarse-grained degrees of freedom and thus reproduces the correct probability distributions of coarse-grained coordinates. All the structural properties display consistency independent of the level of coarse-graining.

Figure 7: Depiction of the three block model mapped onto a random walk. The coarse-grained coordinates are and , which are the bond vectors joining effective coarse grained sites, and the angle, , between them. Note that many underlying atomistic configurations are possible for each configuration of the coarse-grained coordinates.
Figure 8: Comparison of the intramolecular distribution, (top) and intermolecular correlation function, , (bottom) between block centers calculated from both 3-block and united atom simulations of PE225. The solid line depicts theoretical predictions from the coarse-grained model.

It is straightforward to show that the compressibility is also preserved in coarse-graining. The isothermal compressibility, , is calculated from liquid state theory as

(36)

with the limit of the static structure factor . Taking this definition in the two levels of description, for the monomer (m) and for the blob (b) respectively, we have

hence

(38)

for any level of block description. We have used here the property that , as can be seen from Equation 21 given that , as required by the normalization. As it is expected, the compressibility of the liquid is independent of the level of detail that is used in the representation of the polymer chain. In this way, the described coarse-grained representations can be conveniently adopted in molecular dynamic simulations of polymer liquids, together with the related potentials described in the previous section: this coarse grained description allows one to select the most convenient length scale for the molecular representation, which depends on the properties that one wants to study, while being confident that both structure and compressibility are correctly represented. More thermodynamic properties are studied in the following section.

V Equation of State, Free Energy, Internal Energy, and Entropy

This section shows that not only the structure and the isothermal compressibility, but also the thermodynamic quantities of pressure, and excess free energy are conserved quantities across the variable levels of representation. In other words, the multi-scale scheme is self-consistent across multiple block levels of description and once the parameters are correctly chosen reproduces the correct equation of state as those predicted from UA simulations. We also show that other thermodynamics quantities, such as the internal energy and the entropy depend on the level of coarse-graining, as it was expected.

v.1 Equation of State

As a first test of the capability of the coarse-grained potential to provide thermodynamic consistent simulations across multiple levels of representation, we performed a series of simulations of both soft-blobs with variable levels of coarse-graining and UA models under the same set of conditions. The pressure is calculated from each simulation using the standard viral expression, with the total number of monomers

(39)
Figure 9: Comparison of the pressure from united atom simulations of PE (black circles), PE (squares) and PE ( black triangle) with coarse-grained simulations of the soft sphere model (stars) and tri-block model (green triangles) and 5-block model (red squares). Lines are a guide to the eye.

Figures 9 and 10 show the pressure from polyethylene simulations performed at for a variety of chain lengths and densities. In all cases the coarse-grained simulations reproduce quantitatively the trend of pressure as a function of density of the equation of state without any post-optimization scheme or fitting procedure, and across variable levels of coarse-graining. We also consider a set of simulations of increasing chain length performed at a constant monomer density of and temperature, . As united atom simulations of long chains are difficult to equilibrate, a reasonable starting configuration for these simulations were provided to us by Mavrantzas and coworkers.Mavrantzas (); Mav2 () Comparison of theses simulations with various levels of coarse-graining are presented in Figure 11 and shows consistency of the pressure for all the samples.

Figure 10: Comparison of the pressure from united atom simulations of PE (black circles), PE (black squares) with coarse-grained simulations at the soft sphere level (stars) and 3-block level for PE100 (green triangles). Lines are a guide for the eye.

We now present a formal analysis of the equation of state in the variable-level coarse-grained representation. For our model the effective potential between coarse-grained units can be accurately approximated, in the mean-field limit of a polymer liquid, by an analytical expression, which leads for the equation of state toMcCarty1 (); AnthonyPRL ()

(40)

where is the site number density, is the number of monomers per chain, is the density of chains, and is the direct correlation function at zero wave vector, , defined from liquid state theory and discussed earlier on in this paper. This equation of state holds for any level of blob descriptionAnthonyPRL (). All of the non-ideal contributions to the pressure that arise from highly system-specific interactions are contained in this single parameter, . As a result, the general expression of Equation 40 is independent of the model used, but evaluation of other thermodynamic properties from the pressure requires knowledge of the state dependence of for the given system of interest. Previously, we used the PRISM thread model as an informative limiting case, which predicts a linear dependence of on the monomer density, giving the equation of state as,McCarty1 ()

(41)

However, from numerical calculations presented above, the density dependence of is not linear for higher densities, and the analytical result grossly underestimates the pressure for high density systems where excluded volume effects, not included in the thread model, dominate the pressure. We now consider a convenient expression for that allows for the calculation of the pressure and other thermodynamic properties, which can then be tested by comparison to simulation results. When introduced into Equation 40, this provides an expression for the equation of state for the polymer chain.

In simple liquids, the compressibility factor () is often plotted against the packing fraction, , where is the hard sphere diameter. Although real polyethylene chains are not simple hard sphere fluids, a similar analysis can be made by defining an effective packing fraction where is an “effective” hard sphere diameter and is the monomer density. The effective hard sphere diameter will be smaller than the real hard sphere diameter because of the overlaps between intramolecular sites, which decreases the amount of space occupied by the molecules. We assume that for polymer chains, the direct correlation function can be written as a power series in the effective packing fraction,

(42)

Substitution into Equation 40 gives

(43)

Following the intuition that for a simple hard sphere fluid, the series can be written as a linear combination of only the first and second derivatives of a geometric seres, we factor out a polynomial, , such that the remaining terms are given by the binomial coefficients of the form,

(44)

which can be summed to give,

(45)

Next, we assume to be a quadratic polynomial in terms of , such that [], giving the final result for the equation of state for the polymer,

(46)

and the direct correlation function,

(47)

Equation 46 has the form of a Carnahan Starling-type expression, but with numerical pre-factors in front which reflect the chain connectivity and the fact that the real potential is not a hard sphere potential. For large N, the equation of state can be approximated as,

(48)
Figure 11: Pressure vs. degree of polymerization calculated for a hierarchy of soft-block simulations as compared to united atom simulations. All simulations were carried out at constant temperature, T=509K, and density, g/mL. United atom model is depicted with black circles, soft sphere model with blue asterisk, tri-block with green triangles, 5-block with red squares, 10-block with maroon diamonds, and 20-block with orange left-oriented triangle. The line is the theoretical prediction of Eq. 40 with dfined as in the right panel of Figure 5.

Assuming that and are independent of and temperature, a plot of the pressure against the effective packing fraction should yield a universal pressure curve for all polyethylene chains. This plot is shown in Figure 12 where the pressure calculated from united atom simulations is presented as a function of the effective packing fraction. The points nearly fall onto a universal curve of the form of Equation 48. A fit of all the data yields the empirical values of and , which are independent of both temperature and degree of polymerization. The plot suggests that by using Equation 47 we can estimate for any chain length at any temperature for polyethylene.

Figure 12: Pressure data from United Atom simulations plotted as function of the effective packing fraction. The orange stars are the simulations at each with different chain length, and all the other points are simulations at with different chain lengths and different densities. The dot-dashed line is Equation 48, which is independent of the degree of polymerization N.

v.2 Free Energy Estimation

Having a form of the equation of state allows for an expression of the Helmholtz free energy change associated with a change in packing fraction from , to be obtained by integration

(49)

where is the number of polymer chains. Substitution of Equation 46 into Equation 49 gives an expression for the free energy change per monomer

(50)
Figure 13: Estimated free energy change obtained via thermodynamic integration of the pressure, as a function of packing fraction compared to a reference packing fraction of , which is the lowest packing fraction that was simulated. All systems collapse to a universal curve within numerical precision, independent of degree of polymerization. The solid and dashed lines represent equation 50 evaluated at N=44 and 200 respectively.

Taking the limit as gives a closed expression for the excess Helmholtz free energy per monomer, associated with the liquid as distinct from an ideal gas

(51)

In the high density limit, in which we are interested, every particle interacts with a large number of surrounding molecules and the mean-field approximation applies. Since the coarse-grained potential is long ranged and bounded without short range excluded volume interactions, the excess Gibbs free energy in the canonical ensemble is calculated from the mean field approximation,Likos (); Evans ()

(52)

with the blob density and the number of monomers per blob,

(53)

which is independent of the level of coarse-graining as , the number of monomers per chain. The excess Gibbs free energy per monomer is given by . Equation 53 and Equation 50 are different because of the density dependence of the potential and . Equation 50 is derived by integration over the density, thus the density dependence of had to be accounted for. Equation 53 was calculated in the canonical ensemble where the density is constant, and thus the potential is not changing. The main point here is that the excess free energy, in both ensembles, does not depend on the level of coarse-graining, and is a constant, which is consistent with Figure 13 and with the fact that the computed pressure is independent of the level of coarse graining.

As an example, we performed both MS and UA simulations of short chains at constant volume and temperature ( and monomer density sites/). Using the identity,

(54)

the change in free energy was calculated numerically from the simulations, by integration of the pressure as a function of the number of monomers per chain and is shown in Figure 14. As a test for consistency, the numerical values are compared to the predicted expression from using Equation 40,

(55)

for which quantitative agreement between simulation data and analytical theory is observed.

Figure 14: The change in free energy per monomer associated with a change in the number of monomers per chain from a reference chain of 36 monomers as calculated from United Atom simulations (black circles) and Mesoscale simulations (open green circles) by integration of the pressure at constant volume and temperature. The solid line is the predicted free energy change from equation 55

v.3 Potential Energy in the Variable-Level Coarse-Grained Description

We have shown so far how pressure, compressibility, structural distributions (at length-scales larger than ), and free energy of the system are equivalent between UA simulations and our coarse-grained level of representation across all levels of block description. This indicates that the coarse-graining scheme is sound since all directly observable bulk phenomena are captured. We now consider two related, but non-directly observable phenomena: the system potential energy and the entropy.

We first examine the potential energy, which is readily calculated from the molecular dynamics simulation as the difference between the internal and the kinetic energy. In the coarse-grained soft sphere representation, the average potential energy is equal to the free excess energy. The entropic contributions are given by the translational entropy, which does not appear in the excess free energy, where the kinetic energy of the ideal gas is subtracted out. Thus the potential energy of the soft sphere includes only the contribution from the intermolecular pair potential, which we have previously shown to be given by

(56)

Substitution of Equation 47 into Equation 56 gives the potential energy as a function of the packing fraction, , as

(57)

which, as mentioned earlier, is a free energy in the coarse-grained coordinates.

Figure 15: Internal energy calculated from coarse grained simulation of PE200 (top) and PE1000 (bottom) as a function of the number of blocks (). The last point is the internal energy from united atom simulations. The contribution from the kinetic energy has be subtracted. The solid line depicts the theoretical prediction and is extrapolated for large .

In the multi blob coarse-grained representation, and in the united-atom representation, the monomer-level potential energy is not represented by the excess free energy, Equation 56, which corresponds to just the intermolecular contribution to the potential energy. In fact, the potential energy is a sum of intermolecular pair interactions, and intramolecular interactions, . The intramolecular contributions can be decomposed further into separate contributions: bond stretching, angle bending, torsional rotation, and non-bonded pair interactions acting between monomers on the same chain. In this way, the potential energy depends on the level of coarse-graining and is not consistent across representations.

The potential energy in the two simulations is simply calculated as the average total internal energy minus the kinetic energy contribution. Figure 15 show the potential energy per molecule for two different systems, PE200 at at and PE1000 at and as a function of the number of effective sites, . For comparison, the last data point is the potential energy from UA simulations, where the number of effective sites is equal to the number of united atoms. In both sets of simulations, the potential energy changes as the number sites is increased.

The intermolecular component to the potential energy, due to the coarse-grained potential, is generalized from the soft sphere limit as

(58)

which again gives

(59)

Equation 59 is the internal energy due to the intermolecular coarse-grained potential, which is equivalent to the excess free energy of the system in the mean field limit. However, Figure 15 shows an increase in energy with the number of coarse-grained sites. This is due to the addition of intramolecular bonding and angular potentials. Since the bond energy is a harmonic potential with a Gaussian probability distribution, the average bond energy is simply the equipartition result,

(60)

Equation 60 is equivalent to the correction term that must be included in rescaling of the dynamics of coarse-grained simulations to account for the missing entropic degrees of freedom using a freely-rotating-chain model.IVAN1 (); IVAN2 () For the angular contribution to the energy we add an additional contribution per chain. The total predicted energy, calculated from Equations 59, 60, and the angular contribution, is shown in Figure 15 with symbols representing the MS simulations and the UA simulation corresponding to the last symbol, where . The line is the representation of the analytical expression for increasing number of blobs. It is clear that the energy in the coarse-grained simulations is approaching that of the united atom simulations in the limit that .

The theoretical prediction in the long chain limit is, however, approximate because the estimated potential energy was calculated under the mean-field assumption , which becomes increasingly less accurate as the local structure becomes important. Because in the limit of , i.e. very small blocks, short range excluded volume interactions become important, and the pair distribution function presents local peaks due to packing effects, the mean-field approximation breaks down. Nonetheless, Figure 15 shows that the mean field approximation does give a quite reasonable estimate of the total energy when extrapolated to the monomer limit.

v.4 Entropy of the Variable-Lavel Coarse-Graining Model

The basic structure of any coarse-graining procedure is the averaging of the microscopic states that are then represented by effective units, with the consequence that the entropy of the system in a given coarse grained representation is modified with respect to its atomistic description, and that the extent of the change in entropy depends on the extent of the coarse-graining, i.e. on the level of details of the representation. In our model, which allows for the effortless tuning of the level of coarse-graining, it is possible to formalize the calculation of the entropy as a function of the level of details of the representation. As the entropy is a key quantity that needs to be correctly evaluated to properly use coarse-graining procedures, it has been the focus of several studies in relation to coarse-graining modeling. For example, we have argued that the change in entropy, together with the modification of the friction, provides the key correction contributions that need to be evaluated to properly reconstruct the correct atomistic dynamics from the accelerated dynamics measured in mesoscale simulations of coarse-grained systems.IVAN1 (); IVAN2 ()

A different type of entropy has been defined by Shell as the relative entropy, based on the information function that discriminates between coarse-grained configurations sampled in two levels of representationShell (). Rudzinski and Noid have discussed the relative entropy in detail in the context of numerical coarse-graining schemes such as the multiscale coarse-grained potential formalism.<