Modeling elastic properties of polystyrene through coarse-grained molecular dynamics simulations

TL;DR

This study uses extended coarse-grained molecular dynamics simulations with the MARTINI force field to analyze and tune the elastic properties of polystyrene, aligning simulation results with experimental data.

Contribution

It introduces modifications to the MARTINI force field to better replicate the elastic properties of polystyrene in molecular dynamics simulations.

Findings

Simulations successfully reproduce experimental elastic moduli of polystyrene.

Iterative force field modifications improve shear and bulk modulus accuracy.

Extended $s$-scale simulations capture slow relaxation processes.

Abstract

This paper presents an extended coarse-grained investigation of the elastic properties of polystyrene. In particular, we employ the well-known MARTINI force field and its modifications to perform extended molecular dynamics simulations at the $\mu$s timescale, which take slow relaxation processes of polystyrene into account, such that the simulations permit analyzing the bulk modulus, the shear modulus, and the Poisson ratio. We show that through the iterative modification of MARTINI force field parameters it turns out to be possible to affect the shear modulus and the bulk modulus of the system, making them closer to those values reported in the experiment.

Tables

Table 1. Table 1: Parameters of bonded, angular and nonbonded interactions employed in the coarse-grained simulations. Notations of the beads are consistent with rossi2011coarse such that B and R indicate the backbone and ring, respectively. The parameters for the standard MARTINI force field for polystyrene have been adopted from an earlier study rossi2011coarse and are denoted as MARTINI-std in the table and text. Two modifications to the force field are also considered where the non-bonded energy parameter ε 𝜀 \varepsilon of the standard force field is scaled by 1.5 and 2; the corresponding force fields are denoted as MARTINI-1.5 and MARTINI-2, respectively.

Bonded	$r_0$ (Å)	$k_0$ (kcal/mol/Å2)
B–R	2.7	9.5602
R–R	2.7	constrained
Angular	${\vartheta }_0$ (degrees)	$k_{\vartheta }$ (kcal/mol)
B–R–B	52	65.7266
B–B–R	120	2.9876
B–R–R	120	11.9503
R–R–R	60	0.0
Nonbonded pairwise	$r_{min}$ (Å)	$\epsilon$ (kcal/mol)
MARTINI-std
B $\cdots$ B	4.82659	-0.6274
R $\cdots$ R	4.60209	-0.5736
R $\cdots$ B	4.82659	-0.5557
MARTINI-1.5
B $\cdots$ B	4.82659	-0.4183
R $\cdots$ R	4.60209	-0.3824
R $\cdots$ B	4.82659	-0.3705
MARTINI-2
B $\cdots$ B	4.82659	-0.3137
R $\cdots$ R	4.60209	-0.2868
R $\cdots$ B	4.82659	-0.2778

Table 2. Table 2: A summary of the performed coarse-grained simulations reflecting two compression regimes (volumetric and axial) studied. Extensive equilibration was also performed prior to the production simulation, and is summarized in the table.

Simulation type	Box size	Simulation	Temperature	Force fields
	(Å$\times$Å$\times$Å)	length (ns)	(K)	employed
equilibration	167.2$\times$167.2$\times$167.2	30	300	MARTINI-std
equilibration	167.2$\times$167.2$\times$167.2	100	1000	MARTINI-std
equilibration	167.2$\times$167.2$\times$167.2	100	800	MARTINI-std
equilibration	167.2$\times$167.2$\times$167.2	100	600	MARTINI-std
equilibration	167.2$\times$167.2$\times$167.2	100	400	MARTINI-std
equilibration	167.2$\times$167.2$\times$167.2	2800	300	MARTINI-std
volumetric	171.0$\times$171.0$\times$171.0	800	300	MARTINI-1.5,2
volumetric	171.5$\times$171.5$\times$171.5	800	300	MARTINI-std,1.5,2
volumetric	172.0$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2
volumetric	172.5$\times$172.5$\times$172.5	800	300	MARTINI-std,1.5,2
volumetric	173.0$\times$173.0$\times$173.0	800	300	MARTINI-std,1.5,2
volumetric	173.5$\times$173.5$\times$173.5	800	300	MARTINI-std,1.5,2
volumetric	174.0$\times$174.0$\times$174.0	800	300	MARTINI-std,1.5,2
volumetric	174.5$\times$174.5$\times$174.5	800	300	MARTINI-std
axial	170.5$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2
axial	171.0$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2
axial	171.5$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2
axial	172.0$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2
axial	172.5$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2
axial	173.0$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2
axial	173.5$\times$172.0$\times$172.0	800	300	MARTINI-std,1.5,2

