Exploring the role of MgO nanoparticles in the mechanical properties of the MgO/PE composite: a tight-binding and molecular dynamics simulation study

Abstract

Metal oxide nanoparticles have been widely used as reinforcing agents in polyethylene (PE)-based composites but the reinforcing mechanism is not yet fully understood. Here we report a study of the interfacial and mechanical properties of a composite composed of magnesium oxide (MgO) particles and an amorphous PE (a-PE) matrix using ab initio density functional tight binding (DFTB) and classical molecular dynamics (CMD) simulations. Through the DFTB calculations of the MgO/a-PE composite models, we reveal that the oxygen-rich termination of the MgO surface not only exhibits a much stronger attractive interaction between the MgO particles and the PE chains than the Mg-rich ones, due to the formation of interfacial O–H covalent bonds, but also provides a more favourable condition for the cross-linking between the a-PE chains. Furthermore, we demonstrate via DFTB and CMD simulations that the elastic moduli and yield stress of the MgO/a-PE composite models are obviously enhanced compared to those with a-PE, verifying the role of MgO particles as a reinforcing agent in the composite. Highlighting the importance of controlling the size and distribution of MgO nanoparticles, we believe that the present work contributes to getting atomistic insights into interfacial and mechanical properties of MgO/a-PE composites and providing a guide for developing advanced PE-based composites.

---

1 Introduction

Polyethylene (PE) has been widely used for producing various commercial goods, including high voltage insulating materials,^1–4 agricultural films,⁵ hot water pipes,⁶ implantable medical devices^7,8 and food packaging.⁹ In particular, advanced composites using PE as a matrix have attracted significant interest in several fields of industry.^10–13 This is associated with its excellent properties such as high electrical insulation, chemical corrosion resistance, good wear resistance and economical processability.^14–17 To meet the increasing demand in industry, however, advanced PE-based composites with more improved mechanical, thermal and insulating properties are needed. In this context, numerous reinforcing agents have been suggested to be incorporated into the PE matrix, including carbon nanotubes,^18–20 exfoliated graphite platelets¹⁰ and metal oxide nanoparticles (NPs) such as MgO,²¹ ZnO,^22,23 Al₂O₃ (ref. 24–26) and SiO₂ (ref. 27–29) NPs.

Recently, magnesium oxide (MgO) NPs have been given attention as a good nanofiller and reinforcing agent for PE^30–34 and other polymers^35–38 due to their high electrical insulation, high temperature stability, non-toxicity and high mechanical strength.^39,40 Moreover, there is a rich abundance of MgO-containing mineral resources in the Earth's crust and thus its cost is relatively low. El-Khatib et al.³⁰ revealed that the composites prepared with high density PE (HDPE) and MgO NPs of 5 weight (wt.)% showed optimal mechanical properties with a remarkable enhancement in tensile strength compared to pure PE. Lin et al.³¹ reported that the optimal value of MgO concentration in MgO/LDPE (low density) composites was about 5% for the best mechanical properties, revealing also that an addition of a small amount of 0.25 wt% MgO into the PE matrix leads to a significant enhancement of the static yield stress in the MgO/LDPE composites. Other experimental works also reported MgO/LDPE nanocomposites prepared with MgO weight percentages of 3–6%.^32–34 However, the underlying mechanism of mechanical strength enhancement in the PE composites by introducing MgO nanoparticles remains obscure.

Atomistic simulations play a key role in gaining an insight into material properties at the atomic scale. In relation to this, numerous works with classical molecular dynamics (CMD) simulations have been carried out to predict the physical and chemical properties of PE-based materials.^41–50 Ha et al.⁴¹ performed CMD simulations based on ReaxFF (force field) as an interatomic potential to gain an understanding of the hydrothermal gasification process of PE. Xu et al.⁴² also performed ReaxFF-based CMD simulations to investigate the pyrolysis of PE waste under both non-isothermal and isothermal conditions. Zhao et al.⁴³ investigated the dependence of the thermomechanical properties of bulk PE on its chain length using coarse-grained MD simulations. The mechanism of amorphous PE deformation was provided through CMD simulations.^44,48 The interfacial behavior of carbon nanotube fiber/PE composites was modeled and simulated by using a MD approach.⁴⁹ To the best of our knowledge, however, theoretical works for MgO/PE composites have not yet been reported, and there remains a lack of atomistic insight into the enhancement of the mechanical properties of MgO/PE composites and the role of MgO reinforcement.

In this work, we investigate the mechanical properties of composites composed of amorphous PE (a-PE) as a matrix and MgO NPs with different radii as a reinforcement using density functional tight binding (DFTB) and CMD simulations. Through DFTB simulations, we estimate the interfacial properties such as interfacial binding energy, interfacial bonding characteristics and charge transfer at the MgO/a-PE interface, together with the mechanical properties of the MgO/a-PE composite including elastic constants and moduli, as we have already applied to the ZnO/a-PTFE (polytetrafluoroethylene) composite.⁵¹ Using the larger simulation boxes for the MgO/a-PE composite, we evaluated the stress–strain curves and yield stresses through CMD simulations.

2 Methods

2.1 Atomistic modeling

The structural and mechanical properties of the MgO/a-PE composites were investigated by a combinatorial method of DFTB and CMD. For the DFTB simulations, we prepared simulation boxes containing thousands of atoms and calculated the interfacial properties at the quantum level with a reasonable accuracy and less computational cost. The larger simulation boxes, containing tens of thousands of atoms, were constructed for investigating the mechanical properties of MgO/a-PE composite at the atomic scale with CMD simulations.

