Calculation of the work of adhesion of polyisoprene on graphite by molecular dynamics simulations

ABSTRACT Elastomeric compounds are reinforced with fillers such as carbon-black and silica to improve mechanical, dynamical, and tribological properties. The stability and physical properties of these materials are dominated by the intermolecular interactions occurring at the polymer/particles interface that determine the magnitude of the polymer/particles adhesion. Using molecular dynamics simulations, in this work, we evaluate the solid–liquid interfacial tension and the corresponding work of adhesion for a system composed of graphite/Polyisoprene 100% cis-1,4 within a range of molar masses and temperatures. We employ a simulation strategy for estimating the surface tension of fluid/vacuum and fluid/solid interfaces that use directly the local stress fields in the Irving–Kirkwood formalism. Using such procedure, we decompose the stress field into the individual components of the stress tensor and correlate them with the values of the work of adhesion in the different systems analyzed. Abbreviation: PI (100% 1,4); MD: molecular dynamics; PT lateral component of the stress; PN normal component of the stress.


Introduction
Solid/liquid interfacial properties not only govern surface wettability but also drive the development of composite materials in a broad range of industrial [1,2] , biological [3][4][5] , and medical applications. [6]n particular in polymer science, polymer (nano) composites, where (nano)fillers are added into a polymer matrix, have become popular materials thanks to the fact that their properties can be tailormade depending on the filler's geometry and chemical composition. [7]n the automotive industry, the main component of car and truck tyres is the prototypical polymer composite, where the polymer matrix is an elastomer (usually cis-1-4-polyisoprene (PI), but also styrenebutadiene [9] copolymer or polybutadiene can be employed [8] ) and the fillers, added to improve mechanical and rheological properties, are normally either Carbon-Black (CB) or silica. [9]The physico-chemical interactions between the filler and the rubber at the interface play an important role on the dispersibility of the fillers in the polymer matrix and, therefore, the choice of the reinforcing fillers determines not only the final material properties but also the process conditions.The knowledge of the filler-polymer surface energy is of paramount importance to predict and understand the wettability of the filler by the rubber during the mixing process and the (re-)agglomeration (flocculation) of the fillers during the post-mixing stage [10,11] , as well as the formation of chemical linkages between the filler and the polymers during cross-linking reactions (a process commonly known as vulcanization).The enhancement of the mechanical and tribological properties exhibited by the polymer composites is mainly governed by two factors: the individual chemical identity and composition of the filler and polymer matrix, and the polymer/ filler interface that can be considered as a third component.
A useful quantity to estimate the adhesion between filler and polymer is the measure of the liquid-solid work of adhesion, W a , defined as the reversible work required to separate a solid from a liquid up to a distance at which they no longer interact.Such work can be expressed in terms of the surface tension of the individual phases and their interfacial tension as following [12] : where γ sv is the solid-vapor surface tension, γ lv is the liquid-vapor surface tension, and γ sl is the solid-liquid interfacial tension.In this study, we assume that the solid is stiff enough to consider negligible its deformation when separated from the liquid.Moreover, the liquid-vapor pressure is assumed to be low such that the amount of excess vapor adsorbed on the surface is negligible. [13,14]Under such assumptions, the solidliquid interfacial tension is equal to the solid-liquid interfacial free energy; hence in this work, we will use the term interfacial tension and interfacial free energy interchangeably.
The thermodynamic quantities that define the work of adhesion W a are also linked to the contact angle θ of the sessile drop on the solid support.
The experimental measurement of the solid-liquid surface tension and of the Wa relies on the estimation of this contact angle, but the measurements are challenging because interfacial phenomena can be hard to probe.Furthermore, such experiments are affected by the surface conditions (roughness and/or chemical heterogeneities, liquid penetration, surface deformation, underlying coated substrate) and by the effects of airborne contaminants. [15]18] Molecular modeling represents the optimal tool to investigate the solid/liquid interfacial forces and indeed it has been extensively used to evaluate the surface tension. [13,19,20]Over the years several methodologies have been proposed to estimate the interfacial free energies of a solid-liquid interface, γ sl .These methodologies can be roughly casted into a thermodynamic and a mechanical route.
In the thermodynamic route the surface tension is directly computed from the derivative of the free energy of the system.In this context, Leroy and Müller-Plathe have proposed an algorithm, called phantomwall method [19,21] , to obtain the solid-liquid surface free energy.This method has been successfully applied to study the wetting properties of Lennard-Jones (LJ) systems [21] and graphene/water interfaces [22,23] .In this approach, the solid and liquid phases are separated by a repulsive wall, which repels only the liquid.Initially the wall is within the solid substrate, and then is reversibly and perpendicularly moved toward the liquid such that it brings the liquid outside the range of interaction of the solid.The change in Gibbs free energy is associated with this transformation, and the W a is computed by thermodynamic integration. [13,21]The phantom-wall approach is particularly adapted for studying rough surfaces. [22,23]e same authors have developed a similar approach, mainly used for smooth interfaces.In their approach, called dry-surface method [24] , the solidliquid separation is realized by turning off the attractive solid-liquid interactions and maintaining unchanged the solid-solid and the liquid-liquid interactions.The change in the free energy is associated with the work required to transform the actual solid-liquid interface into a repulsive interface and separates the two phases. [24]ther thermodynamic methods have been proposed over the years.In the cleaving method [25] , the solidliquid interfacial free energy is obtained using thermodynamic integration of the reversible work per unit area required to continuously transform separate solid and melt polymer into a single system containing an interface.In the test-area method, [26] the surface tension is computed as a change in the free energy for an infinitesimal change in the surface area in the NVT ensemble.The difference in the energy between the perturbed and the reference states is computed and thermodynamically averaged over a large number of configurations, yielding to the estimation of the surface tension. [24,26,27]he thermodynamic methods, although provide accurate results, are computationally demanding because they require a large set of simulations to be performed into the integration.
Another route to calculate the surface free energy is derived from its mechanical definition, where the surface tension is obtained from the interfacial stress anisotropy.In their seminal work, Kirkwood and Buff [28] were the first who explicitly expressed the component of the stress tensor as a function of the derivative of the intermolecular potential.Such approach has been extensively employed in the computation of fluid/fluid surface tension [20,[29][30][31] , while it results more challenging and less employed in solid/liquid systems. [32,33]hese difficulties arise, for example, from the anisotropy of the solid, that leads to a different definition of the surface free energy and the surface stress tensor according to the Shuttleworth relationship [34] , valid only when the interface is free from transverse stresses.As explained above, due to the assumption that in our work the solid surface is rigid, this limitation can be overcome.Additionally, the calculation of solid/liquid surface tension is strongly affected by the parameter settings of the simulations, for example finite size effects and long range interactions [33] , as shown recently on a graphene/methane system. [35]In our case, electrostatic interactions are neglected since the system is free of charges.Thus, due to the efficiency of the method compared to the thermodynamics one, in this work we follow the mechanical route according to the formalism of Irving-Kirkwood. [28,35]This approach provides a local definition of the stress tensor σ x ð Þ from the knowledge of interaction forces and velocities of the individual atoms.
The local stress tensor σ x ð Þ is the sum of a kinetic term, σ K x ð Þ, describing the flux of momentum due to internal vibrations and a potential term, σ V x ð Þ; resulting from the intermolecular forces: where x refers to the center of a three-dimensional rectangular element of the grid in which the simulation box is discretized. [36,37]n particular, the two contributions to σ x ð Þ are defined as follows: [40] The term r αβ ¼ r α À r β is the inter-particle distance, f αβ is the intermolecular force between the particles α and β and is expressed as the sum of all the site-site forces acting between these two particles.F αβ ¼ P N β¼1 f αβ is the total force acting on particle α.The f αβ term satisfies the relation ÞÞds, the bond function weights the contribution to the stress tensor at the grid point x from the two particle interacting, α and β.B is the integral of the weight function along the line segment connecting α and β. [36,37] For a system with planar symmetry and isotropic surface, such as a planar graphite layer in contact with polymer melt, the tensor σ x ð Þ can be reduced in three components: σ xx and σ yy are the parallel vectors to the surface and σ zz is the perpendicular vector.
If we consider the stress profile along the perpendicular direction to the surface (z) we define P N z ð Þ ¼ Àσ zz z ð Þ which is the normal component of the pressure tensor along this direction, P T z ð Þcorresponds to the tangential pressure and it is given by P . Hence, we can compute the interfacial tension γ through the calculation of the stress profile [35] : The main objective of this work is to employ a reliable simulation strategy, already used for other systems [32,33,41] to calculate the surface tension and the work of adhesion of PI/graphite interfacial tension.We also investigate the effect of the temperature and of the polymer Molecular Weight (M W ) on the surface tension.We perform several simulations within a broad range of temperature (300-413 K) and molar masses (1022-13624 Da).

Methods and Computational Details
The systems studied in this work consist of a thick film of polyisoprene (PI) (100% cis-1,4) in contact with graphite on one side and to vacuum on the other as shown in Figure 1.
The PI chains are modeled using a united-atom force field, according to which each methylene and methyl group along the main chain backbone is considered as a single Lennard-Jones (LJ) interacting site, the details of the model employed can be found in ref. 9 and 47.
The thickness of the PI film, L z film , is chosen to be at least 4 times larger than the radius of gyration R g of the chains in the melt along the z direction. [42]This guarantees that far from the interfaces the structural, dynamic, and thermodynamic properties are identical to those of the bulk polymer.From previous studies, [43] we know that R g ~31 Å for cis-1,4 Polyisoprene melt with 150 monomers per chain (PI-150) and R g ~42 Å for PI-300 at T = 413 K, so we choose a thickness of L z film = 120 Å for all our systems (PI-15, PI50, and PI-200).Such a thickness allows establishing the bulk behavior in the central region of the simulation box.
The graphite surface is defined as two-dimensional planar material, made up of six graphene sheets.We use a fully atomistic representation (see Figure 1), where the carbon atoms are placed in their crystallographic structure of a hexagonal honeycomb lattice, located at z =0.The graphite slab extended for 16.75 Å along z direction.Each graphene layer has lateral dimension of L x = L y = 60 Å.All the graphene sheets are considered to be rigid, i.e., the carbon atoms are not allowed to move during the simulation runs.The distance Carbon-Carbon is fixed at 1.418 Å.
The MD simulations of the cis-1,4 PI/graphite system are performed using GROMACS (v.2016.3) [44]imulation package in the canonical ensemble with periodic boundary conditions in all the three directions.We use a leap-frog integration algorithm [45] with a time step Δt of 1fs.The temperature is kept constant by using the Nose-Hoover algorithm [46] with a time constant of 0.5 ps.
The initial polymer configurations are generated positioning randomly the chains in the simulation box; this step is then followed by a potential energy minimization.The locally relaxed configurations are then equilibrated for up to 200 ns in NVT ensemble and used as starting configurations in the subsequent MD production runs.
The interatomic interactions are modeled by the LJ 12-6 potential (cutoff radius 13 Å) and the long-range electrostatic forces are zero given that the carbon atoms charges are zero.The LJ parameters for the carbon-PI interactions are taken from ref. 47.The geometric mixing rules are used to produce the parameters of the PI-carbon interactions.
The following systems are simulated in this work: 210 chains of PI-15 on graphite surface, 54 chains of PI-50, and 15 chains of PI-200.All these systems are simulated in a range of temperatures above the glass transition temperature of PI (300, 350, and 413 K), at which the experiments are usually conducted. [47]The length of the production runs varies between 200 and 600 ns, depending on the polymer molecular weight and the temperature of the system (see Table 1).

Density Profile of PI Melt
The local density of a polymer melt in the vicinity of a planar smooth, impenetrable solid surface exhibits an oscillatory behavior, this interfacial layering has also been observed by molecular simulations with either coarse-grained [43] or atomistic models. [48]In Figure 2, the spatial arrangement of PI is quantified through the calculation of the density profiles along the normal direction to the graphite surface.In the panel a) the dependency of the density with respect to the temperature is shown, by reporting the density profile of PI-15 at three different temperatures (413 K, 350 K, 300 K).
The local density distribution of the polymer is computed by dividing the polymer region in slabs of thickness 0.35 Å along z direction and averaged over time.The curve matches the result previously obtained in ref. 43 and 52: the density profiles display three consecutive PI structured-layers close to the graphite slab.For the system PI-15 at 300 K (blue line), the first peak in the density profile is located at ~4.25 Å from the outer graphene layer with a height of 1,650 kg/m 3 , and the second peak is located at 9.35 Å with an intensity of about 1,033 kg/m 3 , and the last one at 14.5 Å. [47] As expected at the polymer/vacuum interface the melt density drops to zero, exhibiting the characteristic sigmoidal shape typical of a melt free surface. [42]By increasing the temperature, the polymer/vacuum interface drops to zero less steeply.The slope of the linear fit of the density profile curve at the liquid-vapor interface increases from −108.7 at 300 K to −73.3 at 413 K.This is known in the literature and is connected to the decrease in the surface tension of the melt with increasing the temperature. [49,50]he effect of the graphite slab on the distribution of PI is propagated until about 19 Å from the outer graphene layer in the direction normal to the surface.The density at the central region of the film converges to the value of the density of the bulk PI melts, which is experimentally about 840 kg/m 3 at 413 K and 910 kg/ m 3 at 300 K. [47,51] From Figure 2(a), we observe that the same oscillatory behavior in the interfacial region is present for the system PI-15 at different temperatures (three peaks in density at the same position), but the intensity of the peak decreases with increasing temperature.It can be noticed that, as the temperature increases, the bulk value of the density in the middle regions of the film decreases following the same trend than a bulk PI melt (see Table 2).
Figure 2(b) shows also results for the density profile of PI/Graphite as a function of molecular weight (M W in Da), at T = 413 K.We observe that at low molar mass, the bulk density, ρ bulk , of PI melts increases slightly with the chain length from ρ bulk , PI-15 = 828 kg/m 3 to ρ bulk , PI-200 = 869 kg/m 3 (see Table 2) while for larger molar masses the density saturates to a constant value as it shows in Figure 3.
In Figure 3, we report the results for the density of PI in the bulk as a function of M W at three different temperatures.The open symbols in Figure 3, at low M W (<2,000 Da) for T = 413 K and 350 K are extracted from Ref 52.We observe that the density increases quickly with M W , particularly at low M W (number of monomers < 15) until it reaches the plateau value at high M W.
In Figure 3 each curve is fitted by using the equation below [53] : where ρ T; 1 ð Þcorresponds to the value of the density at infinite M W and V e T ð Þ is the excess free volume of chain ends.By fitting the density data to Eq. 6 we can estimate the temperature dependence of ρ T; 1 ð Þ and V e T ð Þ.Both parameters are function of temperature only, and they depend linearly on T.
The fitting parameters a, b, c, and d in Eqs.7 and 8 result in very good agreement with experimental [54] and computational [48,52] data of previous studies.