Equations

$$m_i{\vec{a}}_i=m_i\frac{d^2{\vec{r}}_i}{d{t}^2}={\vec{F}}_i-\gamma {\vec{v}}_i+\sqrt{2k_{\mathrm{B}}T\gamma }{\vec{R}}_i(t),\forall i=1\dots N_{beads}.$$

$$U$$

$$L={\left(\frac{M_{sty}{N}_0{N}_{str}}{\rho }\right)}^{1/3},$$

$$⟨{\rho }_{cell}⟩=\frac{1}{N_{conf}}\sum _{t=1}^{Nconf}\sum _{i\in cell}{m}_i/V_{cell},$$

$$P(V)=\frac{3K}{2}\left[{\left(\frac{V_0}{V}\right)}^{7/3}-{\left(\frac{V_0}{V}\right)}^{5/3}\right]\left\{1+\frac{3}{4}(K^{\prime }-4)\left[{\left(\frac{V_0}{V}\right)}^{2/3}-1\right]\right\}$$

$$E(V)=E_0+\frac{9V_{0}K}{16}\left\{{\left[{\left(\frac{V_0}{V}\right)}^{2/3}-1\right]}^3{K}^{\prime }+{\left[{\left(\frac{V_0}{V}\right)}^{2/3}-1\right]}^2\left[6-4{\left(\frac{V_0}{V}\right)}^{2/3}\right]\right\}.$$

$$P=\frac{N{k}_{B}T}{V}+\frac{1}{3V}⟨\sum _{i=1}^N{\vec{r}}_i{\vec{F}}_i^{int}⟩.$$

$$P=\frac{N{k}_{B}T}{V}-\frac{1}{3}\mathrm{Tr}\hat{\sigma },$$

$${\sigma }_{ij}=\sum _{kl}{C}_{ijkl}{\epsilon }_{kl},$$

$$C_{ijkl}=(K+\frac{2}{3}G){\delta }_{ij}{\delta }_{kl}+G\left({\delta }_{il}{\delta }_{jk}+{\delta }_{ik}{\delta }_{jl}\right),$$

$$\frac{1}{3}\mathrm{Tr}\hat{\sigma }=K\mathrm{Tr}\hat{\epsilon },$$

$$\mathrm{volumetric}\mathrm{stress}$$

$$\mathrm{volumetric}\mathrm{strain}$$

$$\frac{V-V_0}{V_0}=\frac{L_x-L_0}{L_0}.$$

$${\sigma }_{\alpha \beta }=(P_0{)}_{\alpha \beta }-P_{\alpha \beta },$$

$${\sigma }_{xx}=\left(\frac{4}{3}G+K\right){\epsilon }_{xx},{\sigma }_{yy}={\sigma }_{zz}=\left(-\frac{2}{3}G+K\right){\epsilon }_{xx}.$$

$$\frac{2{\sigma }_{xx}-{\sigma }_{yy}-{\sigma }_{zz}}{4}=G{\epsilon }_{xx}.$$

$$\nu =\frac{3K-2G}{6K+2G}.$$

Full text

I Introduction

Polymers and polymer-based materials are widely used as components of various products and goods nowadays and are employed in daily life. Such an extensive use is mainly due to the wide range of physical and mechanical properties of polymers, which can be enhanced by numerous additives camargo2009nanocomposites ; Mark1999 ; Mark2013 ; Harper2002 . Modern technologies have made it possible to create new types of polymer composites using nanoparticles of various nature and composition fawaz2015synthesis . The manufacturing process is, however, often associated with specific technological difficulties, including the monitoring of the aggregation-disaggregation processes to obtain the uniform distribution of nanoparticles in the polymer matrix supova2011effect .

Preparation, characterization, and applications of polymeric nanocomposites with various nanoparticles has influenced multiple research areas nasir2015review ; armstrong2015introduction . Nanocomposites with silica nanoparticles have particularly attracted substantial academic and industrial interest zou2008polymer ; mathioudakis2016multiscale ; rahman2012synthesis ; polymer/silica composites are seemingly the most commonly reported ones, which have also been successfully employed in a variety of applications, including tire industry, food industry, civil engineering, and recyclable materials Mark2013 ; Harper2002 .