To study the mechanical properties of the MgO/a-PE composites, we constructed the cubic simulation boxes, in which a MgO nanoparticle with a certain radius was placed at the centre and the surrounding region was filled with amorphous PE chains. Here, the MgO particles were made from the unit cell of the MgO crystal in the cubic phase with a space group Fm3̄m (Fig. 1(a)), using the building tool nanocluster in Materials Studio (MS) 2023 with setting the radius as 7, 8 and 9 Å, respectively. The shape of the MgO particle was selected exclusively as a sphere, although other shapes including simple box, cylinder, cone, frustum, tetrahedron and pyramid were also available. We note that the functionalities of nanoparticles are solely dependent on the size and shape.⁵² Meanwhile, PE was known to be crystallized in the orthogonal phase (c-PE) with a space group Pnma (Fig. 1(b)). The the a-PE chain was made from the ethylene monomer, we used the building tool of homopolymer in MS 2023 with a chain length parameter of 20.

Fig. 1 Ball-and-stick view of (a) the bulk MgO unit cell in the cubic phase with a space group of Fm3̄m, (b) the c-PE unit cell in the orthogonal phase with a space group of Pnma, and (c) the MgO/a-PE cubic simulation box composed of the MgO particle and surrounding a-PE chains.

Then, the simulation boxes containing one MgO particle and multiple a-PE chains were constructed by using the Amorphous Cell module in MS 2023. For the DFTB simulations, a cube box with a 25 Å edge was created, the MgO particle with a radius of 7 or 8 Å was placed at the centre of box, and then the a-PE chains were packed around the particle by setting the density parameter to be the well known PE density of 0.8 g cm⁻³.^18,42 Fig. 2 shows the constructed simulation boxes with isolated MgO particles with radii of 7 and 8 Å. For the CMD simulations, the cell length of the cube box was selected to be larger at 60 Å, and the radius of the MgO particle was chosen as 7, 8 and 9 Å. Hereafter, the composite models were named “Model r-a”, where r = 7, 8, 9 and a = 25, 60. We note that Models r-60 with the MgO particle radius of r = 7, 8 and 9 Å correspond to the reasonable MgO weight percentages of 3.17, 4.75 and 6.33 wt%, respectively, in line with experiments.^30,32–34 During the packing process, the geometries of the packed structures were optimized until the energy and force converged to 10⁻⁴ kcal mol⁻¹ and 0.005 kcal (mol⁻¹ Å⁻¹) respectively, using the smart algorithm and the condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) III forcefield.^53,54

Fig. 2 Initial configurations of the MgO/a-PE composite models, (a) Model 7-25 and (b) Model 8-25, where the MgO particle radii are 7 and 8 Å respectively while the cell length of the cubic simulation box is 25 Å. The enlarged MgO particles are shown in the bottom panel.

2.2 Computational details

For the DFTB calculations, we used the DFTB+ package (version 21.2)⁵⁵ with the 3ob-3-1 set for the electronic Slater–Koster parameters.⁵⁶ The ELPA (Eigenvalue SoLvers for Petaflop-Applications) solver provided from the ELSI (ELectronic Structure Infrastructure) library (version 2.6.1)⁵⁷ was used for the electronic structure calculation. The structural optimizations were carried out using the self-consistent charge density functional tight-binding (SCC-DFTB) method as implemented in the DFTB+ package. The special k-point meshes for Brillouin zone (BZ) integration were set to 5 × 5 × 5 and 2 × 3 × 6 for the MgO and c-PE bulk unit cells, respectively, whereas the only Γ point was used for all the calculations of the MgO/a-PE composite models. The van der Waals (vdW) dispersion between PE molecules was considered by using the Tkatchenko–Scheffler method⁵⁸ adapted for DFTB. In the structural optimizations, the threshold for the Hellman–Feynman force acting on each atom was chosen as 0.1 mHa bohr⁻¹. The mechanical properties such as elastic constants and moduli were calculated by applying the stress–strain as implemented in the ElaStic code⁵⁹ in connection with DFTB+.

We performed DFTB-MD simulations for the smaller size MgO/a-PE models shown in Fig. 2, i.e., Model 7-25 and Model 8-25. Firstly, NPT simulations for the initial configuration models were performed for 3 ps with a time step of 1 fs under the condition of 0 atm external pressure and room temperature (298 K). During the NPT simulations, we adopted the Nose–Hoover thermostat⁶⁰ and Berendsen barostat⁶¹ to control temperature and pressure, respectively. Then, NVT simulations were subsequently performed at the same temperature using the Nose–Hoover thermostat for 3 ps with a time step of 1 fs. Finally, NVE equilibrations were carried out for 3 ps with a time step of 1 fs. To speed up the MD simulations, we applied the XLBOMD (eXtended Lagrangian Born–Oppeheimer Molecular Dynamics) method,⁶² as implemented in the DFTB+ package. The mechanical properties of the MgO/a-PE models, including elastic moduli such as Young's (E), bulk (B) and shear (G) moduli, were calculated using the ElaStic code⁵⁹ in connection with the DFTB+ code. An annealing treatment was performed to account for the influence of MgO clusters on the temperature variation in the MgO/a-PE models. The annealing simulation was divided into 3 different steps; (i) the temperature was increased from 298 K to 503 K for 500 steps, (ii) the temperature was maintained at 503 K, being much higher than the melting point of PE,³⁰ for 1000 steps, and (iii) the temperature was decreased back to 298 K for 500 steps, with a time step of 1 fs.

