Physically constrained random contact-based model construction
In the CVD method, the skeleton of GF is constructed by stacked graphene layers. As the composite is fabricated by immersing GF in PDMS solution, the connection between GF and PDMS mainly relies on Van der Waals’ force4.
The physical model utilized in this simulation is based on the experimental preparation of the GF/PDMS composite. Limited by the computational capability of the Vienna Ab initio Simulation Package (VASP), PDMS models are constructed by repeating units ranging from 3 to 8. The molecular structures of PDMS with different units are presented in Fig. 1a. The GF model mimics the layer-by-layer structure. Additionally, a supercell is employed to guarantee structural uniformity and the non-randomness of cross-plane heat transfer.
a Chain structure of PDMS molecules with 4, 6 and 8 repeating units; b Point, edge and face contacts between PDMS with graphene; c GF with varying PDMS doping rates in MD simulations.
To comprehensively explore all possible integration of PDMS with GF without compromising computational efficiency, a stochastic contact algorithm was implemented through Python code. The primary interaction sites between PDMS and graphene are CH–π interactions, where methyl hydrogen atoms in PDMS engage with the π-electron cloud of graphene. These interactions, supplemented by potential crosslinking between the CH3 side groups of PDMS and graphene, form the basis of the interfacial modeling strategy26,27.
Given the localized nature of these interactions, our method strategically places PDMS chains via three distinct contact modes-point, edge, and face contacts- anchored at randomly selected carbon atoms on the graphene surface, as shown in Fig. 1b. Each attachment is geometrically optimized to prevent atomic overlaps, with molecular orientations sampled from a uniform distribution. This approach effectively captures the distribution of interfacial configurations observed in experimental composites, as validated by the dynamic behavior of contact points during MD equilibration.
Following stochastic placement, energy minimization is performed to relax the system into stable adhesion structures. The resulting supercell models (60Å × 60Å × 300Å) systematically vary PDMS doping from 0% to 10%, with all PDMS molecules bound to the GF framework through non-covalent interactions, as depicted in Fig. 1c. Detailed parameters for the five representative models are provided in Table 1.
Neuroevolution potential training
Traditional empirical potentials for MD simulations exhibit limited transferability to the continuous compression behavior of GF/PDMS composites. This limitation arises from the heterogeneous organic-inorganic interfaces and multiscale structural evolution inherent to such hybrid systems. To address this gap, we developed a machine-learning potential (MLP) trained on quantum mechanical datasets derived from the DFT calculations, as shown in Fig. 2. This approach leverages the Neuroevolution Potential (NEP-4), which integrates the unique advanced Seperable Natural Evolution Strategy (SNES) for the global optimization of atomic interactions. The NEP-4 architecture achieves sub-meV/atom accuracy in energy predictions while showing higher computational efficiency comparable to classical MD potentials28,29. Hence, the NEP-4 model was adopted for potential training in our research.
The process involves: (i) obtaining configurations through AIMD or MD; (ii) screening configurations that cannot be adequately described by the potential; (iii) training the potential function through active learning to achieve a comprehensive description of nearly all structural features of the system; (iv) conducting thermodynamic and mechanical property analyses based on the optimized potential function.
The interatomic potential was trained on a multiscale dataset encompassing defective and pristine CVD-fabricated graphene, pure PDMS chains with tunable crosslinking densities, and GF/PDMS nanocomposite interfaces. Figure 3a presents the scanning electron microscope (SEM) characterization of GF/PDMS. The initial model trained on the dataset accurately reproduces structural features at the microscopic scale. The number of atoms in all configurations ranges roughly from 40 to 150. Ab Initio Molecular Dynamics (AIMD) generated initial configurations for GF/PDMS systems. The VASP software suite (DFT-based) provided reference energies, forces, and virial tensors. Using the NPT ensemble, we perturbed and simulated configurations in the training set to filter out cases inadequately described by the potential. Representative structures were selected via farthest-point sampling30,31, with their energy, force, and virial data fed into training and testing sets to iteratively improve model accuracy. The Bussi-Donadio-Parrinello thermostat and Bernetti-Bussi barostat maintained stable temperature and pressure conditions during simulations, while GPUMD (a GPU-accelerated MD tool) enabled efficient training32,33.
a SEM images of GF/PDMS composite45 (Reproduced with permission of Elsevier, 2016); b Loss functions of training and testing datasets with the number of training steps; Plots comparing NEP predictions against DFT reference values for energy (c), atomic forces (d), and virial stresses (e); f Convergence check results in atomic forces; g The comparison of Reaxff and DFT; h Comparison of computational speeds between AIMD and MLP-NEP; i Size dependence of thermal conductivity in GF/PDMS composites with a 5% doping rate at different sizes of 30 Å, 45 Å, 60 Å, 75 Å, 300 Å and 600 Å.
The fitting results of the loss function are shown in Fig. 3b. Convergence behavior ensures consistency between computed values and fitted data. Comparative analysis of DFT versus NEP predictions, shown in Fig. 3c–e), revealed sub-meV/atom energy root-mean-square errors (RMSEs) and sub-0.1 eV/Å force deviations, consistent with state-of-the-art NEP benchmarks23,33. The diagonal clustering of data points confirms the potential accuracy in modeling GF/PDMS interfacial dynamics.
To validate the robustness of the potential, we performed convergence checks during the final 10 ns of MD simulations, compared the differences in energy, force, and virial obtained by NEP and DFT under extreme deformation conditions, ensuring that the prediction errors in energy, forces and virials remained within the bounds established by the training data. Testing configurations (uniformly sampled) showed that NEP-predicted energy, force, and virial errors remained within training data bounds. As shown in Fig. 3f, the root-mean-square error (RMSE) for forces between NEP and DFT was 291.3 meV/Å, confirming the potential’s accuracy and stability23,24,25.
To further assess the accuracy of our model, we used part of the training dataset to compare the atomic forces predicted by the trained NEP-MLP with those from a reactive force field34 (Reaxff) (Fig. 3g) based on previously reported force field. The Reaxff-predicted atomic forces exhibited a significantly higher RMSE of 7170.9 meV/Å compared to NEP-MLP’s 384.2 meV/Å (Fig. 3g), clearly demonstrating the machine learning potential’s superior accuracy and reliability for capturing complex interactions in GF/PDMS systems.
With comparable accuracy and identical computational resources, NEP-MLP achieves a speedup of ~30,000,000 times over traditional AIMD (Fig. 3h), providing the computational foundation for the discoveries presented in this work. Due to the substantial computational variability observed in small systems (as discussed in Section “Strain-dependent mechanics of GF/PDMS composites”), simulations involving large-scale systems—on the order of tens of thousands of atoms—are essential. In this context, the extremely low time complexity of NEP provides a significant advantage, enabling efficient and accurate calculating of large systems.
Size-convergence validation
In the calculation of thermal conductivity based on MD, it has been commonly observed that the size effect tends to induce an underestimation of the thermal conductivity values of materials. This phenomenon indicates that the size of the material must be greater than the average free path of phonons. To mitigate finite-size effects in non-equilibrium molecular dynamics (NEMD) thermal conductivity calculations, the GF/PDMS composite with a 5% doping rate is systematically evaluated across six system sizes ranging from 30 Å to 600 Å, as shown in Fig. 3f. Sub-75Å model exhibits atomic crossovers between group layers during the relaxation process due to the intergroup displacements less than 6 Å, distorting temperature gradient calculations.
Above 75 Å, the proportion of atoms involved in crossovers decreases exponentially (R2 = 0.97) with system size, stabilizing thermal conductivity values within 2% variance beyond 300 Å, as depicted in Fig. 3i. Limited by computational resource constraints, the 300 Å model is selected in the calculations here.
Strain-dependent mechanics of GF/PDMS composites
To explore the mechanical and thermal response of GF/PDMS composites under uniaxial strain, the compression and stretching are imposed along the z-axis direction on the GF with different PDMS doping rates. The deformation features of the GF/PDMS composite are depicted in Fig. 4a–e. Fixed boundary conditions in x-y planes eliminate lateral expansion artifacts. In the process of compression, the squeezing effect between atoms accumulates continuously until the box cracks irreversibly. Terminating compression at 50% prevents cell fracture. It shows that the graphene layers experience bending and wrinkling during compression. The inter-layer space is squeezed until the graphene layers stack together. In the process of stretching, the strain increases with the stretching length. Percolative fracture networks emerge at 40% tensile strain, which leads to the disruption of continuous thermal pathways. Therefore, the stretching length is confined to 40%.
a Pristine structure; b Structure 25% compression ratio; c Structure with 50% compression ratio; d Structure 20% stretching ratio; e Structure with 40% stretching ratio; f Tensile stress variation with varying PDMS doping rates under stretching; g Maximum stretchability rate with varying PDMS doping rates; h Tensile stress variation with varying PDMS doping rates under compressions; i The relationship between PDMS doping rates and tensile modulus, tensile strength, and fracture strain.
Five models selected from Table 1 serve as base models, within 0.5 ns stretching and compressive deformations applied along the z-axis. Figure 4f–i illustrate their evolution from the initial state to fracture under varying PDMS doping rates. As doping levels increase, the maximum stretching ratio gradually improves, indicating that PDMS effectively enhances the tensile strength of pristine GF. A critical doping threshold emerges between 2.5% and 5%, where the maximum tensile strain increases by 4.5%. In the 2.5% doped model, insufficient PDMS quantity results in inadequate interlayer connections between graphene layers. Consequently, the graphene undergoes fracture at the unbonded regions where PDMS is absent, and thus the stretching curve almost mirrors that of pure GF. At 5% doping, PDMS bridges nearly all graphene layers, enabling load transfer significantly that elevates peak tensile strength to 1.7 GPa at 10% doping. Force curve fluctuations indicate graphene layer fracture: incomplete force drops to 0 GPa signal sustained interlayer connections via PDMS chains. Additionally, the frequency of stress drops also increases with doping levels. Figure 4i quantifies mechanical parameter improvements of 20% rise in tensile modulus and tensile strength in the z-direction with PDMS doping rates. While fracture strain improvements are subtle, consistent enhancements are observable. The behavior during compression follows a trend that is observed during stretching. As illustrated in Fig. 4h, at equivalent strain, higher doping reinforces compressive pressure (negative force), attributed to the rigid PDMS framework restricting the global deformation. Notably, the 5% PDMS compressive stress curve lies significantly below 0% and 2.5% counterparts, underscoring its role in forming a stable rigid structure that stabilizes thermal conduction pathways consistent with its mechanical reinforcement function.
Thermal transport modulation in GF/PDMS composites
While PDMS enhances the mechanical properties of GF, its role as a dopant in modulating phonon-mediated heat transport remains critical for dynamic thermal management applications. Here, we measured the thermal conductivity and thermal conductance of five initial composite models using NEMD simulations. The simulation box is partitioned into 13 axial groups along the heat-flow direction, with a temperature gradient imposed along the z-axis by defining a hot plane (330 K) and a cold plane (270 K). Heat transfer is induced via kinetic energy exchange between particles in these planes, such that energy removed from the hot region is exactly balanced by energy added to the cold region. A linear increase in this energy exchange over time indicates the establishment of steady-state thermal transport.
Figure 5a illustrates the temporal evolution of energy in the hot (blue line) and cold (red line) regions. During the initial 0.5 ns, non-uniform heat distribution prevented steady-state conditions. After this transient period, the thermal transfer of the system gradually stabilizes, and the heat flux between hot and cold regions approaches equilibrium. Figure 5b shows the temperature profile of each axial group from the onset of steady state (0.5 ns) to 2 ns, demonstrating uniform thermal gradient development across the system.
a Temporal evolution of energy gain (red) and loss (blue) in cold and hot regions; b Z-direction temperature profile across axial groups from the establishment of steady-state thermal transport (0.5 ns) to 2 ns; c Variation of thermal conductivity and thermal conductance of GF with different doping rates of PDMS (T = 300 K); d Equivalent thermal circuit of pure GF (without PDMS); e Comparison of heat-flow pathways with and without PDMS; f Equivalent thermal circuit of GF doped with PDMS; g Thermal conductivity variation with the varying deformation; h Thermal conductance variation with the varying deformation.
The thermal properties of GF doped with PDMS are investigated under both compression and stretching. Figure 5c shows the increasing doping rate of PDMS directly results in the decline of the total thermal conductivity of the composite. This is because the thermal conductivity of PDMS (0.14 W/(m·K)) is an order of magnitude lower than that of pure GF at the identical temperature (T = 300 K) and density (ρ = 0.49 g/cm3). Here, thermal resistance analysis is employed to reveal the underlying mechanism of thermal property variation with doping rates. The thermal circuit of pristine GF, composed of stacked graphene layers, is displayed in Fig. 5d. Each of the graphene layers acts as a parallel thermal resistance in the heat-transfer process. When PDMS molecules are introduced, regardless of the contact form, additional branch thermal resistance arises at graphene-PDMS interfaces. Figure 5e compares the heat-flow pathways with and without PDMS, while Fig. 5f depicts the equivalent thermal circuit after PDMS doping. It shows PDMS incorporation increases the total thermal resistance, which is the primary reason for the decline in thermal conductivity. However, the thermal resistance increases moderately for limited PDMS loadings, leading to an insignificant reduction in the overall thermal conductivity.
The thermal conductivity and the thermal conductance of GF under compression are measured, as shown in Fig. 5g, h. Results indicate that across all PDMS doping rates, the thermal conductivity λ decreases as the model height reduces, aligning with increased contact thermal resistance from gradual graphene layer stacking during compression. The simulated thermal conductivity values are in good agreement with experimental results, exhibiting the same order of magnitude as the literature-reported data: 0.94 W/(m·K) at 8vol.% PDMS doping and 0.56 W/(m·K) at 0.7 wt.% PDMS doping. This close correspondence validates the reliability of our simulation framework19,35,36. Within the Tight Binding model framework37, electrons are confined near atomic nuclei with interactions primarily between adjacent atoms. Morphological changes induced by PDMS incorporation and composite compression enhance van der Waals forces between graphene layers, introducing stronger hindrances to atomic thermal motion. It amplifies phonon-phonon scattering during propagation, shortening phonon lifetimes and reducing thermal conductivity.
Thermal conductance variations across models with differing PDMS doping rates exhibit an intriguing trend. Over the deformation spectrum from 50% compression to 40% stretching, the thermal conductivity scales approximately linearly with model length. Since thermal conductance depends on the ratio of thermal conductivity to model length, thermal conductance is governed primarily by the relative shifts in thermal conductivity. Below a 20% stretching ratio, the thermal conductance shows a minor initial decrease followed by an increase, a behavior rationalized by a spring model38 rooted in the effective medium theory.
During the compression from the initial state to 50%, the limited formation of new internal bonds, especially in low-density CVD models, is offset by increased thermal pathways, driving a continuous rise in overall thermal conductance consistent with computational results. When the stretching ratio is below 20%, the gradual alignment of graphene chains equalizes thermal resistances across pathways, enhancing the uniformity of heat flow and thermal efficiency, which are key factors for the concurrent increase in both thermal conductivity and thermal conductance. Once the stretching ratio exceeds 20%, although thermal conductivity continues to rise, the overall thermal conductance declines across all doping rates. The decline stems from the fracture of some graphene layers, which reduces the effective thermal pathways and consequently leads to a drop in thermal conductance. However, if the elongation of the model has a more pronounced effect on thermal conductivity than the reduction in thermal pathways, the overall thermal conductance may still maintain an upward trend.
Notably, the thermal conductance remains within the range of 0.026 GW/(m2·K) and 0.033 GW/(m2·K) across all models, irrespective of doping or deformation mode. This stability arises because the applied strains remain below the critical threshold for effective medium behavior, thereby preventing significant fluctuations of thermal conductance. It is crucial to note that when the stretching ratio exceeds 20%, the dampening effect of PDMS doping on thermal conductivity diminishes. It results from the stochastic nature of thermal pathways fracturing during stretching, which introduces variable impacts on thermal resistance. Consequentially, models with different PDMS doping rates show increased intercomputation variability in thermal conductance, reflected in larger standard deviations across multiple simulations.
Heat-transfer mechanisms in GF/PDMS composite
VDOS provides insights into the vibration modes of a material at different frequencies, typically visualized as the distribution of acoustic frequency. The phonon density and the thermal conductivity are generally positively correlated. In other words, a higher peak in the VDOS spectrum signifies a larger phonon density, driving enhanced thermal transport.
Heat transfer in GF/PDMS composites occurs via phonons across a wide frequency range from 0 THz to 60 THz, where nearly all phonons act as heat carriers. Figure 6a shows the VDOS spectra of the composites with different doping rates under a 40% compression ratio. This deformation ensures robust connectivity between PDMS and graphene layers, facilitating heat transfer through PDMS molecular structures.
a VDOS plots for composites with varying PDMS doping rates under 50% compression; b VDOS plots for the composite with 5% PDMS doping rate under different compression ratios; c VDOS plots for the composite with 5% PDMS doping rate under different stretching ratios; d Relationship between phonon frequency functions and corresponding thermal conductivity of GF composites with different PDMS-doping rates.
In the HNEMD calculations, an uncompressed model was initially employed to conduct an in-depth analysis of the VDOS and phonon contributions in models with doping rates of 0%, 5%, and 10%. For clarity in the analysis, phonons are categorized into three distinct frequency ranges: the relatively low-frequency (0 ~ 13 THz), the relatively mid-frequency (13 ~ 26 THz), and the relatively high-frequency (26 ~ 60 THz). The integral of the phonon frequency function within specified limits quantifies its contribution to the thermal conductance. Figure 6b depicts the contributions of phonons in different frequency ranges at 300 K, where GF/PDMS composites exhibit a wide range of phonon characteristics. Results are discussed collectively as both Fig. 6a, d correspond to models with identical heights (300 Å). For the pristine GF, conspicuous peaks in the spectral thermal conductivity appear across all frequency domains, indicating that phonons across the entire spectrum drive thermal transport. Nevertheless, in the GF doped with 5% PDMS, the contributions of both low- and high-frequency phonons to the overall thermal conductivity exhibit a diminishing trend. The decline intensifies with higher doping, as PDMS introduces interfacial thermal resistance that impedes heat flux and increases composite structural complexity, enhancing phonon scattering and reducing transport efficiency. The synergistic effects lower the composite’s overall thermal conductivity.
The VDOS analysis provides strong evidence to corroborate this trend. As the doping rate increases, the VDOS curve shifts downward overall, with the most pronounced changes in high- and low-frequency regions. Notably, the contribution of low-frequency phonons decreases continuously, whereas the high-frequency phonons initially experience a slight reduction upon PDMS addition but stabilize above 5% doping. The suppression of low-frequency phonons is identified as the dominant factor contributing to the reduction in thermal conductivity of the composite with the increasing PDMS content. This phenomenon can be elucidated by the fact that, although the changes in high-frequency phonons are more significant in the VDOS spectrum, the low-frequency phonons, owing to their higher group velocity, possess stronger thermal conductivity capabilities. Consequently, the low-frequency phonons play a more pivotal role in the variation of thermal conductivity.
Deformation-induced changes in model dimensions profoundly affect phonon contributions to thermal conductivity in GF/PDMS composites. During compression, reduced internal spacing suppresses high-frequency phonon propagation due to dense packing and increased interatomic forces, which enhance scattering. Conversely, new ligament connections formed under compression create additional thermal pathways for low- and mid-frequency phonons, partially offsetting the heat transport losses. However, interfacial thermal resistances from PDMS limit the overall recovery of thermal conductivity.
During stretching, expanded internal space activates high-frequency phonons, while diminishing the contribution of low-frequency phonons. Reduced interatomic forces decrease phonon scattering, favoring high-frequency transport. Considering the effects of group velocity, increased high-frequency phonons slightly outweigh decreased low-frequency phonons, leading to a modest rise in thermal conductivity. These dimensional-dependent phonon dynamics highlight the interplay between structural geometry and thermal transport in deformable composites.
Under cyclic tensile-compressive deformation, the pristine GF exhibits a thermal conductivity decrease from 1.18 W/(m·K) at 40% tensile strain to 0.49 W/(m·K) at 50% compression. In contrast, the thermal conductivity of the GF doped with 10% PDMS shows a smaller decline from 1.02 W/(m·K) to 0.43 W/(m·K). Across all doping concentrations, the changes in thermal conductivity are modest relative to structure deformations, ranging from 2.38 to 2.55 times. Thermal resistance calculations demonstrate that within the elastic range, the 5% PDMS-doped model undergoes the largest change in thermal resistance, increasing by a factor of 7.13. It highlights its promise for thermal tunable materials. While higher PDMS doping could further improve the mechanical robustness, computational limitations restrict the maximum simulated compression to 50%, lower than the real-world scenarios. This implies that a greater thermal resistance modulation could be achieved in practical applications, reinforcing the utility of GF/PDMS composites for adaptive thermal management.