Different polymers have been used for polymer/silica composites. In particular polypropylene bracho2012functionalization , polyimide musto2004polyimide , polyamide sarwar2008polyamide , poly(ethylene terephthalate) liu2004preparation , polymethylmethacrylate yang2004pmma , polyurethane petrovic2000structure have been employed, while polystyrene was also a popular focus of both experimental liu2005thermal ; vaziri2011thermophysical ; bartholome2005viscoelastic and theoretical mathioudakis2016multiscale ; ndoro2011interface ; ghanbari2011interphase research. Polystyrene is a general-purpose, well-known frequently utilized, light weight, famous plastic with worthy dimensional stability, noble chemical resistance, easy processing, and low cost Mark1999 . It has been extensively utilized in packing materials, electronics, household appliances, etc.

Nowadays, computer simulation becomes a promising and powerful tool to investigate the structure and properties of polymers, thereby providing an in-depth understanding of experimental data. Historically, polymers have been widely studied through theoretical continuous models Rivlin1948 ; Mooney1940 ; LL7 as well as discreet microscopic models WallFlory1951 ; JamesGuth1943ThrElstPropRub , that permit establishing bulk properties of the materials. The well-known hyperelasticity model of Mooney and Rivlin (MR) Rivlin1948 ; Mooney1940 , and first microscopical models for polymer networks elasticity were suggested in the 50th. Since then generalization of the MR model was suggested, but the first successful microscopical / molecular models, that take entanglements effects into account, are available since 80th Edwards1967StatMechTopologyI ; Gaylord1982 ; MergellEveraers2001TubeModels ; RubinPanyuk2002ElastPolymNetw . Molecular dynamics (MD) was also seen as an important tool to study polymer properties. The most popular approach exploited in MD simulations of polymer dynamics is the generic bead-spring model of Kremer and Grest (KG) GrestKremer1986MDSimPolymHeatBath , originally invented for polymer melts. The KG model is actively applied until nowadays to study the dynamics of melts KremerGrest1990DynEntanglLinPolyMelt ; SvaneborgKarimiHojdisFleckEveraers2016MultiscaleApproach , solid rubbery polymers GrestKremer1990StatPropRandXLnkdRubb ; EveraersKremer1996ElastPropPolyNetw ; SvaneborgEveraersGrestCurro2008StressContributions and glassy polymers after proper modification Bennemann1998 .

Recently, due to the emerging supercomputer facilities, it also became feasible to study polymers through advanced MD simulations, including atomistic MD and unprecedentedly long simulation on experimentally accessible timescales. Atomistic MD has been widely applied to study biological polymers, such as, e.g., proteins Dubey2018 ; Kimo2018 ; Kattnig2018a ; Akimov2018 ; Friis2017 and DNA Klecka2018 ; Jepsen2017a ; Zou2012a , and additionally it was also applied to synthetic polymers, such as, e.g. polystyrene rossi2011coarse ; Chakravarthy2000 . Advanced MD simulation protocols have been developed for both explicit atom (EA) models, treating every atom as a separate interaction site Yoon1993 , but also united atom (UA) models where several atoms are grouped into a single bead (e.g., CH2 or CH3) barrat2010molecular . The reduction of the number of interaction sites in the UA model leads to significant simplifications of the complex polymeric systems, and, therefore, permits performing significantly longer simulations capable of predicting many universal polymer properties independent of local chemical structure.

UA model has been successfully employed for studying polystyrene lyulin2002correlated ; lyulin2002molecular ; lyulin2003molecular ; lyulin2005strain ; lyulin2007time . In these coarse-grained simulations the typical simulation time was about 2 ns, and it was observed that polystyrene chains experience several relaxation processes with very different relaxation times lyulin2007time . Some of the established relaxation times turned out to be comparable to the overall simulation time, which indicates the necessity of significantly longer simulations. There are alternatives to the UA approach for studying polystyrene on a coarse-grained level, which includes the so-called 2:1 model harmandaris2009dynamics and the MARTINI force field Marrink2007 ; Marrink2013 ; rossi2011coarse .

The paper presents extended coarse-grained MD simulations of polystyrene and aims to establish its various elastic properties, such as the bulk modulus, the shear modulus, and the Poisson ratio computationally. The simulations rely on the well-established MARTINI force field Marrink2007 ; Marrink2013 ; rossi2011coarse , which is revealed to be somewhat inaccurate in predicting all the first order elastic moduli of polystyrene simultaneously. Investigations reported here demonstrate that through iterative modification of the force field parameters it turns out to be possible to affect the shear modulus and the bulk modulus, making them closer to those values reported in the experiment. Through these modifications, it is revealed that the MARTINI approximation (and its modifications) can provide a qualitative force field for computing of the first-order elastic constants. The key aim of this study is to provide a model and a force field that have been justified to be qualitatively applicable for studying pure polystyrene, such that this model could be further extended through doping with nanocomposites such as, e.g. the SiO2 nanoparticles.