To estimate the binding strength at the interface, we calculated the binding energy between the MgO particle and a-PE per surface area using the following equation,

where E_tot is the total energy of the MgO/a-PE model, and E_MgO and E_PE are the total energies of the MgO particle and a-PE, respectively. Here, A is the surface area of the MgO particle calculated by 4πr², which is assumed to be a radius used when constructing the MgO particle (r = 7 and 8 Å),

For the CMD simulations of the larger systems, we used the Forcite module in MS 2023 (see Table S2, ESI†). The COMPASS III (version 1.2) forcefield^53,54 was used as an interatomic potential, since it takes into account all the terms of bond (or valence), valence cross and non-bonded interactions. The forcefield parameters were determined using the ab initio density functional theory (DFT) calculations and optimized to give a good agreement with experimental data derived from large-scale databases. Firstly, the structural optimizations were performed using the smart algorithm with the fine convergence tolerance. Then, the NPT simulations were performed at 300 K and 0 GPa using the Berendsen thermostat and barostat⁶¹ with a relaxation time of 0.1 ps to equilibrate the systems (see Fig. S4–S6, ESI†). The Verlet velocity algorithm⁶³ was used to integrate the equations of motion for 1 ns with a time step of 1 fs. To obtain the stress–strain curves and thus the yield stress, the StressStrain.pl script was used as implemented in the MS software. The number of zero stress equilibration cycles and time steps were set to 5 and 100 000, respectively. The XX component of stress tensor was only considered, being gradually increased from 0.03 to 0.4 GPa with an interval of 0.03 GPa (total of 13 points). The number of balancing calculation steps corresponding to each stress value was set to 500 000.

3 Results and discussion

As a preliminary step, the structural optimizations for bulk MgO and c-PE unit cells, shown in Fig. 1(a) and (b), were carried out using the DFTB+ code (Table S1, ESI†). For the case of the bulk MgO unit cell, the lattice constant optimized with DFTB+ was determined to be a = 4.35 Å in reasonable agreement with the experimental value of 4.21 Å.^64,65 The relative error in the lattice constant is 3.33%. For the c-PE unit cell, the lattice constants were optimized to be a = 7.31, b = 4.74 and c = 2.55 Å, which also agreed well with the experimental values of a = 7.38, b = 4.94 and c = 2.54 Å,⁶⁶ with the relative errors of −0.95%, −4.05% and 3.94% respectively. These results indicate that the computational settings in the DFTB calculations are reasonable for simulations of MgO and PE compounds with the maximum relative error below 4.05%.

Then, the mechanical properties of bulk MgO and c-PE were calculated using the stress–strain method as implemented in the ElaStic code⁵⁹ in connection with the DFTB+ code. Table 1 lists the calculated elastic constants (C_ij) and elastic moduli such as bulk (B), shear (G) and Young's modulus (E) of bulk MgO and c-PE. For the case of the MgO crystal, our calculated values of elastic constants, C₁₁ = 305.9, C₁₂ = 123.8 and C₄₄ = 110.2 GPa, are found to be in good agreement with the experiment data of C₁₁ = 297, C₁₂ = 99.6 and C₄₄ = 151.5 GPa.⁶⁷ The elastic moduli were calculated to be B = 184.5, G = 102.5 and E = 259.5 GPa, in reasonable agreement with the experimental values of B = 165, G = 127.9 and E = 305.0 GPa.⁶⁷ For the c-PE bulk, the elastic constants were calculated to be C₁₁ = 7.7, C₂₂ = 11.4 and C₃₃ = 336.7 GPa, in good agreement with the experimental values of C₁₁ = 8.0, C₂₂ = 9.9 and C₃₃ = 315.9 GPa.⁶⁸ Such agreement with the available experimental data again indicates that our computational settings for DFTB calculations are reasonable for MgO and PE. When comparing between c-PE and a-PE, the former was found to have larger elastic constants and moduli than the latter in accordance with the general view that amorphous material is mechanically weaker than crystalline material. Meanwhile, the MgO crystal shows the larger values of elastic constants and moduli, except in the C₃₃ constant, than the crystalline PE as well as the amorphous PE, indicating that the bulk MgO can act as a reinforcement when making composites with PE.

Table 1 Elastic constants (C_ij) and moduli such as bulk (B), shear (G) and Young's modulus (E) of bulk MgO, crystalline PE (c-PE) and amorphous PE (a-PE), calculated using the ElaStic code in connection with DFTB+, in comparison with the available experimental data (unit: GPa)

	C11	C12	C44	B	G	E
a Ref. 67.b Ref. 65.c Ref. 68.d Ref. 18.
MgO
DFTB+	305.9	123.8	110.2	184.5	102.5	259.5
Exp.^a	297.0	99.6	151.5	165.0	127.9	305.0
Exp.^b	297.08	95.36	156.13

---

	C11	C22	C33	B	G	E
c-PE
DFTB+	7.7	11.4	336.7	23.2	13.2	33.1
Exp.^c	8.0	9.9	315.9
a-PE
DFTB+	6.0	8.8	7.2	2.9	1.7	4.3
Cal.^d						4.26

---