Profile of the Tangential and Normal Stress Components
The calculation of the interfacial tension of PI interacting with the graphite substrate is calculated by using the mechanical route through the IK method, as described in the previous section.
To compute the local stress tensor we use the GROMACS-LS code of Vanegas et al. [37,[55][56][57] Through this code we obtain P T and P N and we define the lateral   [52] .
pressure profile as π z ð Þ ¼ P T z ð Þ À P N z ð Þ [58,59,60] In order to compute the profiles of the normal and the tangential stress, we discretize the simulation box into a threedimensional rectangular grid of cell size L grid ¼ 1Å and compute the average stress over each cell of the grid.
In Figure 4 on the left axis, we show the lateral pressure profile π(z) which arises from the local forces acting on the polymer in the direction of the graphene plane for the system of PI-15 at 413 K.
At equilibrium, due to mechanical stability, the integrated lateral pressure profile in the bulk is zero.By comparing the stress profile and the density profile on the right axis of Figure 4 we see a perfect correspondence between the peaks of π(z) and the density profile, located, respectively, at z = 4:25 Å and 9.25 Å from the graphite surface.The local stress oscillations die out within ~19 Å of the graphite walls in both cases.
Figure 5 shows the profile of interfacial tension obtained by IK method.The black solid line displays ðP N z ð Þ À P T z ð ÞÞΔz and the red dotted line the profile of the integral γ(z) (according to Eq. 5).Statistical fluctuations of interfacial tensions are estimated using the block averages, i.e. the calculation of the interfacial tension is performed over 200,000 configurations.Standard deviations of the interfacial tensions are calculated by breaking the trajectories into 4 block averages.
The value of the solid-liquid interfacial tension, γ sl , for the system PI-15 at 413 K simulated here is 25.4 mN/m, in agreement with experiments of CB/rubber where γ sl is in the range of 19-30 mN/m for different CB structures. [8]n order to understand what is the relative contribution of the intermolecular forces to the surface tension value, we decompose the stress profile for the system PI-15 at T = 413 K.In Figure 6, the components to the tangential stress, σ xx (= σ yy for symmetry), which are the main contributions to π(z), are reported.The components include the kinetic and potential stress parts with the latter containing the contributions from pairwise non-bonded interactions (van der Waals), bond stretching, bond angle, improper, and proper dihedral (Ryckaert-Bellemans potential).From Figure 6, it appears clear that among these contributions, the improper dihedral component (red solid line), which is defined by the atoms connected through the double bond, displays the largest positive and negative values.This contribution nearly coincides with the total lateral stress profile and therefore, the other parts nearly balance each other in this region.In Figure 7, we compare σ xx from the dihedral contribution of the system PI-15 at different temperatures.We observe that by increasing the temperature, the intensity of the first peak decreases by 10% at 350 K and 17% at 413 K with respect to the highest peak at 300 K. Interestingly, we notice that at the liquid-vapor interface the peaks are shifted to the highest z values and the intensity decreases by increasing the temperature.This behavior can be correlated with the density profile at the liquid-vapor interface where the structure of the polymer is expanded at high temperatures.This result confirms that the surface tension is directly linked with the structural properties of the polymer.
In Table 3, the values of the interfacial tension for different systems at T = 413 K are reported.
From Table 3, we can see that by increasing the M W , the liquid-vapor surface tension at a fixed temperature increases until it reaches a plateau at 3,406 Da (PI-50).The solid-liquid surface tension instead slightly decreases with increasing the M W until a plateau value, again at 3,406 Da.This mild dependency on the M W is also in agreement with the structural properties of coarse-grained and atomistic models of polyisoprene/graphite [43,61] systems where it has been observed that the average length of trains (i.e. the length of the polymer chain directly adsorbed on the surface) ranges from 5 to 6 monomers independently from the M W . [43,61] This indicates that the number of the adsorbed segments does not change with the M W and therefore, since the main contribution to the W a is due to the dihedrals, which are strictly linked with the conformation of the adsorbed chains, we expect a similar trend for the solid-,liquid surface tension.
In Tables 4 and 5, we can infer the effect of temperature on the surface tension for the systems PI-15 and PI-50 (the results obtained for the longest polymer chain of 200 monomers were affected by a high statistical error due to the long relaxation time required at low temperatures and therefore we omit them from the table).We observe that for both systems, by increasing the temperature, the surface tension solid-liquid and liquid-vapor decreases in agreement with what has been observed for polyethylene. [41]This behavior reflects a decrease in the cohesive energy density of the polymer.
The computation of the local stress profile allows for the calculation of the interfacial free energies via Eq. 5, which can provide a deep insight into the wetting properties such as the W a .By knowing the values of γ sl, γ lv and γ sv we can calculate the value of the W a according to Eq. 1, reported in Table 6.The value of γ sv has been recently measured by Perkin et al. [15] for a single layer and few layers of graphene as 115 and 119 mN/m, respectively.By using the value 115 mN/m and the calculated γ sl and γ lv , the values of W a are computed for PI-15 and PI-50 at three different temperatures (413 K, 350 K, and 300 K) and for PI-200 at 413 K (see Table 6).From the values reported in Table 6, we notice that the W a decreases with increasing the temperature, while there is a negligible dependency with respect to the M w.

Summary and Conclusions
In the current work, we employ and optimize a simulation strategy for the estimation of interfacial energies and wetting properties of polymer/solid interfaces.By performing  MD simulations, we calculate the structure and surface tension of the solid-liquid graphite-polyisoprene interface.
From the analysis of the density profile, we observe an oscillatory behavior in the region adjacent to the graphite.By fitting the data of density vs M W we are able to estimate the volume of chain free ends from the free volume of the system and the density at infinite M W which are in good agreement with experimental and computational data. [48,54]he lateral pressure profile offers an appealing way to gain an insight into the relationship between polymer chemical structure and its interfacial thermodynamics properties.We observe a perfect match in the peaks of the lateral stress and density profiles.By analyzing the different contributions to the stress profile, we find that the main influence on the tangential stress is due to the dihedral components.
Furthermore, we estimate the solid-liquid and liquidvapor interfacial tension for the systems composed of 15 and 50 monomers at three different temperatures (413 K, 350 K, and 300 K) and we find that the surface tension decreases with increasing temperature.By increasing the molar mass, the liquid-vapor and solid-liquid surface tensions reach a plateau already for chain composed by 50 monomers.Finally, the values of the surface tensions allow for the calculation of the W a which decreases with increasing the temperature, instead there is a negligible dependency with the M w .

Figure 1 .
Figure1.MD snapshot of a thin PI melt film confined between a semi-infinite graphite (six layers) phase on its one side (bottom) and vacuum on its other side (top).The z-direction is normal to the surface.Periodic boundary conditions apply only in the x and y directions of the coordinate system, the surface area is defined by A = L x Â L y (with L x = L z = 51 Å and L z = 60 Å).The different colors distinguish the several chains.

Figure 2 .
Figure 2. Density profile of PI as a function of the distance from the outer graphene layer of the graphite slab (located at z= 0) at various chain lengths and temperatures.(a) Shows results for PI-15 at several temperatures at 413 K (green dotted line), 350 K (red dashed line), 300 K (blue line).The inset shows the first peak in the density profile.(b) Results from PI-15 (black line), PI-50 (red dotted line), and PI-200 (blue dashed line) at 413 K.

Figure 3 .
Figure 3.The density of the cis-1,4 PI melt, in the bulk, as a function of molecular weight.The open symbols at low molecular weights (M W < 2000 Da) are taken for T = 315 and 413 K from ref.[52] .

Figure 4 .
Figure 4. Later pressure profile (π) (blue circles) and density profile (red solid line) of PI melt a function of z the normal direction to the solid substrate for the system PI-15 at 413 K.

Figure 5 .
Figure 5. Profiles of ðP N z ð Þ À P T z ð ÞÞΔz (black solid line) and of its integral (red dashed line) as a function of z normal direction to the graphite surface for PI-15 at 413 K.

Figure 6 .
Figure 6.The different contributions of angles (black), bonds (dashed blue line), van der Waals interactions (orange dotted line), dihedrals (magenta solid line), improper dihedral (red solid line), and the kinetic component (green solid line) on the tangential stress profile (σ xx ¼ σ yy ) for the system of PI-15 at T = 413 K.

Figure 7 .
Figure 7. Improper dihedral contribution on the tangential stress profile (σ xx ¼ σ yy ) for the system of PI-15 at T = 413 K (dotted green line), 350 K (dashed red line), and 300 K (solid blue line).The insert shows the first peak of the improper dihedral component of the tangential stress profile.

Table 1 .
Computational details of the systems simulated in this work.

Table 2 .
Values of ρ bulk at various temperatures and molecular weight.

Table 3 .
Surface energy values of systems with different chain length for the PI at same T = 413 K.

Table 4 .
Interfacial tensions of PI-15 at different temperatures.

Table 5 .
Interfacial tensions of PI-50 at different temperatures.