The paper delivers a general recipe on how to use MD to obtain elastic properties of polymers. It discusses several theoretical approaches that can be used to determine the elastic moduli and argues about their accuracy and applicability. Specifically, the analysis of potential energy and internal pressure of the polystyrene under strain has been analyzed, and it has been demonstrated that both approaches may yield significantly different results in terms of the elastic moduli values and polystyrene relaxation time.

II Theoretical methods

This section discusses the key theoretical and computational methods used for studying the elastic properties of polystyrene. Below we describe the general approach of the computational method used in the investigation, explain how the simulations were set up in the coarse-grained representation and how the interactions between coarse-grained beads were modeled. Finally, a table summarizing all the performed simulations is presented.

II.1 System preparation: all atom representation

The bulk polystyrene studied here was initially constructed in the all-atom representation. The system was assumed to contain 216 polystyrene strings consisting of 120 monomers each. The strings were placed in a simulation box of 167.2 Å$\times$167.2 Å$\times$167.2 Å to ensure the volumetric density of the system being equal to 0.96 g/cm3, being consistent with the literature Hocker1971 ; Zoller1995 ; rossi2011coarse . Initially, styrene monomers were placed randomly from the two sides of the parent polymer chain, as illustrated in Fig. 1A. The studied model assumed several short polystyrene strings instead of a single long chain to ensure faster relaxation of the system. The studied model is consistent with earlier investigations rossi2011coarse ; Mathioudakis2016 , where bulk properties of polystyrene were probed on a system consisting of multiple shorter polymeric chains.

Atomistic MD simulations of the constructed bulk polymer were performed to ensure initial equilibration of the system prior to coarse-graining. MD program NAMD 2.12 PHIL2005 was used in the computations, while interactions between the atoms of the styrene monomers were modeled with the CGenFF force field Vanommeslaeghe2010 ; the partial charges of one monomer are indicated in Fig. 1B. Analysis of the simulations results and visualization of molecular structures were accomplished with VMD 1.9.2 HUMP96 . Simulations were carried out with the 1 fs integration time step. The cut-off distances for the van der Waals and Coulomb interactions were set to 12 Å , where the long-range electrostatic interactions were treated using the PME method darden1993particle ; Solovyov2012c , employing periodic boundary conditions. The NVT statistical ensemble was used for the simulation, with a temperature of 300 K by applying Langevin forces with a damping coefficient of 5 ps*-1*. The system was simulated in all-atom representation for 100 ns to ensure a random and uniform polystyrene configuration prior coarse-graining, while the elastic properties were computed within the coarse-grained approximation. This study aims to fine-tune the parameters of the coarse-grained potential in such a way that it can reproduce experimental observables directly, i.e., without transferring the coarse-grained model to an all-atom one, as done in some other studies ZhangMoreiraStuehnDaoulasKremer2014EquilHierarchStrat ; SvaneborgKarimiHojdisFleckEveraers2016MultiscaleApproach .

II.2 Coarse-graining the system

The polystyrene chains were coarse-grained to permit longer simulations and ensure equilibration of the system. Each styrene monomer was substituted with four coarse-grained beads, as illustrated in Fig. 2. Such a representation is consistent with earlier studies rossi2011coarse ; the backbone carbon atoms of the polystyrene are modeled through one bead of type B, while three beads of type R represent the styrene side chain.

The coarse-grained system of bulk polystyrene was constructed from the last frame of the atomistic MD simulation, such that each of the atomistic polystyrene chains was represented by a coarse-grained analog, using an in-house conversion script.

II.3 Coarse-grained molecular dynamics and MARTINI force field

The coarse-grained bulk polystyrene was simulated dynamically employing NAMD 2.12 PHIL2005 in the NVT statistical ensemble. In production simulations, i.e., those that were used for data sampling, the temperature $T$ was assumed 300 K, the number of coarse-grained particles was 103,680, while the volume of the simulation box was different, depending on the simulation; all performed simulations are listed in the summarizing Table 2. Prior to the production simulations, an extensive equilibration simulation of the system was carried out to ensure the uniform density distribution of the polymers in the system. The temperature was gradually changed in the equilibration simulations from 1000 K to 300 K as indicated in Table 2 to anneal the system to thermal equilibrium.