Upon proving the confidence of our calculations for the bulk, we performed the structural optimizations of MgO/a-PE composite models, i.e., Model 7-25 and Model 8-25 (see Fig. S1 for optimized configurations, ESI†). Then, the binding energies between the MgO cluster and a-PE chains in the optimized models were calculated using eqn (1). It was revealed that the interaction between the MgO particle and a-PE chains in the composite models was attractive since the binding energies were calculated to be negative as −51.25 and −13.35 meV Å⁻² for Model 7-25 and Model 8-25, respectively. The binding energy in Model 7-25 was found to be much larger in magnitude than that in Model 8-25, i.e., nearly four times larger, indicating that the MgO–PE interaction in the model with \ smaller radius of MgO particle is much stronger than that in the model with larger radius of MgO particle.

The reason for stronger interaction in the smaller radius MgO model is that some hydrogen atoms are dissociated from the a-PE chain and adsorbed on the surface of the MgO particle, as shown in Fig. 3. On the contrary, no hydrogen atoms bound to the MgO particle were observed in the Model 8-25. In Fig. 3(a), the transparent blue-colored circle shows the interface part of the MgO and a-PE where hydrogen atoms form O–H chemical bonds on the MgO side and as the result some C–C bonds of the PE chain change from single bonds (σ bonds) to double bonds (σ and π bonds). In the case of Model 7-25, one can observe that four hydrogen atoms are released from the carbon backbone of the PE chain, resulting in formations of two double bonds in the PE chain and four O–H covalent bonds on the MgO particle. Why did the dissociation of H atoms and formation of new bonds happen in the MgO-7 particle but not in the MgO-8 particle? This is related to the different termination characteristics of MgO particles. In the case of the MgO-7 particle, the number of O atoms is larger than the number of Mg atoms, while vice versa in the case of MgO-8 particle (see Table S2, ESI†). Therefore, the MgO-7 particle exhibits O-rich termination on its surface, whereas the MgO-8 particle shows Mg-rich termination on its surface. Based on such analysis, it can be also said that the O-rich termination of MgO particles induces a stronger binding between the MgO particle and a-PE chains than a Mg-rich termination. We note that through the same DFTB calculations the attractive interaction between the oxide and polymer was also revealed due to the negative binding energy for the ZnO/a-PTFE interfaces⁵¹ but only the interfacial covalent bond was found for them in contrast to the complicated cross-linking reactions for the MgO/a-PE interface.⁶⁹

Fig. 3 (a) Structure of Model 7-25 optimized with DFTB+. The transparent blue-colored circle shows the interface part between MgO and a-PE, where four H atoms are dissociated from the PE chain and adsorbed on the MgO surface, resulting in a change into two double bonds in the PE chain and formation of four O–H covalent bonds on the MgO particle. Green, red, dark gray and white colors represent Mg, O, C and H atoms. (b) Schematic diagram summarizing the interfacial bonding mechanism.

Upon the formation of the interface between MgO and a-PE, some electron redistribution or electron exchange occurs for bond dissociation and creation. In order to estimate the electron redistribution, we calculated the electron density difference using the formula, Δρ(r) = ρ_tot(r) − [ρ_MgO(r) + ρ_PE(r)], with the electron densities of MgO/a-PE model, MgO particle and a-PE part while fixing the atomic coordinations. Fig. 4 depicts the isosurface view of the electron density difference at the value of 0.0015|e| Å⁻³ for Model 7-25 and Model 8-25. For both of the MgO-7 and MgO-8 models, the electronic charge accumulation was found around the oxygen atoms on the surface of MgO particle, whereas the charge depletion was observed around the nearest hydrogen atoms of the a-PE chains near the MgO particle. This indicates that some electrons are transferred from the a-PE chains to the MgO particle. When comparing between the two models, the MgO-7 model exhibits distinctly higher degree of electronic charge rearrangement than the MgO-8 model. This is probably due to the complex chemical process occurring at the interface between the MgO particle and PE chains in the case of the MgO-7 model, such as chemical bond dissociation at the PE side and strong covalent bond creation at the MgO side, in good agreement with the analysis of binding energy mentioned above.

Fig. 4 Isosurface view of the electron density difference upon the MgO/a-PE interface formation in (a) Model 7-25 and (b) Model 8-25. Yellow (cyan) color represents the charge accumulation (depletion) at the value of 0.0015|e| Å⁻³.

The mechanical properties of the MgO/PE composite models, such as bulk, shear and Young's moduli, were also calculated using the ElaStic code in connection with DFTB+. Table 2 summarizes the calculated elastic moduli for the MgO/a-PE composite models. As can be expected from the elastic moduli of bulk MgO and PE in Table 1, the elastic moduli of the MgO/a-PE composite models were found to be larger than those of a-PE but lower than those of bulk MgO. Moreover, when increasing the MgO percentage, the elastic moduli of the composites were increased, such as the Young's moduli for the MgO-7 and MgO-8 models which were calculated to be 13.9 and 23.7 GPa, respectively. These indicate again that MgO can act as a reinforcing agent in the MgO/a-PE composite.

Table 2 Calculated elastic moduli, such as bulk (B), shear (G) and Young's modulus (E), for MgO/a-PE composite models, Model 7-25 and Model 8-25. The MgO percentage in the models is also presented

Model	MgO portion (wt%)	Elastic moduli (GPa)
		B	G	E
Model 7-25	31.66	15.0	4.7	13.9
Model 8-25	43.18	21.1	5.5	23.7

---