The equation of motion in the case of the coarse-grained dynamics is the same as in the case of atomistic simulations and reads

[TABLE]

Here ${\vec{F}}_{i}$ is the force acting on a coarse-grained bead $i$, $k_{\mathrm{B}}T$ denotes the thermal energy in the system, $\gamma$ is the damping coefficient that is assumed equal 1 ps*-1* in all coarse-grained simulations and ${\vec{R}}_{i}(t)$ represents a delta-correlated stationary Gaussian process with zero mean. Equation (1) is called as the Langevin equation Solovyov2012c ; PHIL2005 and is used to control the temperature in the system. The difference from the atomistic simulations is that it is applied to every bead instead of every atom in the system. Due to the significantly higher mass of one bead, Eq. (1) can be solved efficiently with a higher integration time step, as compared to the atomistic simulations; a time step of 15 fs was used in the present investigation.

In the performed simulations, the forces acting between the beads are determined by the parametric MARTINI potential Marrink2007 ; Marrink2013 ; rossi2011coarse , which reads

[TABLE]

This simple parametrization mimics the non-bonded van der Waals interactions between the beads through the Lennard-Jones potential and the covalent interactions through a number of parameters and fit functions. The first term in Eq. (2) is the non-bonded term. The second and the third terms describe the potential energy arising due to stretching of bonds between pairs of bound beads, and the change of angles between every topologically defined triple of beads in the system. The force field parameters $r_{0}$, $k_{0}$, $\vartheta_{0}$, $k_{\vartheta}$. $r_{min}$, and $\varepsilon$ for styrene monomers are available from the literature rossi2011coarse and, therefore, have also been utilized in the present study. For convenience, Table 1 lists these parameters for all the relevant pairs and triples of the coarse-grained beads R and B, see Fig. 2. The table includes the parameters for the standard MARTINI potential (MARTINI-std), while at the same time it also summarizes two sets of non-bonded parameters (MARTINI-1.5 and MARTINI-2), which have been employed in the present investigation to fine-tune the interactions.

II.4 Simulation protocol

Equations (1)-(2) were solved numerically, where the simulation protocol was following the standard guidelines advised by NAMD PHIL2005 . For the coarse-grained simulations, the cutoff interaction distance was set to 12 Å, the switch distance to 9.0 Å, the pair list distance 14.0 Å, and the dielectric constant was assumed equal to 15.

A series of simulations were carried out to determine the elastic properties of the bulk polystyrene and to establish the optimal interactions within the system; the complete list of simulations is summarized in Table 2. Simulations were performed for different sizes of the simulation box, as indicated in the table, and with different force fields for treating the non-bonded interactions, see Table 1. Two compression regimes were considered: uniform volumetric compression and axial compression. These perturbations of the bulk polystyrene were modeled by changing the size of the simulation box while keeping the same number of particles inside. The resulting changes in energy and internal pressure were then recorded and analyzed, as discussed in the following sections.

III Results and Discussion

This section presents and discusses the results of the simulations. Bulk polystyrene is first equilibrated in a cubic simulation box to ensure a uniform density distribution. Elastic properties are computed next for the equilibrated system for the different variants of the MARTINI coarse-grained potential, studied here. The obtained results are compared with earlier investigations obtained both experimentally and computationally.

III.1 Equilibration of polystyrene

The density of polystyrene varies depending on its type and manufacturing technique. A characteristic value is expected to be around 960 kg/m3 Hocker1971 ; Zoller1995 ; rossi2011coarse , which is also assumed in the present investigation. The edge length of the simulation box can be determined from the expected density value $\rho$ as:

[TABLE]

where $M_{sty}=8(m_{H}+m_{C})$ is the mass of one styrene monomer $N_{0}=120$ is the number of monomers in one polystyrene string, while $N_{str}=216$ is the number of strings in the simulation box. For the desired value of the polystyrene density, the size of the cubic simulation box becomes 167.2$\times$167.2$\times$167.2 Å3, see Table. 2. Extensive equilibration has been performed following a protocol that includes temperature annealing, as outlined in the table. The equilibration simulation was 3.23 $\mu$s long and allowed to ensure a uniform density distribution of bulk polystyrene. Figure 3 shows a snapshot of the system after the equilibration and illustrates how periodic boundary conditions were effectively used to eliminate the possible boundary artifacts. In other words, the simulated system represents an infinite sample of polystyrene that has been severely intertwined following the annealing procedure. Indeed, Fig. 3 features an utterly random orientation of the polystyrene monomers without any preferred orientation.

Figure 4 shows the average density of the bulk polystyrene inside the simulation box. Here the simulation box was subdivided into 6$\times$6$\times$6 cells, and the average density of the system was computed in each of these cells following the definition:

[TABLE]

where $m_{i}$ is the mass of an atom inside a cell (which is calculated from the corresponding coarse-grained bead), while $V_{cell}$ is the volume of the cell. $N_{conf}$ is the number of snapshots of the simulation box used for averaging. In the present calculations we have assumed $N_{conf}=20$.

Figure 4 features a three-dimensional density distribution that does not show any specific regions with significantly different density. The distribution is rather uniform and illustrates that small local fluctuations of density are possible due to the small and finite size of the sampling cells. Figure 4 demonstrates that the extended equilibration was sufficient to achieve a uniform distribution of the system, making it appropriate for the study of elastic properties.

III.2 Volumetric compression

The first compression regime to consider is the volumetric compression of the bulk polystyrene. In this case, the equilibrated system is subjected to a uniform compression applied to all sides of the simulation box. This is achieved by varying the size of the simulation box, which effectively results in the change of the polystyrene density, and internal stress, while the performed change in the box size is easily related to strain. As the applied strain is finite, and the compression of the system is hydrostatic, it is possible to employ the Birch-Murnaghan isothermal equation of state Birch1947 ; Murnaghan1944 to describe the relationship between the volume of the simulation box and the internal pressure to which it is subject. Specifically, the third-order Birch–Murnaghan isothermal equation of state can be used and reads as:

[TABLE]

where $P$ is the pressure, $V_{0}$ is the reference volume of the simulation box, $V$ is the deformed volume, $K$ is the bulk modulus, and $K^{\prime}$ is the derivative of the bulk modulus with respect to pressure. Equation (5) allows to compute the internal (potential) energy of bulk polystyrene as the function of volume, by integration:

[TABLE]

Here $E_{0}$ is the integration constant which has the meaning of the reference energy of the polystyrene sample. This equation of state can be used to determine the bulk modulus $K$ if the change of the internal energy of the system is known as the function of its volume, which in turn can be established from simulations.

Figure 5A-C shows the time evolution of the internal (potential) energy of the entire simulation box with polystyrene computed for the three different force fields employed, see Table 2. Figures 5A-C were obtained for the different sizes of the simulation box, as indicated with color and feature evolution of the internal energy of the system during the interval of 800 ns, which has been chosen to (i) allow polystyrene to adjust to the new box size, i.e. re-equilibrate, and (ii) sample the energy in order to acquire the $E(V)$ dependence (production simulation). The production simulations were assumed to last 375 ns in all cases.

The average energies from Figs. 5A-C, sampled over the production simulations, allow plotting the relative potential energy of bulk polystyrene computed upon volumetric compression, see Fig. 6A. The results of the simulations are shown with symbols, while solid lines correspond to the theoretical model of Birch and Murnaghan, Eq. (6). Fitting the simulation points with Eq. (6) allows establishing the values of the bulk moduli for the three force fields considered, as shown in the inset. The figure shows that the simulated data matches the theoretical curve closely and reveals that the bulk modulus of polystyrene decreases if the van der Waals interaction between the coarse-grained beads is decreased, i.e., the value of $K$ appears lower for the MARTINI-2 force field, as compared to the value obtained in MARTINI-std. This happens because the stiffness of the polystyrene chains in the considered model is controlled through the parameter $\varepsilon$ of the MARTINI force field equation, Eq. (2). This parameter is steadily decreased, see Table 1, for the MARTINI-std$\rightarrow$MARTINI-1.5$\rightarrow$MARTINI-2 force fields shifts, thereby making bulk polystyrene in the case of MARTINI-2 have a lower bulk modulus than in the case of MARTINI-std force field. Note that the value of $K$ obtained for MARTINI-1.5 force field from fitting the average energies with Eq. (6) appears to be lower than in the case of MARTINI-2. This can be attributed to the significant uncertainty of energy, as follows from Fig. 5, which makes the obtained volumetric stress values obtained in Fig. 6A rather qualitative.