To take into account the effect of temperature on the binding characteristics in the MgO/a-PE composites, we performed DFTB-MD simulations as described in subsection 2.2 for the optimized models, checking the number of hydrogen atoms (N_H) forming O–H covalent bonds on the MgO side (see Fig. S1, ESI†). We first performed DFTB-MD NPT simulations of the two models at room temperature, T = 298 K. In the case of Model 7-25, N_H was found to be increased from 4 at absolute zero temperature to 8 at room temperature (Fig. S2(b), ESI†), meaning that the number of unsaturated bonds in the PE chains in the vicinity of the interface with the MgO particle was increased and thus the binding strength was enhanced. On the contrary, for the case of Model 8-25, N_H was found to be still zero at room temperature (see Fig. S2(e), ESI†), indicating that the binding characteristics are not affected by increasing temperature. We then performed annealing dynamics simulations by setting the initial temperature as 300 K and the maximum temperature (mid-cycle temperature) as 500 K, which is higher than the melting temperature of PE (407 K), to ensure that a-PE chains are completely melted. After the annealing treatment, N_H in the MgO-7 model was increased by 12, meaning there was an increase of the unsaturated (double) bonds in the PE chains (see Fig. S2(c), ESI†) and thus an enhancement of binding strength. For the case of the MgO-8 model, however, N_H was still zero, meaning that no hydrogen atoms were dissociated from the PE chain and thus the O–H covalent bonds were not created on the surface of the MgO particle (see Fig. S2(f), ESI†). From these results, it can be said that the MgO nanoparticles with O–rich surfaces provide favourable conditions for a formation of cross-links between MgO and PE chains even at finite temperature.

In order to directly verify the enhancement of mechanical strength by making composites with MgO particles and an a-PE matrix, we performed CMD simulations of larger simulations boxes, i.e. Model r-60, where the radii of the MgO particles were r = 7, 8, 9 Å and the cell length was 60 Å, as shown in Fig. 5. As a preliminary step, the NPT simulations were performed for 1 ns with a time step of 1 fs, using the COMPASS III force field as an interatomic potential. It was confirmed that the simulation time of 1 ns was long enough to equilibrate the temperature, pressure and density of the simulation boxes (Fig. S3 and S4, ESI†). The densities of the composite models were found to gradually increase with the increasing the radius of the MgO particle from 0.850 to 0.866 and 0.874 g cm⁻³ for r = 7, 8, 9 Å in accordance with expectations.

Fig. 5 Configurations of (a) Model 7-60, (b) Model 8-60 and (c) Model 9-60 with the cell length of 60 Å.

The structural characteristics of the composite models were analyzed by calculating the radial distribution function (RDF). Fig. 6 shows the total RDF for the MgO particles and the RDF for O–H intermolecular interactions in Model r-60 (r = 7, 8, 9 Å). As shown in Fig. 6(a), the total RDFs for MgO particles with radii of 7, 8 and 9 Å exhibited overall similar features with almost similar positions of the main peaks. However, the intensity of the RDF was found to gradually decrease with increasing the radius of the MgO particle. Moreover, it was found that the MgO-7 model showed a distinct peak at a radius of 2.65 Å for the intermolecular O–H interaction as shown in Fig. 6(b). In accordance with the DFTB+ calculation, the outermost O atoms of the MgO particle were found to move towards the outside of the MgO particle in the Model 7-60 while they moved towards the inside of the MgO particle in the MgO-8 and MgO-9 models after NPT equilibration (see Fig. S5, ESI†).

Fig. 6 Radial distribution function for (a) the MgO particle and (b) the O–H interaction in Model r-60 with r = 7, 8, 9 Å.

After equilibrations of the simulation boxes, we performed the StressStrain.pl script to obtain the stress–strain curves and the yield stress as a check of mechanical strength. Fig. 7 shows the stress against the applied strain, averaged over the 5 production runs for the PE-60 and the three Model r-60 (r = 7, 8, 9 Å). The calculation data were fitted to a logarithmic function of y = a ln(x + b) + c with the fitting parameters of a, b and c providing estimates for the yield stress. For all the simulation models, it was found that by increasing the strain value the stress value increased rapidly at the beginning and arrived at the saturation value with strain values of over 0.3%. One can find that for the Model r-60 the gradients of the tangent lines at the origin and the saturated stress values are clearly larger than those for the pure PE-60 model. Meanwhile, the derivative of the stress–strain curve gave the strain-dependent Young's modulus (Fig. S6, ESI†), which were fitted to a power function of y = ax^b with the fitting parameters of a and b (Table S3, ESI†). With these fitting functions and parameters, we could determine the Young modulus (E), yield stress (Y) and critical strain (ε_c) with increasing the portion of MgO reinforcement in the MgO/a-PE composites, as listed in Table 3. It is worth noting that MgO crystals are known to be hygroscopic in nature and the hygroscopic feature of MgO should cause a reduction of mechanical strength of MgO/PE composites, because the ingress of water in composites increases separation between the molecular chains, leading to expansional strain and thus degradation.^70,71

Fig. 7 Stress–strain diagrams of (a) the PE-60 model, (b) Model 7-60, (c) Model 8-60 and (d) Model 9-60 for MgO/PE composites, calculated at the room temperature of 298 K. Red-colored lines represent the fitting to a logarithmic function of y = a ln(x + b) + c with the fitting parameters of a, b and c.

Table 3 Calculated Young's modulus (E), yield stress (Y) and critical strain (ε_c) for the pure PE-60 model and composite models, Model r-60 with r = 7, 8, 9 Å, together with the MgO weight percentage (wt%)