It is important to stress that the density of the bulk polystyrene used in the computational study of volumetric compression appears to be lower than the equilibrium value of 960 kg/m3. This is evident from Table 2 and is also featured in Fig. 6A, where the minimum of the bulk polystyrene potential energy appears at a volume that is somewhat lower than 167.2$\times$167.2$\times$167.2 Å3. The approximate density of bulk polystyrene at the minimum of the potential energy is 880 kg/m3, which is still an acceptable value for polystyrene Zoller1995 ; rossi2011coarse . The reason for this lowering of density is because the original long-equilibration was performed at a fixed polystyrene density of 960 kg/m3, which, however, does not match the configuration with the minimum potential energy for the employed force fields. Once the simulation box size changes, it is, therefore, important to perform an additional equilibration to ensure that the system accommodates to the change of its dimensions and redistributes the density uniformly. This has been achieved through the 425 ns long equilibration simulations, see Fig. 5.

The obtained bulk moduli could now be compared to experimental values. For polystyrene, one expects to have a volumetric compressibility modulus of about 3 GPa Mark1999 . This value is, however, notably different from the fitted values shown in Fig. 6A. Although the bulk modulus for MARTINI-1.5 and MARTINI-2 is closer to the experimental data, the computed points are expected to have rather significant error bars. As was already pointed out before, the fluctuations of energy in Figs. 5A-C make it very difficult to distinguish between the individual energy curves, leading to a conclusion that measuring the energy as the function of deformation gives qualitatively a reasonable description of polystyrene volumetric compression; however, it can not be used for quantitative analysis.

An alternative approach is to compute the hydrostatic pressure of the system, which can be defined as:

[TABLE]

Here $N$ is the total number of particles in the simulation box, $V$ is its volume, $T$ is the temperature of the system, $k_{B}$ is the Boltzmann factor, $\vec{r}_{i}$ is the position of a particle with the index $i$, and $\vec{F}_{i}^{int}$ is the internal force acting on it. The hydrostatic pressure is related to the stress tensor $\hat{\sigma}$ in the following way

[TABLE]

where $\frac{1}{3}\mathop{\mathrm{Tr}}\hat{\sigma}=\frac{1}{3}(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})$ is the volumetric stress. The linear elasticity theory suggests LL7 that the stress tensor $\hat{\sigma}$ is proportional to the strain tensor $\hat{\varepsilon}$, such that

[TABLE]

where $C_{ijkl}$ is the stiffness tensor. For an isotropic medium, it has the form

[TABLE]

where $K$ and $G$ are the bulk modulus and the shear modulus, respectively. Thus, for any deformation tensor $\hat{\varepsilon}$, the volumetric stress has the form

[TABLE]

where $\mathop{\mathrm{Tr}}\hat{\varepsilon}=\delta V/V$ is the volumetric strain.

Figures 5D-F shows the time evolution of the hydrostatic (internal) pressure of the studied system at different volumetric compression regimes. The plots illustrate that in all simulations, for the three different force fields employed, the internal pressure saturates quickly with fairly small fluctuations. It is especially remarkable that in the case of the MARTINI-2 force field, it is still easily possible to distinguish between the individual curves of the internal pressure vs. simulation time dependence. Through computing the average values of the internal pressure it is ultimately possible to establish the stress-strain dependency for bulk polystyrene in the case of its volumetric compression. Here one defines:

[TABLE]

where $V_{0}$ and $P_{0}$ are the reference values of the simulation box volume and the associated internal pressure.

Figure 6B shows the dependence of the stress-on-strain for the case of a volumetric compression of the studied polystyrene system, calculated from the internal pressures shown in Fig. 5D-F. In the case of small deformations, it is expected that stress will depend linearly on strain Slaughter2002 , and the tangent of this dependency will yield the bulk modulus. Figure 6B illustrates that the obtained values for the bulk moduli for MARTINI-std, MARTINI-1.5, and MARTINI-2 force field are different from those shown in the inset to Fig. 6A, and appear to be much closer to those expected from experiment Mark1999 . In particular, the value of $K=3.7$ GPa gives a good agreement, and suggests the at the MARTINI-2 modification of the standard MARTINI force field may lead to some improvements in the computation of the bulk modulus.

III.3 Uniaxial compression

The second compression regime to consider is uniaxial compression. Equation (11) holds for any deformation and can also be applied to uniaxial deformation of the polystyrene system, yielding the bulk modulus as well. Such a uniaxial compression can also yield other elastic constants, for example, the shear modulus. In the case of the uniaxial deformation, the volumetric strain is equal to axial strain, namely

[TABLE]

Here $L_{x}$ and $L_{0}$ are the deformed and the reference values of the polystyrene box along the compression x-axis.

The results of volumetric stress-strain dependence obtained in the uniaxial compression simulation are shown in Fig. 7A. In this case, the bulk moduli are 8.2 GPa, 5.8 GPa, and 4.1 GPa for MARTINI-std, MARTINI-1.5, and MARTINI-2 force fields, respectively. These values are 7–17% higher than those obtained from volumetric compression. This difference can be caused by the local heterogeneity and anisotropy of the simulated polystyrene. It is however remarkable that in the case of a uniaxial compression the value of the bulk modulus steadily decreases for the MARTINI-std$\rightarrow$MARTINI-1.5$\rightarrow$MARTINI-2 force fields, as the fluctuations of the pressure components are smaller, providing more distinct values.

In order to obtain the shear modulus $G$, it is important to consider the pressure as a tensor with different components. The pressure tensor is related to the stress tensor as follows:

[TABLE]

where $\alpha=(x,y,z)$ and $\beta=(x,y,z)$. In the case of the uniaxial compression along the $x$ direction of the simulation box, there is only one nonzero component of the deformation tensor, $\varepsilon_{xx}$. The stress tensor, in this case, has a diagonal form with the diagonal elements given as

[TABLE]

Combining the diagonal components of the stress tensor one obtains:

[TABLE]

This equation gives an expression for the shear stress, that can be computed from the components in its left-hand side. The results for the shear stress are given in Fig. 7B. The following values for the shear stress are obtained after fitting the simulation data with the linear approximation: 1.5 GPa, 0.42 GPa, 0.007 GPa for MARTINI-std, MARTINI-1.5 and MARTINI-2 force fields, respectively. It is remarkable that the shear modulus for the system simulated with the MARTINI-2 force field becomes close to zero. This happens because the melting temperature of the polystyrene system becomes below 300 K. In this case, the interaction between individual parts of the polymeric chains becomes too weak, and the polystyrene model behaves like a liquid with zero shear modulus.

Additionally, to the bulk modulus and shear modulus, an important elastic parameter of any material is the Poisson ratio. It can be determined using the known values for the bulk modulus $K$ and the shear modulus $G$ as:

[TABLE]

The simulation results discussed above allow calculating the Poisson ratio for MARTINI-std, MARTINI-1.5 and MARTINI-2 force fields being equal to 0.4, 0.46, 0.499, respectively. All the values are a bit higher than the typical value for the polystyrene, being 0.34 Harper2002 ; Zoller1995 . It is also remarkable that the Poisson ratio value obtained for MARTINI-2 is close to 0.5, which is the expected value for the liquid state. This result suggests that once the interactions between polymeric chains are significantly reduced (as in the case of the MARTINI-2 force field) the melting temperature decreases.

IV Conclusion

The present investigation documents an approach on how to obtain elastic properties of polymer-based systems from molecular dynamics simulations. Polystyrene is considered here as a specific case study, which is investigated through the coarse-grained MARTINI force field. The reported results indicate that the conventional MARTINI force field overestimates elastic moduli (both bulk and shear moduli) as compared to experiment. We propose here that the elastic constants could be efficiently tuned through the modifications of the MARTINI force field by reducing the inter-particle van der Waals forces. Two such modifications were studied here, and it was demonstrated that significantly different values of the bulk and shear moduli for the system could be obtained. The simulations have, however, also revealed that the reduction of the van der Waals interaction leads to decreasing of the polystyrene melting temperature: once the van der Waals interactions between the coarse-grained are reduced by a factor of 2 or more, the shear modulus drops to 0, which indicates the melting of the polymeric material for the simulation temperature is below 300 K.

Based on the obtained results we, therefore, conclude that coarse-grained simulations permit performing significantly long molecular dynamics simulations of polymeric systems (here we have considered nearly 1 $\mu$s long simulations), which take slow relaxation processes into account. On the other hand, the physical correctness of the coarse-grained system turns out to be very parameter dependent, and could likely be only considered qualitatively. Indeed, we have demonstrated that both the bulk and shear moduli of polystyrene could be obtained, that would be comparable to the values known from experiment, however, it is not entirely clear how the derived force fields would apply to the modified polystyrene systems, e.g., those doped with nanoparticles. Here a followup investigation is called for, which can root upon the results presented in the present manuscript, but deserves a separate validation and justification.

Acknowledgements

Financial support by the Russian Science Foundation(RSF), grant no. 17-72-20201, is greatly acknowledged. Computational resources for the simulations were provided by the DeiC National HPC Center, University of Southern Denmark.