Model	MgO portion (wt%)	E (GPa)	Y (Gpa)	εc (%)
PE-60	0.00	4.563	0.178	0.45
Model 7-60	3.17	21.083	0.203	0.30
Model 8-60	4.75	19.031	0.262	0.34
Model 9-60	6.33	17.946	0.193	0.27

---

For the pure a-PE model, the yield stress was estimated to be 0.178 GPa in reasonable agreement with the previously reported experimental value of 0.145 GPa (ref. 44) at the critical strain value of 0.45%, and the Young's modulus was determined to be 4.563 GPa in good agreement with the previously calculated value of 4.26 GPa (ref. 18) and our present DFTB+ calculation value of 4.3 GPa. When including the MgO particle into the a-PE matrix, both the yield stress and Young's modulus were obviously found to be increased, while the critical strain values were decreased, indicating the enhancement of the mechanical strength by making the MgO/a-PE composite. It was found that with increasing the MgO portion in the composite (i.e. the radius of the MgO particle) the modulus decreased from 21.083 GPa for 3.17 wt% to 19.031 GPa for 4.75 wt% and to 17.946 GPa for 6.33 wt%. This tendency is in accordance with the binding characteristics between the MgO particle and the a-PE matrix discussed above. On the other hand, the yield stress was found to have a maximum value of 0.262 GPa for Model 8-60, with smaller values of 0.203 and 0.193 GPa for Model 7-60 and Model 9-60, respectively. We note that such a tendency was also experimentally found in other composites, for example, the multi-wall carbon nanotube (MWCNT)-enhanced poly(methyl methacrylate) (PMMA) composites.⁷² From our results for mechanical strength, we can conclude that a critical MgO content exists for effectively enhancing the mechanical strength of the composite, at least from the yield stress point of view. For this reason, we suggest that the reduction of mechanical strength with increasing the MgO content is mainly associated with the interfacial binding strength between the MgO particle and PE chains. When increasing the MgO content, the interfacial binding energy was decreased in magnitude, indicating the decrease of interfacial binding strength related to the termination of the MgO particle. In fact, the MgO-7 particle shows an O-rich termination, causing a cross-linking reaction at the interface, whereas the MgO-8 and MgO-9 particles have a Mg-rich termination, causing weaker interfacial binding and leading to poor MgO-PE wetting. For to high a MgO content, the MgO particles induce stress concentrations during deformation, probably causing the strength reduction. This indicates that the control of the size and distribution of the MgO reinforcing agent in the composite plays an important role in enhancing the mechanical strength.

4 Conclusions

In this work, we have performed ab initio DFTB and classical molecular dynamics simulations of MgO/a-PE composite models to reveal the role of MgO nanoparticles in the interfacial and mechanical properties of the composite. For the simulations, we constructed the MgO/a-PE composite models by introducing a MgO nanoparticle with different radii of r = 7, 8 and 9 Å into the cubic box with different cell edges of a = 25 and 60 Å containing the a-PE chains with an initial density of 0.8 g cm⁻³. From the binding energy calculations and electronic density difference analysis using the optimized geometries with DFTB, we revealed that the O-rich termination of the MgO particle surface induces dissociation of H atoms from the adjacent a-PE chains and the formation of new O–H covalent bonds, resulting in much stronger attractive binding between MgO and PE compared with the Mg-rich ones. Through the DFTB-MD simulations at room temperature and with annealing treatment, we demonstrated that the MgO nanoparticles with O-rich termination provided favourable conditions for the formation of the cross-linking between MgO particle and PE chains. Furthermore, we calculated the mechanical properties such as elastic moduli and yield stress, finding that the mechanical strength of MgO/a-PE composites can be significantly enhanced by introducing MgO particles into a-PE matrix. Highlighting the importance of controlling the size and distribution of MgO nanoparticles in the MgO/a-PE composite for enhancing the mechanical strength, we believe that the present work contributes to gaining atomistic insights into interfacial and mechanical properties of MgO/a-PE composites and providing a guide for developing advanced PE-based composites.

Data availability

The data supporting this article have been included as part of the ESI.†

Author contributions

Chol Ryu performed the DFTB calculations and drafted the first manuscript. Song-Mu Kim, Il-Ung Kim and Jun-Gi Ri assisted with the post-processing of calculation results and contributed to the useful discussions. Chol-Jun Yu performed CMD simulations, drafted the first manuscript and supervised the work. All authors reviewed the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work is supported as part of the fundamental research project “Design of New Energy Materials” (no. 2021-5) from the State Commission of Science and Technology, DPR Korea. Computation in this work was done on the HP Blade System C7000 (HP BL460c) that is owned and managed by Faculty of Materials Science, Kim Il Sung University.

References

1. F. Nilsson, M. Karlsson, U. W. Gedde, R. Kádár, K. Gaska, C. Müller, P.-O. Hagstrand, R. T. Olsson, M. S. Hedenqvist and T. Gkourmpis, Composites, Part B, 2021, 204, 108498–108506.
2. R. Gupta, L. Smith, J. Njuguna, A. Deighton and K. Pancholi, ACS Appl. Electron. Mater., 2020, 2, 1880–1891.
3. Y. Wang, C. Wang, Z. Zhang and K. Xiao, Nanomaterials, 2017, 7, 320–332.
4. G. Mazzanti, G. C. Montanari, F. Palmieri and J. Alison, J. Appl. Phys., 2003, 94, 5997–6004.
5. Y. Chen, S. Chai, H. Tian, Y. Chai, Y. Li, L. Chang and H. Cheng, Agric. Manage. Water, 2019, 211, 142–151.
6. J. C. Montes, D. Cadouxa, J. Creusb, S. Touzainb, E. Gaudichet-Maurinc and O. Correc, Polym. Degrad. Stab., 2012, 97, 149–157.
7. K. Ishihara, Polym. J., 2015, 47, 585–597.
8. B. L. Armstrong, A. F. Senyurt, V. Narayan, X. Wang, L. Alquier and G. Vas, J. Pharm. Biomed. Anal., 2013, 74, 162–170.
9. A. Nasiri, E. Gastaldi, N. Gontard and S. Peyron, Appl. Clay Sci., 2020, 198, 105803–105812.
10. Y. Liu, J. Zhang, X. Wang, Y. Liu, X. Hu, C. Cao, X. Qu and B. Abdel-Magid, RSC Adv., 2024, 14, 27948–27956.
11. Y.-F. Huang, J.-Z. Xu, J.-Y. Xu, Z.-C. Zhang, B. S. Hsiao, L. Xu and Z.-M. Li, J. Mater. Chem. B, 2014, 2, 971–980.
12. Y. Gao, X. Huang, D. Min, S. Li and P. Jiang, ACS Sustain. Chem. Eng., 2019, 7, 513–525.
13. P. Khodaparast and Z. Ounaies, Smart Mater. Struct., 2014, 23, 104004.
14. F. D. Stefano, M. Baur, C. D. Rosa and S. Mecking, Macromolecules, 2024, 57, 1072–1079.
15. C. Zou, Q. Wang, G. Si and C. Chen, Nat. Commun., 2023, 14, 1442.
16. M. Andersson, J. Hynynen, M. Andersson, V. Englund, P. A. Hagstrand, T. Gkourmpis and C. Müller, ACS Macro Lett., 2017, 6(2), 78.
17. M. Sugimoto, A. Shimada, H. Kudoh, K. Tamura and T. Seguchi, Radiat. Phys. Chem., 2013, 82, 69–73.
18. R. Chawla and S. Sharma, Compos. Sci. Technol., 2017, 144, 169–177.
19. Z. M. Dang, S. H. Yao and H. P. Xu, Appl. Phys. Lett., 2007, 90, 012907.
20. S. P. Bao, G. D. Liang and S. C. Tjong, Carbon, 2011, 49, 1758–1768.
21. L. Pallon, A. Hoang, A. Pourrahimi, M. Hedenqvist, F. Nilsson and S. Gubanski, et al., J. Mater. Chem. A, 2016, 4, 8590–8601.
22. J.-G. Gao, X. Li, W.-H. Yang and X.-H. Zhang, Crystals, 2019, 9, 481.
23. B. P. Chang, H. M. Akil, R. B. M. Nasir, I. M. C. C. D. Bandara and S. Rajapakse, J. Reinf. Plast. Compos., 2014, 33, 674–686.
24. J. P. Jose and S. Thomas, Phys. Chem. Chem. Phys., 2014, 16, 14730–14740.
25. S. Zhang, X. Y. Cao, Y. M. Ma, Y. C. Ke, J. K. Zhang and F. S. Wang, eXPRESS Polym. Lett., 2011, 5, 581–590.
26. V. Bertolino, G. Cavallaro, G. Lazzara, G. Lazzara, S. Milioto, F. Parisi, E. V. Rozhina and R. F. Fakhrullin, Ind. Eng. Chem. Res., 2016, 55, 7373–7380.
27. G. Cavallaro, A. A. Danilushkina, V. G. Evtugyn, G. Lazzara, S. Milioto, F. Parisi, E. V. Rozhina and R. F. Fakhrullin, Nanomaterials, 2017, 7, 199.
28. K. Y. Lau, A. S. Vaughan, G. Chen, I. L. Hosier and A. F. Holt, J. Phys. D, 2013, 46, 095303.
29. D. Panaitescu, F. Ciuprina, M. Iorga, A. Frone, C. Radovici, M. Ghiurea, S. Sever and I. Plesa, J. Appl. Polym. Sci., 2011, 122, 1921–1935.
30. A. M. El-Khatib, M. M. Gouda, M. S. Fouad, M. Abd-Elzaher and W. Ramadan, Sci. Rep., 2023, 13, 9945.
31. X. Lin, Y. Wu, L. Y. Tang, M. H. Yang, D. Y. Ren, J. W. Zha and Z. M. Dang, J. Appl. Polym. Sci., 2016, 133, 43038–43045.
32. A. Samad, K. Y. Lau, I. A. Khan, A. H. Khoja, M. M. Jaffar and M. Tahir, J. Phys. Chem. Solids, 2018, 120, 140–146.
33. F. Nilsson, M. Karlsson, L. Pallon, M. Giacinti, R. T. Olsson, D. Venturi, U. W. Gedde and M. S. Hedenqvist, Comput. Sci. Technol., 2017, 152, 11–19.
34. L. Zhang, Y. Zhou, J. Tian, Y. Sha, Y. Zhang, H. Wu and Y. Wang, J. Electrost., 2014, 72, 252–260.
35. E. Natarajan, T. Sekar, K. Markandan, S. Nesappan, A. D. A. Selvaraj and S. M. Sekar, Emergent Mater., 2024, 7, 1283–1294.
36. V. N. Hegde, T. M. Pradeep, V. V. Manju and N. C. Sandhya, Mater. Sci. Eng. B, 2024, 310, 117775.
37. S. Kumari, R. K. Mishra, S. Parveen, S. K. Avinashi, A. Hussain, S. Kumar, M. Banerjee, J. Rao, R. Kumar, R. K. Gautam and C. Gautam, Sci. Rep., 2024, 14, 2128.
38. A. M. E. Sayed, A. M. Abdelghany and A. A. Elfadl, Braz. J. Phys., 2022, 52, 150.
39. M.-A. Gatou, E. Skylla, P. Dourou, N. Pippa, M. Gazouli, N. Lagopati and E. A. Pavlatou, Crystals, 2024, 14, 215–261.
40. M. S. S. Danish, A. Bhattacharya, D. Stepanova, A. Mikhaylov, M. L. Grilli, M. Khosravy and T. Senjyu, Metals, 2020, 10, 1604.
41. D. T. Ha, H. D. Tong and T. T. Trinh, Sci. Rep., 2024, 14, 18771–18787.
42. Z. Xu, Q. Xie, C. Chen and X. Jiang, Polym. Degrad. Stab., 2023, 208, 110249.
43. J. Zhao, S. Nagao and Z. Zhang, J. Mater. Res., 2010, 25, 537–544.
44. D. Hossain, M. A. Tschopp, D. K. Ward, J. L. Bouvard, P. Wang and M. F. Horstemeyer, Polymer, 2010, 51, 6071–6083.
45. A. Henry and G. Chen, Phys. Rev. Lett., 2008, 101, 235502.
46. A. Henry and G. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 144305.
47. S. Buell, K. J. Van Vliet and G. C. Rutledge, Macromolecules, 2009, 42, 4887–4895.
48. C. Peng and F. Zeng, Comput. Mater. Sci., 2017, 137, 225–232.
49. Y.-B. Lu, Q.-S. Yang, X.-Q. He and K.-M. Liew, Comput. Mater. Sci., 2016, 114, 189–198.
50. A. Moyassari, T. Gkourmpis, M. S. Hedenqvist and U. W. Gedde, Polymer, 2019, 161, 139–150.
51. C. Ryu, J.-G. Ri, Y.-S. Kim, C.-H. Rim, C.-I. Kim and C.-J. Yu, RSC Adv., 2024, 14, 35097–35103.
52. S. S. Singh, V. Parameswaran and R. Kitey, Composites, Part B, 2019, 164, 103–115.
53. R. L. C. Akkermans, N. A. Spenley and S. H. Robertson, Mol. Simul., 2021, 47, 540–551.
54. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364.
55. B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitric and A. Dominguez, et al., J. Chem. Phys., 2020, 152, 124101.
56. M. Gaus, A. Goez and M. Elstner, J. Chem. Theory Comput., 2013, 9, 338–354.
57. V. Yu, F. Corsetti, A. García, W. P. Huhn, M. Jacquelin, W. Jia, B. Lange, L. Lin, J. Lu, W. Mi, A. Seifitokaldani, A. Vázquez-Mayagoitia, C. Yang, H. Yang and V. Blum, Comput. Phys. Commun., 2018, 222, 267–285.
58. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
59. R. Golesorkhtabar, P. Pavone, J. Spitaler, P. Puschnig and C. Draxl, Comput. Phys. Commun., 2013, 184, 1861–1873.
60. G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Mol. Phys., 1996, 87, 1117–1157.
61. H. J. C. Berendsen, J. P. M. Postma, W. F. V. Gunsteren, A. Dinola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
62. B. Aradi, A. M. N. Niklasson and T. Frauenheim, J. Chem. Theory Comput., 2015, 11, 3357–3363.
63. L. Verlet, Phys. Rev., 1967, 159, 98–103.
64. R. R. Reeber, K. Goessel and K. Wang, Eur. J. Mineral., 1995, 7, 1039–1047.
65. CRC Handbook of Chemistry and Physics, ed. D. R. Lide, CRC Press, Boca Raton, FL, 2005, http://www.hbcpnetbase.com.
66. G. T. Davis, R. K. Eby and J. P. Colson, J. Appl. Phys., 1970, 41, 4316.
67. A. J. Cinthia, G. S. Priyanga, R. Rajeswarapalanichamy and K. Iyakutti, J. Phys. Chem. Solids, 2015, 79, 23.
68. K. Tashiro, M. Kobayashi and H. Tadokoro, Macromolecules, 1978, 11, 914–918.
69. H. Ahmad and D. Rodrigue, Polym. Eng. Sci., 2022, 62, 2376–2401.
70. J. Garcia-Espinel, D. Castro-Fresno, P. P. Gayo and F. B.-M. noz, Mater. Des., 2015, 66, 46–50.
71. V. B. Mora, M. Mieloszyk and W. Ostachowicz, Mech. Syst. Signal Process., 2018, 99, 534–549.
72. F. Sun, U. Wiklund, F. Avilés and E. K. Gamstedt, Compos. Sci. Technol., 2017, 153, 95–102.

---

Footnote

† Electronic supplementary information (ESI) available: Tables for optimized lattice constants of bulk MgO and c-PE, simulation cell parameters and fitting parameters, and figures for equilibrated structures of MgO/a-PE models, variations of temperature, pressure and density, and strain-dependent Young's modulus. See DOI: https://doi.org/10.1039/d5ra02394g

---