Transferable machine learning model for multi-target nanoscale simulations in hydrogen-carbon system from crystal to amorphous

Abstract

Carbon materials especially with hydrogenation have attracted wide attention for their novel physical and chemical properties and broad application prospects. A systematic theoretical simulation method accurately describing atomic interactions for hydrogen-carbon systems is crucial for the design of carbon-based materials and their industrial applications. Multiphases of hydrogenated carbon materials, from crystal to amorphous, with covalent network and diverse chemical reactions bring huge difficulties to construct a general interatomic potential under various conditions. Here, we demonstrate a transferable active machine learning scheme with separated training of sub-feature spaces and target-oriented finetuning, and construct a general-purpose pre-trained machine learning potential (MLP) for hydrogen-carbon systems. The pre-trained MLP is further efficiently transferred to three target spaces of deposition, friction and fracture with scale reliability. This work provides a robust tool for the theoretical research of hydrogen-carbon systems and a general scheme for developing transferable MLPs in multiphase systems across compositional and conditional complexity.

Introduction

Carbon materials, which include diverse morphologies such as graphene, diamond, carbon nanotubes, amorphous carbon^{1, 2, 3–4}, etc., have been widely investigated for their potential applications in various fields spanning electronic devices, energy storage, lubrication and wear resistance, owing to the intrinsic electrical conductivity, thermal stability, mechanical strength and tribological property^{5, 6, 7–8}. Hydrogen, the most abundant element in nature, has been found to significantly affect the performance of carbon materials^{9, 10–11}. Through hydrogenation modification, the superior electrical and tribological properties of carbon materials can be achieved with the alterable band-gap⁹ and electrostatic repulsion of hydrogen¹⁰, while the hydrogen embrittlement can also deprave the adhesion strength and wear resistance of carbon coating¹¹. Establishing a systematic research method for hydrogen-carbon system is crucial for the investigation of optimized design of carbon-based materials and their industrial applications.

Carbon materials generally present in various forms, ranging from crystal to amorphous. The rehybridization of the carbon atoms in the covalent network from sp¹ to sp³, could lead to the transition between different carbon materials under certain conditions, e.g., graphite-to-diamond transformation with coherent interfaces during the recovery from static compression², shear induced mechanical amorphization of diamond¹², graphitized transformation of amorphous carbon (a-C) film during friction process^8,13, etc. In particular, the amorphization behavior in carbon materials has been widely observed under multiple conditions, e.g., defect existence¹⁴, fracture and healing¹⁵, friction and wear^12,16, etc., and can significantly influence the electrical, mechanical and tribological properties of carbon materials. Accurately capturing the structural evolution and hybrid bonding transformation of amorphous carbon is crucial to investigate the physical and chemical properties of carbon materials. However, it still meets challenges as the existing characterization methods are difficult to quantitatively identify the amorphous nanostructure evolution¹⁷.

Theoretical simulations, which can provide atomic-level structural information, have been widely used in performance research and material design of carbon systems. Ab initio molecular dynamics (MD) simulations, with quantum-mechanical calculation of the electronic structure, can accurately describe atomic interactions and chemical bonding reactions in carbon systems^{18, 19, 20–21}. Nevertheless, the computational expense and unfavorable scaling inherent to the ab initio method preclude its use for complex, large-scale dynamic processes. Instead, empirical potentials are employed, which typically exhibit a functional form based on physical/chemical intuition and convenience. However, these potentials are not always sufficiently flexible to accommodate multiple properties simultaneously²², and they are limited in accuracy and transferability with respect to phase transitions, defective structures, and stress-induced structure evolution in carbon materials^21,23,24.

Machine learning (ML) methods, which are capable of fitting functions in high-dimensional spaces without specific parameterized functional forms²², have been introduced to construct interatomic potential with density functional theory (DFT) level accuracy for MD simulations^{25, 26, 27–28}. Many efforts have been made to develop the machine learning interatomic potential in carbon materials for their complicated structural responses under various environmental conditions^{27, 28, 29–30}, e.g., growth mechanisms of graphene and carbon nanotube on metal surface^27,28, graphite–diamond phase transition²⁹, deposition process of tetrahedral amorphous carbon³⁰. A general machine learning potential for amorphous carbon³¹ was developed at 2017 which was further extended to describe numerous carbon allotropes³² at 2020. Recently, some efforts were made to describe the atomic interactions for hydrogen-carbon systems^33,34, including hydrocarbon molecules, even to hydrogenated amorphous carbon. However, a general interatomic potential aiming at larger scale and various conditions with accurate description of the atomic interactions for hydrogen-carbon systems, including equilibrium and non-equilibrium states³⁵, still needs to be developed systematically.

The development of MLPs hinges on the fundamental aspect: generation of an efficient training dataset that encompasses the full structural and chemical space³⁶. For the hydrogen-carbon systems with multiphases and phase transition processes, the comprehensive coverage of structural and chemical space, including the complex hybridization transitions in covalent networks and active hydrogen chemical bonding reactions³⁷, meets a huge challenge. Especially in the amorphous carbon systems, the atomic chemical environment of each atom can exhibit significant variation. However, considering high-quality ab initio calculations with high computational cost, building a training dataset through exhaustive exploration of the chemical space is intractable.

Active learning, which aims at construction of general-purpose models with minimal scope of training data, can help to select the most informative data from unlabeled structures and provides a possibility to construct the MLPs with more efficient and robust dataset^38,39. Uncertainty quantification (UQ) is a centerpiece for active learning to mitigate epistemic uncertainty at extrapolative region, especially for the sampling of some rare events^{39, 40–41}. Many UQ schemes have been utilized for the facilitation of active learning, e.g., Bayesian neural networks (BNNs)⁴², Monte Carlo dropout⁴³, D-optimality⁴⁴, and neural networks (NN) ensemble³⁸, among which, BNNs and NN ensemble have been widely recognized and promoted to the sampling of rare events with enhanced sampling methods^45,46, even though these methods could still lead to high computational cost for large and complex systems⁴¹.

In this work, we demonstrate a transferable active machine learning scheme for the hydrogen-carbon systems, covering the multiphases from crystal to amorphous structures, with the strategy of separated pre-training of the sub-feature spaces and the target-oriented finetuning. The dataset with the dimensionality reduction of atomic environments was separated to multiple chemical spaces to improve the training efficiency and transferability of MLP model. The transferability and reliability of the model are inspected under three different conditions: growth of hydrogenated amorphous carbon (a-C:H) film on diamond substrate, shear-induced hybrid bonding transformation and structure distortion in a-C:H during friction process and fracture behavior of crystal and amorphous carbon, all of which showing a DFT-level accurate description of atomic interactions. This transferable active machine learning model can be easily applied to any other conditions in hydrogen-carbon systems with selective target-oriented finetuning. This work provides a general scheme for the development of transferable MLPs across any compositional and conditional complexity.

Results

Transferable machine learning model with pre-training and target-oriented optimization

The training strategy of our transferable MLP model for hydrogen-carbon systems is described below and schematically depicted in Fig. 1. The whole training procedure involved a transferable pre-training process and a finetuning optimization at specific target space, and the training processes were carried out with the neural network architecture via DeePMD-kit package⁴⁷ (see “Methods”). For the pre-train process, the feature space of hydrogen-carbon systems was separated to three major sub-feature spaces: a-C, a-C:H and crystal space, based on the atomic environmental distribution of data structures (further discussion in Fig. 2) to reduce the interference of complex structure information, improving the efficiency of dataset construction (Supplementary Note 1). The initial structures in the sub-feature space were obtained from Material Project database and the melt quenching process in ab initio molecular dynamics (AIMD) simulations (Supplementary Note 2). In each sub-feature space, an independent two-stage active learning^38,48 workflow with NN ensemble for dataset expansion was conducted with the simulation scale transferring from the order of magnitude 10² to 10⁴ atoms to increase the reliability of MLPs over different length scales. During the parallel active learning process, the MLP-driven MD simulations were conducted under various conditions, e.g., extreme temperatures (from 0.1 to 5000 K), large-range pressures (from 1 to 15 GPa) and structure deformation process (Supplementary Table 1), and the enhanced sampling method⁴⁹ was also utilized to capture the transition state of chemical reactions (e.g., hydrogen transfer reactions) (see “Methods” for more details). Notably, the secondary sub-feature spaces, in which the structure features are partial intersection, can also be utilized for enhanced training of specific atomic environment, which, in this work, were the bulk, slab and homogeneous/heterogeneous interface structures.

Fig. 1 Training strategy of transferable MLP for hydrogen-carbon system and target-oriented finetuning. [Images not available. See PDF.]

a Schematic visualization of the training strategy of our transferable MLP for hydrogen-carbon system, which concluded a pre-training process for separated sub-feature spaces and the target-oriented finetuning. b Schematic diagram of convergence process of MLP in structure space during active learning workflow. There are three major sub-feature spaces in this work: a-C, a-C:H and crystal space with the data points represented in different shapes in Fig. 1b, while the secondary sub-feature spaces, e.g., bulk phase and different surface states, are not shown in this schematic diagram.

Although the pre-trained MLP model could accurately describe the atomic interactions from crystal to amorphous hydrogen-carbon structures, majorly in the stable state, some complex physical and chemical processes, which might include local high stress state or complex stress field distribution and lead to the distortion of hybrid covalent network and active chemical bonding reactions in hydrogen-carbon systems⁵⁰, could be beyond the describing ability of the pre-trained model, and the specific model optimization for target space with active learning method was introduced. In this work, we considered three complex physical and chemical processes: deposition, friction and fracture, to evaluate the transferability and scale reliability (from small to large) of our MLP model. Basing on the pre-trained MLP model, the similar two-stage active learning (see “Methods”) was also conducted at each target space. During the iterative active learning, the target processes with various physical and chemical reactions, which could be the rare events for the pre-training dataset, were performed in MD simulations with MLP model (more details in “Methods”). For the deposition target space, considering the influence of substrate scale during the depositing, the melt quenching process, which was widely utilized for emulating the diffusion process of local thermal spike during the deposition^{51, 52–53}, was conducted in small-scale active learning rather than the deposition process, and for large-scale MD simulations, the deposition and melt quenching process were both carried out for adequate sampling.

For the local high stress state or complex stress field distribution in large-scale simulations, the high error region and state for pre-trained MLP model could emerge only at certain local structures and timepoints, and the high throughput calculations for large-scale structures filtered from active learning method were inefficient and computational costly. Here we introduced an effective local configuration identification method during active learning in large-scale simulations, in which the local high deviation regions were identified and extracted in a large system for labeling (more details in “Methods”). Basing on the uncertainty quantification of NN ensemble^38,48, these local configurations with high atomic force deviation were beyond the capability of the MLP model and could contained the unlearned atomic environment information, which therefore were the efficient sampling for the MLP model.

The convergence process of MLP model during the active learning workflow could be represented by a schematic diagram, as shown in Fig. 1b. During the pre-training process with parallel active learning, the sufficient atomic environment information in each sub-feature space was collected efficiently until the MLP model could accurately describe the subspace independently (colored rectangle in Fig. 1b). When the sub-feature spaces were merged to the full feature space (Supplementary Note 3), the interference of different subspaces to model training should be finetuned with limited active learning finetuning to achieve high accuracy description in full feature space. However, there were still other spaces with high uncertainty (shadow region in Fig. 1b) for pre-trained MLP, e.g., some complex physical and chemical processes, which might include local high stress state or complex stress field distribution. These rare events for pre-training dataset would be further finetuned at target space with the two-stage active learning until the target process could be well described.

The final datasets for pre-training process and three different target spaces were shown in Supplementary Table 1, which contained 29,885 structures with ~9.6 million atoms totally, and was randomly divided into the training set (90%) and the testing set (10%) for the model training and validation test, respectively. Further, a similarity metric of atomic environment, Smooth Overlap of Atomic Positions (SOAP)^54,55, and a non-linear dimensionality reduction technique, Uniform Manifold Approximation and Projection (UMAP)⁵⁶, were used to visualize the chemical diversity of atomic environments in our datasets, which we called the SOAP-based UMAP embedding method in the following text, as shown in Fig. 2. About 300 structures in pre-training and target-space datasets including more than 45,000 atomic environments were randomly selected to create the map, in which each point represented an atomic environment, and different colors in Fig. 2 represented the different categories of structures in the insets, and the atomic energy and relative energy error of these atomic environments were shown in the Supplementary Fig. 3. As shown in SOAP-based UMAP map, the atomic environments from target space of friction, deposition, and fracture were partly covered by the sub-feature datasets, which indicated that the pre-trained model was general-purpose and could be applied to various applications. Notably, the atomic environments of amorphous crack structures were also covered by the crystal crack structures, indicating that the information learned from crystal structures could be transferred to related amorphous structures.

Fig. 2 SOAP-based UMAP map of atomic environments in the datasets. [Images not available. See PDF.]

Points in this map were the two-dimensional UMAP embedding of high-dimensional SOAP vectors to visualize and compare the atomic environments in our datasets, colored by different categories of the structures, and their related typical configurations were inserted in the gray circles. The typical configurations in target spaces, which were deposition, friction, and fracture space, were highlighted in the gray squares, represented in color of scarlet, yellow, and orange. In all the snapshots inserted, the gray atoms represented the carbon atoms and the light blue atoms represented the hydrogen atoms.

Validation of MLP models: from crystal to amorphous

Firstly, the validation of the pre-trained model was verified in calculating the atomic interactions in crystal hydrogenated carbon systems, as shown in Fig. 3. The bonding energy curves of hydrogen atom at surface of graphene and diamond (Fig. 3a, b) were calculated, showing high accuracy near the equilibrium position compared to DFT method, while a slight deviation at extreme short-distance⁵⁷ or far away from equilibrium position which needed to be further optimized. Notably, the test configurations for bonding energy of hydrogen atom were not included in the pre-training dataset, indicating the transferability of the pre-trained MLP model. Calculations of energy barriers of hydrogen transfer process and van der Waals interactions between the layers of graphite and hydrogenated graphite in AB stacking mode were also highly consistent with the DFT results with the average error of the energy per atom to be 2.2, 2.3, 5.2 and 2.1 meV/atom, respectively (Fig. 3c–e, g). The total energy change versus the distance between created surfaces was calculated by separating the ideal crystal structure to create new diamond (111) surfaces, which also showed a high consistency between MLP results and DFT calculations with the average error of the energy per atom to be 21.0 meV/atom (Fig. 3f). The MLP results were also compared with other methods and the numerical error were listed (Supplementary Fig. 4, Supplementary Table 3). Further, we applied the pre-trained MLP model to the shear-induced diamond-to-graphite transfer process (Supplementary Note 4, Supplementary Fig. 5), showing the same structure transformation as the DFT calculations²¹. The critical shear stresses for the transition process in our work also showed the higher consistency with the DFT results compared with the empirical potentials³⁰ (Supplementary Table 4).

Fig. 3 Validation of the pre-trained MLP model for atomic interactions in crystal hydrogenated carbon systems. [Images not available. See PDF.]

Bonding energy curves of (a) hydrogen atom at graphene surface and (b) hydrogen atom at diamond (100) surface. Energy barriers for (c) hydrogen atom transferring along armchair at graphene surface with distance of 1.1 Å and (d) hydrogen atom transferring along [010] direction at diamond (100) surface with distance of 1.1 Å. The total energy change versus separation distance of (e) layers of graphite in AB stacking mode, (f) decohesion simulations for diamond (111) surface, (g) layers of hydrogenated graphite in AB stacking mode and hydrogenated diamond (111) surfaces. The ‘D’ in (a) and (b) represented the distance between the atomic centers of hydrogen and carbon atoms, in (e), (g), (h) represented the minimum distance between carbon atoms in the out-of-plane direction between adjacent atomic layers, and in (f) represented the separated distance between the created surfaces. The ‘d’ in (c) and (d) represented the distance of atomic centers of hydrogen atom moving in the in-plane direction. The snapshots for different processes were inserted with the gray and blue spheres representing the carbon and hydrogen atoms. The ‘Distance’ axis represented the ‘D’ perpendicular to the surface in (a), (b), (e–h), and the ‘d’ parallel to the surface in (c, d). The Energy error represented the absolute error between DFT and MLP calculations.

The reliability of pre-trained MLP model was further verified for hydrogenated amorphous carbon systems. The properties of amorphous carbon, including the Young’s modulus, bulk modulus and shear modulus, were calculated over a wide scaling of hydrogen content and density (Supplementary Table 5), which also showed the DFT level accuracy, indicating the general adaptability of the pre-trained MLP model in hydrogen-carbon systems, from crystal to amorphous structures. Besides, the same calculations were also conducted in large scale with MLP model to confirm its scale reliability, as shown in Table 1.

Table 1. Properties of amorphous carbon with different hydrogen contents in density of 2.4 gcm⁻³

Amorphous Carbon	Young’s modulus (GPa)			Bulk modulus (GPa)			Shear modulus (GPa)
H content (atomic%)	DFT	MLP (Small)	MLP (Large)	DFT	MLP (Small)	MLP (Large)	DFT	MLP (Small)	MLP (Large)
0	291.71 (7.04)	298.03 (15.97)	291.72 (7.04)	193.76 (5.81)	215.83 (22.57)	199.78 (11.98)	115.57 (6.46)	118.82 (3.37)	118.29 (3.41)
10	281.27 (3.23)	262.40 (10.46)	278.42 (13.22)	216.91 (4.15)	198.63 (14.21)	193.01 (6.72)	112.11 (4.74)	102.17 (6.78)	106.36 (9.23)
20	267.97 (3.05)	254.69 (15.74)	240.76 (10.15)	199.03 (4.02)	199.47 (8.48)	191.64 (1.53)	106.52 (3.12)	108.56 (2.85)	102.57 (10.70)
30	216.43 (4.76)	196.67 (23.19)	194.63 (13.48)	136.62 (5.15)	130.44 (7.24)	127.15 (9.57)	87.71 (1.85)	85.59 (3.50)	81.22 (6.76)
40	231.13 (3.59)	238.98 (14.69)	234.01 (12.78)	198.47 (4.76)	196.26 (0.84)	195.16 (3.87)	90.17 (3.83)	90.37 (5.60)	87.76 (7.38)

The structures of the a-C:H with different hydrogen contents were constructed through melt quenching process (see “Methods”) in AIMD simulations and MLP-based MD simulations, which, respectively, were calculated with DFT methods and pre-trained MLP model for their modulus. The structures for DFT and MLP (Small) contained the order of magnitude 10² atoms and that for MLP (Large) contained the order of magnitude 10⁴ atoms. The numbers in parentheses are the standard deviation of the corresponding values calculated in three independent MLP-based MD simulations and AIMD simulations.

Although the pre-trained MLP model showed the general reliability of the hydrogen-carbon systems, some complex physical and chemical processes including local high stress state or complex stress field distribution could be beyond the describing ability of the pre-trained MLP model, for example, the deposition, friction and fracture processes. After active learning finetuning of the pre-trained MLP model at target spaces, the MLP models presented highest accuracy of the corresponding test datasets (Table 2 and Supplementary Fig. 6 and Supplementary Note 5). Notably, the target-oriented finetuning processes could be achieved based on the target datasets with an order of magnitude fewer structures than pre-training datasets for all three different target conditions (Supplementary Table 1), indicating the high transferability of pre-trained MLP model and high efficiency of the target-oriented finetuning process based on the active learning methods. The accuracy of MLP models in target spaces were also evaluated compared with the empirical potentials (Table 3 and Supplementary Note 6), and the relative errors were provided (Supplementary Table 6). Since the Tersoff’s empirical bond order potential is not suitable for deposition of hydrocarbons molecules (C₂H₂ and CH₄ used in this work), the corresponding value in Table 3 is vacant. The MLP models with best performance in three target spaces were tested, respectively, and the orders of magnitude improvement of accuracy were investigated with the computational efficiency comparable with the empirical potentials (Supplementary Fig. 7).

Table 2. Validation of the pre-trained MLP model and the target-oriented MLP model

System		Pre-trained model		MLP for target space
				Deposition		Friction		Fracture
		∆E	∆F	∆E	∆F	∆E	∆F	∆E	∆F
Sub-feature space	a-C-bulk	0.032 (0.010)	0.63 (0.39)	0.023 (0.007)	0.71 (0.45)	0.020 (0.006)	0.32 (0.13)	0.027 (0.031)	0.66 (0.60)
	a-C-slab	0.031 (0.010)	0.52 (0.34)	0.029 (0.010)	0.72 (0.42)	0.015 (0.004)	0.36 (0.20)	0.042 (0.018)	0.48 (0.43)
	a-C:H-bulk	0.015 (0.007)	0.29 (0.19)	0.013 (0.007)	0.33 (0.25)	0.010 (0.003)	0.22 (0.11)	0.021 (0.010)	0.34 (0.30)
	a-C:H-slab	0.029 (0.016)	0.4 (0.22)	0.038 (0.015)	0.42 (0.31)	0.026 (0.009)	0.26 (0.15)	0.038 (0.029)	0.38 (0.34)
	Crystal	0.027 (0.018)	0.27 (0.21)	0.208 (0.169)	0.55 (0.34)	0.015 (0.007)	0.24 (0.11)	0.027 (0.015)	0.32 (0.25)
Target feature pace	Deposition	0.115 (0.055)	4.66 (4.36)	0.019 (0.007)	0.78 (0.44)	0.098 (0.035)	1.35 (1.41)	0.047 (0.018)	0.86 (0.66)
	Quench	0.110 (0.073)	1.13 (0.61)	0.021 (0.010)	0.93 (0.54)	0.025 (0.010)	0.88 (0.54)	0.050 (0.033)	0.91 (0.60)
	Shear	0.027 (0.019)	0.64 (0.39)	0.085 (0.041)	1.21 (1.03)	0.008 (0.004)	0.50 (0.31)	0.025 (0.015)	0.77 (0.51)
	Friction	0.110 (0.058)	2.85 (1.55)	0.097 (0.046)	1.55 (0.55)	0.017 (0.010)	0.42 (0.29)	0.059 (0.025)	0.79 (0.57)
	Crack-tip	0.106 (0.070)	0.80 (0.62)	0.173 (0.073)	0.76 (0.43)	0.172 (0.093)	0.92 (0.62)	0.034 (0.020)	0.54 (0.52)

The energy (ΔE) and force (ΔF) root mean square error (RMSE) values compared with DFT calculations in eV/atom and eV/Å were collected on different test datasets at sub-feature space and target feature space. The numbers in parentheses are the standard deviation of the ΔE and ΔF in three times random sampling of the test datasets. The bold numbers represented the highest accuracy of the corresponding test datasets.

Table 3. Performance of the target-oriented MLP models compared with other potentials

Target feature space	Empirical potentials										GAP		MLP (Best performance)
	REBO		Tersoff		REBO-s		Tersoff-s		Reaxff
	∆E	∆F	∆E	∆F	∆E	∆F	∆E	∆F	∆E	∆F	∆E	∆F	∆E	∆F
Deposition	1.114 (0.35)	12.71 (22.3)	—	—	1.320 (0.316)	27.33 (26.53)	—	—	0.127 (0.078)	7.64 (6.61)	—	—	0.019 (0.007)	0.78 (0.44)
Quench	0.990 (0.140)	5.84 (4.89)	1.617 (0.354)	6.97 (5.37)	1.719 (0.656)	18.35 (16.09)	2.186 (0.479)	7.88 (5.37)	0.403 (0.289)	6.07 (4.49)	0.198 (0.033)	0.68 (0.58)	0.021 (0.010)	0.93 (0.54)
Friction	0.816 (0.273)	2.81 (2.11)	1.159 (0.051)	5.33 (4.14)	0.921 (0.293)	3.76 (13.36)	1.502 (0.045)	7.48 (5.41)	0.138 (0.095)	5.10 (3.33)	0.135 (0.021)	4.43 (0.52)	0.017 (0.010)	0.42(0.29)
Crack-tip	0.649 (0.048)	2.03 (1.80)	0.757 (0.213)	2.93 (2.584)	0.606 (0.046)	2.98 (2.65)	0.775 (0.241)	3.33 (2.96)	0.315 (0.113)	4.56 (3.84)	0.093 (0.028)	1.38 (2.58)	0.034 (0.020)	0.54 (0.52)

The energy (ΔE) and force (ΔF) root mean square error (RMSE) values compared with DFT calculations in eV/atom and eV/Å were collected on different test datasets at target feature space. The numbers in parentheses are the standard deviation of the ΔE and ΔF in three times random sampling of the test datasets. The bold numbers represented the highest accuracy of target-oriented MLP models corresponding to different target datasets.

Performance in multi-target conditions and scale reliability

Based on the target-oriented MLP models, we carried out the MD simulations for the three target conditions to evaluate their performance and scale reliability. Firstly, the well-trained MLP model was used for the construction of a-C:H structure via both melt quenching and the deposition process (see “Methods”), as shown in Fig. 4. The radial distribution functions (RDFs) of C–C and C–H in both melting liquid structures and amorphous carbon structures with atoms from 200 to 10,000 showed good agreement with the AIMD simulations (Supplementary Fig. 8), indicating the high reliability of MLP models in large scale. Further, the melt quenching process was carried out in large scale (10,000 atoms) to construct the a-C:H with a series of hydrogen content, which also showed consistent RDFs with the structures in small scale AIMD simulations (Fig. 4a). The simulation setting of MLP based MD and AIMD In this way, the ideal a-C:H structures could be constructed in any required scale with our MLP model, overcoming the multi-periodic repeating of the a-C:H in small cell from AIMD simulations. As the a-C:H could be used as the protective and lubricating coating^8,20, the growth of a-C:H film on diamond (100) from deposition process was also carried out with the MLP model (see “Methods”). With a gradient area of structure between the substrate and surface, the stable growth of a-C:H was achieved with the almost constant sp³ content, relative density and hydrogen content (shadow area in Fig. 4b), which was also observed in the growth of ta-C with a machine learning potential developed for carbon-only system³⁰. The stable growth structure of a-C:H also showed the DFT level accuracy with the consistent RDFs with the structures in AIMD simulations, as shown in Fig. 4c.

Fig. 4 Properties of hydrogenated amorphous carbon structures constructed with target-oriented MLP model in large-scale. [Images not available. See PDF.]

a Radial distribution functions of C–C and C–H in a-C:H with different hydrogen contents constructed with the MLP model for deposition target space and AIMD simulations, respectively, in the same melt quenching process except for timesteps. The RDF was calculated from the configurations in NVT simulations under 300 K for 3 ps. b Atomic configuration of a-C:H film constructed through the deposition process with the MLP model and its properties along the growth direction of film. The gray and blue spheres representing the carbon and hydrogen atoms. The density is the relative value normalized by diamond. c Radial distribution functions of C–C and C–H in a-C:H constructed in AIMD melt quenching process and MLP based deposition process.

In order to evaluate the capability of the target-oriented MLP model to describe the passivated interface in stable friction and the hybrid bonding transformation during the running-in process in hydrogen-carbon systems, the friction tests of hydrogen passivated diamond (111) surface and a-C:H with 30 at.%H were conducted with the MLP model, as shown in Fig. 5. Firstly, the sliding potential surface energy (PES) corrugations ΔE of hydrogen passivated diamond (111) interface with respect to the surface stacking mode with minimum energy was calculated, as well as the corresponding lateral force (Supplementary Note 7). A highly consistent PES corrugations during sliding was observed between MLP and DFT results with the average error of 0.05 meV/Ų (Fig. 5b and Supplementary Fig. 9), as well as the lateral force with the average error of 0.24 meV/ų. The interfacial adhesion of the hydrogen passivated diamond (111) interface was calculated to be 0.37 J/m², consistent with the DFT results and the previous work⁵⁸, showing a better performance than other methods (Supplementary Table 7).

Fig. 5 Friction test of hydrogen passivated diamond and a-C:H with 30 at.%H with target-oriented MLP model. [Images not available. See PDF.]

a, b Frictional energy landscape of hydrogen passivated diamond (111) interface. a Sliding PES corrugations ΔE calculated by DFT. b The potential profile (ΔE, red) and lateral force (F, blue) along the sliding pathways calculated by DFT and MLP. The step lengths along x ([ $1\bar{1}0$ ]) and y ([ $11\bar{2}$ ]) axis for the sliding systems are 0.25 Å and 0.435 Å. The inserted snapshots were the bidimensional representations of hydrogen passivated diamond (111) interface and the sliding directions, respectively. The gray and blue spheres representing the carbon and hydrogen atoms. c–e Typical bonding reorganization in friction process of a-C:H with 30 at.%H with target-oriented MLP model. c Snapshots of typical bonding reorganization in MD simulations of friction process with MLP model in friction target space. The black arrows at bottom left represented the direction of x and z axis in MD simulation. The gray and blue spheres represent the carbon and hydrogen atoms. d Change of the hybridization of C₃ in the four bonding states of (c). The yellow and green square represented the sp² bond and sp³ bond area. The MBA and SDBA represented the mean of bond angles and the standard deviation of bond angles (red and blue dot, respectively). e Frequency of sp³ jumping during the friction process.

Then, the MLP based MD simulation was carried out for the friction process of a-C:H with 30 at.%H (see “Methods”), and a typical bonding reorganization process during the friction was extracted with four bonding states, as shown in Fig. 5c–e. As the reliability of dynamic evolution during the friction process was confirmed with the highly consistent evolution curve of energy and friction force with DFT calculations (Supplementary Fig. 10) and the chemical desorption of hydrocarbons on surface dangling bonds could be well described (Supplementary Fig. 11), we focused on the capability of MLP model to describe the local bonding reorganization and the formation and release of distorted structures during the friction.

As shown in Fig. 5c, the shear-induced C–C bond stretch and twist was observed, followed with the bonding reorganization. The bond angles between central carbon and its surrounding atoms were calculated and the mean value and standard deviation of bond angles for the four target carbon atoms were extracted in four bonding states to investigate the hybrid bonding transformation and the bonding distortion, respectively. As shown in Fig. 5d and Supplementary Fig. 12, the C₁ and C₂ atom could maintain the sp² and sp³ hybridization states, respectively, however, with the shear-induced bonding distortion, while the hybrid bonding transformation was observed for C₃ and C₄ with the release of the distorted bonds. Notably, the hybrid bond shift jumping from sp²-sp³-sp² was observed (sp³ jumping bridge in Fig. 5d), which significantly released the sp² bonding twist, indicating the crucial role of hybrid bonding transformation for the release of distorted structures in carbon network. Further, the frequency of sp³ jumping during the friction process was calculated and the occurrence of the jumping event showed the dependence with the friction force during the sliding, as shown in Fig. 5e, indicating a stress-assisted hybrid bonding transformation at the sliding interface. The similar effect of bonding reorganization in a-C:H system including the hybridization transition and hydrogen transfer reactions was also observed in the MLP based large-scale MD simulations for shear-induced structure ordering transformation of a-C:H in our previous work⁵⁰, indicating the scale reliability of the MLP model.

Finally, the atomic-scale quasi-static fracture and static decohesion MD simulations were carried out with the MLP model for fracture target space (see “Methods”), as shown in Fig. 6a, c. During the propagation of diamond [001] crack front on the (110) plane and graphene in-plane crack front along the zigzag edge, the change of bond length on the site of crack-tip could be well described with DFT level accuracy. The cleavage behavior of other crystal carbon cracks was also well described, as shown in Supplementary Fig. 13, consistent with the previous works and DFT results, and the corresponding accuracy plotted in Supplementary Fig. 14 showed the error in the magnitude of 10⁻² eV/atom^23,26,59. The effect of hydrogen on the diamond [001] crack front on the (110) plane was also tested with acceptable accuracy (Supplementary Fig. 15), and it could be seen that the absorption of hydrogen on the crack edge lead to break of the bond at the crack-tip, which can assist the propagation of the crack. Physically, as the fracture was a process of creating new surfaces, decohesion simulations of diamond and graphene were performed for further validation^23,24. The total energy change versus the distance between created surfaces was calculated by separating the ideal crystal structure to create new diamond (110) surfaces or new graphene zigzag edges, and the cohesive stress was the first derivative of total energy curves normalized by the created surface area, which also showed a high consistency between MLP results and DFT calculations (Fig. 6b, d).

Fig. 6 Validation of the MLP model for fracture target space in crystal carbon systems and simulations of fracture toughness of a-C:H. [Images not available. See PDF.]

Bond length change versus normalized stress intensity factor by the value of fracture toughness ( ${K_{\mathrm{{\rm I}}}/K}_{\mathrm{{\rm I}}\mathrm{C}}$ ): (a) for diamond [001] crack front propagation on the (110) plane and (c) for graphene in-plane crack front propagation along the zigzag edge. The ‘L’ in (a) and (c) represented the distance between the atomic centers on the site of the crack-tip. The energy and force change versus separation distance extrapolated from decohesion simulations: (b) for diamond (110) and (d) for graphene zigzag edge. The ‘D’ in (b) and (d) represented the separated distance between the created surfaces or edges. The snapshots for the processes above were inserted with the gray spheres representing the carbon atoms, and the blue spheres representing hydrogen atoms. Fracture toughness ( $K_{\mathrm{{\rm I}}\mathrm{C}}$ ) extrapolated from the MD simulations: (e) for a-C with different density (2.0–3.2 gcm⁻³) and (f) for a-C:H with different hydrogen content (10–40 at.%H).

Different from the perfect cleavage and lattice-trapping effect of fracture of crystals, the fracture process zone of amorphous carbon was an ensemble of highly-stressed and distorted atomic environments, with cross-linked carbon networks containing hundreds of atoms at least, and the crack-tip propagated bond by bond. With the fracture-targeted model, we performed fracture simulations of a-C with different density and a-C:H with different hydrogen content to evaluate the fracture toughness (Supplementary Fig. 16). As shown in Fig. 6e, the fracture toughness of a-C increased with density, while for a-C:H, it decreased with hydrogen content. Four specimens with different amorphous structure were used in our simulations for each density and H content. The R curves of each specimen were shown in Supplementary Fig. 17 (Supplementary Note 8), which represent the amplitude of load (i.e., stress intensity factor) versus change of crack length, and the vertical coordinate of the inflection point when the value of crack length turns from zero to positive is the critical stress intensity factor ( $K_{\mathrm{{\rm I}}\mathrm{C}}$ ), i.e., fracture toughness of the specimen⁶⁰. The range of calculated fracture toughness was in typical range of reported value: 0.8–4.4 $\mathrm{MPa}\sqrt{\mathrm{m}}$ ^{61, 62–63}, measured by micro-tests of a-C:H. Toughness of hydrogenated structures was significantly lower than pure carbon, indicating that hydrogenation would cause brittleness of a-C^11,64, which also needed further investigation.

Discussion

In summary, a transferable machine learning model for the hydrogen-carbon systems from crystal and amorphous structures was demonstrated. The separated training of the sub-feature spaces followed with a target-oriented finetuning was introduced accompanied by the dimensionality reduction of atomic environments for the efficient generation of training datasets with a multi-scale active learning scheme. The training paradigm of the transferable active machine learning model can be concluded as follows: (1) Clarify the structures of the research system and construct the initial datasets. The initial structures can be from the existing datasets, e.g., Material Project database, or a standard scheme for structure construction; (2) Dimensionality reduction of atomic environments in initial datasets and dividing the sub-feature spaces; (3) Training and iterative optimization of the MLP in sub-feature spaces separately with multi-scale active learning scheme; (4) Merging the sub-feature datasets based on the target spaces and building the pre-trained MLP model; (5) Training and iterative optimization of the MLP in target spaces based on the merged pre-training datasets with the target-oriented active learning. This transferable active machine learning paradigm provides a scheme for the development of MLP in any other systems across compositional and conditional complexity.

The transferability and scale reliability of the model was further evaluated at deposition, friction and fracture target spaces, showing the excellent performance. The ideal a-C:H structures could be constructed in any required scale with the deposition-target MLP model, overcoming the multi-periodic repeating of the a-C:H in small cell from AIMD simulations. The shear-induced complex hybrid bonding transformation and structure distortion in carbon network could also be well described with our MLP model as confirmed in friction test. Finally, the cleavage behavior in crystal carbon fracture could be correctly described with our MLP model and the fracture toughness of a-C(:H) was consistent with experimental results. This transferable MLP provides a general model for hydrogen-carbon systems and can be easily expanded to any other target conditions, providing a powerful tool for the theoretical research and atomic understanding of hydrogen-carbon systems.

Methods

Construction of pre-training dataset

The pre-training dataset was constructed with three parts: AIMD based initial dataset, data from active learning in small scale and data from active learning in large scale. The initial dataset (gray part in Fig. 1a) for the zero iteration MLP model before the sampling process with two-stage active learning was constructed under various conditions (e.g., temperature, pressure, cell scaling and shifting) with AIMD simulations (see Supplementary Table 1 for detail setting), and the iterative active learning driven by MLP was further conducted. During the first-stage active learning, the a-C:H structures in small scale (order of magnitude 10² atoms) were considered in the same range of temperature and pressure with NVT and NPT ensemble. The NN ensemble (Strategy of uncertainty quantification in “Methods” for more details) was utilized for the uncertainty quantification of the active learning. The biased MD simulation was performed only in the regions where the uncertainty indicator was low⁴⁹. The configurations with high deviation of atomic forces in NN ensemble were labeled directly in DFT calculations considering the acceptable computational cost, and then added to the dataset (Refine data in Fig. 1a). For the second-stage active learning, the structures were expanded to the order of magnitude 10⁴ atoms (Scale expansion in Fig. 1a) and only the high deviation regions of atomic forces in a whole configuration were extracted to be labeled in DFT calculations during the active learning. The strategy of labeling region extraction in large-scale active learning would be discussed in the next part.

Strategy of local configuration extraction for large-scale active learning

Taking the friction process as an example, we introduced the strategy of local configuration extraction for large-scale active learning, as shown in Fig. 7. The large-scale MD simulations for friction processes of crystal and a-C(:H) with various densities and hydrogen contents were conducted with the pre-trained MLP model, and the atomic forces in the configurations during the friction were calculated for uncertainty quantification (Strategy of uncertainty quantification in “Methods” for more details). The deviation of atomic force (DAF) between three parallel MLP models were calculated and mapped to the configurations during the friction, as shown in Fig. 7b. The high error regions were extracted as the center of a cubic cell in (12 × 12 × 12) Å with a vacuum of 10 Å (Supplementary Note 9 and Supplementary Fig. 19), and the structures within the central sphere with a diameter of 7 Å (the dotted circle in Fig. 7c) were fixed to maintain the local atomic environment to be calibrated while the surrounding atoms (the shadow region in Fig. 7c) were fully relaxed in DFT calculation. All the labeled local configurations were added to the datasets for the finetuning of the MLP model at target space.

Fig. 7 Effective configuration identification during active learning finetuning in large scale simulations for friction target space. [Images not available. See PDF.]

a Friction evolution with time of self-mated a-C:H with 40 at.%H. b Snapshots of sliding interfacial atomic structure of self-mated a-C:H with 40 at.%H at different timepoints during the friction and the corresponding distribution of atomic force error. The purple and gray atoms represented the carbon atoms and the yellow and light blue atoms represented the hydrogen atoms in the snapshots. c Local high error region extraction and labeling during active learning finetuning. Take the snapshots of sliding interfacial atomic structure at 400 ps as an example. The high error local structures were extracted from large friction system with the high error center fixed and the surrounding atoms were fully relaxed in DFT calculation.

Strategy of uncertainty quantification

In this work, the neural networks (NN) ensemble^38,48 was utilized for the uncertainty quantification of the active learning. The MD simulations were conducted with a MLP model, and the atomic forces in the configurations during the simulations were recalculated with other two MLP models, which were trained based on the same pre-training datasets and with same neural network parameters. As the three parallel MLP models shared the same datasets and learning ability, they should have the similar descriptive accuracy for the atomic environment covered in the datasets³⁸ and the high deviation region could contain the unlabeled atomic structures in target space, which should be further finetuned. For the small-scale active learning, the deviation of model (DM) was defined as the maximum standard deviation of the predictions for the atomic forces following the formula (1):

1

$$DM=\underset{N}{\max }\sqrt[]{⟨{\left(f_i-\overline{f_i}\right)}^2⟩}$$ where i is the ith atom in the structure with N atoms, $\overline{f_i}$ is the mean atomic force value and the ensemble average ⟨…⟩ is taken over the ensemble of models³⁸. The configurations with DM above a certain threshold (0.8–2 eV/Å) were extracted for further labeling in DFT calculation. Meanwhile, an upper threshold (3–5 eV/Å) was also introduced, above which the configuration shown high uncertainty and could be unfavorable for DFT convergence and model fitting⁶⁵. The determination of the thresholds was further discussed in Supplementary Note 10.

For the large-scale active learning, the deviation of atomic force (DAF) between three parallel MLP models were calculated and the deviation of atomic force was defined with the following formula:

2

$$DAF=\sqrt{{\left({\left(f_i\right)}_{max}-\overline{f_i}\right)}^2}$$ where the ${\left(f_i\right)}_{max}$ and $\overline{f_i}$ are the maximum and mean atomic force value calculated with three parallel MLP models, respectively. The similar thresholds (lower threshold (0.8–2 eV/Å) and upper threshold (3–5 eV/Å)) were also defined for the identification of local high error regions. Notably, the DAF of all the atoms in the extracted high error regions should be lower than the upper threshold, not only for the central atoms.

Training of machine learning potentials

The MLPs was trained with neural network architecture via DeePMD-kit package⁴⁷ to learn the atomic environment information of the structures in training set, which was treated by an environment matrix and an embedding matrix to construct the descriptors and map the atomic structures with atomic forces and energies. The descriptor was set up with both angular and radial atomic configurations, and the embedding and fitting networks both included 3 hidden layers with the neural of (25, 50, 100) and (240, 240, 240) neurons, respectively. In the model, the energy was evaluated in forward propagation, while the atomic force which is the gradient of the potential energy, was calculated in the backward propagation⁶⁶. The cutoff radius for neighbor researching was set to 7 Å to consider the long-range interaction in hydrogen-carbon system. The starting and final learning rates were set to 1.0 × 10⁻³ and 1.6 × 10⁻⁸, respectively. A total of 5 × 10⁶ steps was set to train the MLP model for each iteration.

Atomic environment mapping

The similarity metric of atomic environment to generate Fig. 2: Smooth Overlap of Atomic Positions (SOAP) is a descriptor based on the idea of a local-basis expansion of the atomic neighbor density and subsequent construction of a rotationally invariant power spectrum⁵⁴. To efficiently capture the atomic neighbor environment of structures, hyper-parameters related to the radial resolution were set to: r_cut = 7 Å, n_max = 12, and the angular resolution is determined by l_max = 6. The Gaussian function width σ was proportional to the corresponding rcut by σ = r_cut/8. The DScribe python library was used to generate the SOAP descriptors of atoms⁵⁵. To visualize the similarity between high-dimensional SOAP vectors, the uniform manifold approximation and projection (UMAP) was used to reduce the dimension, and represent the similarity of SOAP vectors as the distance between points in a two-dimensional map.

Density functional theory calculations

DFT calculations were performed with Vienna ab initio simulation package (VASP)^67,68. The Perdew–Burke–Ernzerhof (PBE)⁶⁹ generalized gradient approximation (GGA) functional was used for all hydrogen-carbon systems in this work. The optB86b-vdW exchange–correlation functional⁷⁰ was utilized to approximately describe the dispersion interaction (van der Waals forces). For the static calculations, high precision (PREC = Accurate) was employed with a plane wave cutoff energy of 600 eV (ENCUT = 600) and no symmetry constraints applied (ISYM = 0). The electronic convergence criterion (EDIFF) was set at 10⁻⁶ eV, and the smallest allowed spacing between k points (KSPACING) was set at 0.14 Å⁻¹. Gaussian smearing (ISEAR = 0) was utilized with a smearing width of 0.05 eV (SIGMA = 0.05) to assist in the convergence of the calculations. For the DFT-based (“ab initio”) molecular dynamics (AIMD) simulations, the k points were set to 1*1*1 for the systems with more than 100 atoms otherwise set to 2*2*2, and the time step was 1 fs. The quenching process adopted micro canonical ensemble (NVE), while the melting and relaxing processed were conducted with regular ensemble (NVT).

MLP based molecular dynamics simulations

The MLP based molecular dynamics simulations were performed with Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)⁷¹. The time step of 0.2 fs was used with the Velocity Verlet algorithm to integrate the equations of motion.

NVT/NPT relaxation

Simulations took place in the NVT and NPT ensemble using a Nosé-Hoover chain thermostat with a temperature damping parameter of 0.02 ps. The temperature for relaxation process ranged from 0.1 to 5000 K with the interval of 50 K under 500 K and the interval of 500 K at higher range. The pressure was set from 1 to 15 GPa with the interval of 2 GPa.

Melt quenching process

The NVT ensemble was used for the melt quenching process. The carbon and hydrogen atoms were generated at a series of densities ranging from 1.6 to 3.5 gcm⁻³ and hydrogen contents ranging from 0 to 40 at.%H in the specific cell lattice and relaxed at a constant temperature of 5000 K for 3 ps, followed by a quenching process with temperature decreasing from 5000 to 300 K within 0.5 ps to obtain the amorphous carbon structures. The quenched structures were further relaxed under temperature of 300 K for 3 ps to improve their stability.

Deposition and growth of a-C:H film

Diamond (100) surface was constructed with the dimension of 28.5 * 28.5 * 15 Å to be the substrate of a-C:H film, and the carbon source for the deposition was set to C₂H₂ molecules. The bottom two atomic layers were fixed. The remaining region was coupled to Langevin thermostats at a temperature of 300 K except for a cylindrical volume, which centered at the incident direction of C₂H₂ molecules, and the atoms inside the cylinder could move freely. Periodic boundary conditions were applied along both x and y directions. The incident energy was set to 30 eV per carbon atom. The time interval between two sequential incident was set to 1 ps to achieve full relaxation of substrate and film.

Friction

The structure of a-C:H with 30 at.%H with 432 atoms was constructed for friction process with melt quenching process in AIMD simulation, and the front views of the atomistic model are presented in Supplementary Fig. 10. The models were divided into three parts along the z direction as fixed layer (considered as rigid body), Langevin thermostat layer for the temperature control and free layer where the atoms were free of constraints so that they moved according to the interatomic forces. Periodic boundary conditions were applied to both x and y directions to mimic laterally infinite surface. A constant velocity of 20 m/s was endued to the up fixed layer along the x direction to drive the sliding process.

Fracture

A constant increment of stress intensity factor $\mathrm{\Delta }{K}_{\mathrm{I}}=0.05\mathrm{MPa}\sqrt{\mathrm{m}}$ was used for fracture loading (details of the fracture loading method are introduced in Supplementary Note 11), and the results of different increments in Supplementary Fig. 20 showed that the selected $\mathrm{\Delta }{K}_{\mathrm{I}}=0.05\mathrm{MPa}\sqrt{\mathrm{m}}$ was proper. The a-C:H plate specimen for fracture simulations were generated by melt quenching with densities ranging from 2.0 to 3.2 gcm⁻³ and hydrogen contents ranging from 0 to 40 at.%H. The front view of a typical specimen is shown in Supplementary Fig. 16, with 3 nm or 10 nm radius and about 2 nm thickness, containing about 6000–10,000 or 80,000–100,000 atoms. The outmost layer with thickness of about 0.5 nm was set rigid during the simulation, and the atoms within R_free were fully relaxed at 0.1 K with NVT ensemble after each load step.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant Nos. 52225502, 52305200), New Cornerstone Science Foundation through the XPLORER PRIZE, National Key Research and Development Program of China (2024YFB3410201, 2018YFB0704300), Scientific and Technological Project of Yunnan Precious Metals Laboratory (YPML-20240502088), Key Research and Development Program of Yunnan Province (202203ZA080001).

Author contributions

W.C. and Z.X. contributed equally to this work. W.C. and T.M. conceived and designed the project. W.C. conducted the construction of the transferable active machine learning scheme, training of MLP models, MD simulations, and DFT calculations. Z.X. conducted the dimensionality reduction of atomic environments in datasets, finetuning of MLP model in fracture target space, and test of empirical interatomic potentials. W.C., Z.X., A.S., and T.M. contributed to the writing, review and editing of the paper. K.W. and L.G. contributed to the discussion and results analysis. W.C., Z.X., and T.M. contributed to writing the manuscript.

Data availability

The constructed machine learning potential in this work has been shared online: https://github.com/TianbaoMa/MLP_CH_deepmd. The datasets used and/or analyzed during the current study available from the corresponding author on reasonable request.

Ethics declarations

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41524-025-01629-y.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Li, X et al. Large-area synthesis of high-quality and uniform graphene films on copper foils. Science; 2009; 324, pp. 1312-1314.1:CAS:528:DC%2BD1MXms12gtbY%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19423775][DOI: https://dx.doi.org/10.1126/science.1171245]

2. Luo, K et al. Coherent interfaces govern direct transformation from graphite to diamond. Nature; 2022; 607, pp. 486-491.1:CAS:528:DC%2BB38XhslentrnP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35794481][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9300464][DOI: https://dx.doi.org/10.1038/s41586-022-04863-2]

3. Utsumi, S. et al. Giant nanomechanical energy storage capacity in twisted single-walled carbon nanotube ropes. Nat. Nanotechnol. https://doi.org/10.1038/s41565-024-01645-x (2024).

4. Bhaskaran, H et al. Ultralow nanoscale wear through atom-by-atom attrition in silicon-containing diamond-like carbon. Nat. Nanotechnol.; 2010; 5, pp. 181-185.1:CAS:528:DC%2BC3cXislaisLw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20118919][DOI: https://dx.doi.org/10.1038/nnano.2010.3]

5. Biswas, C; Lee, Y. Graphene versus carbon nanotubes in electronic devices. Adv. Funct. Mater.; 2011; 21, pp. 3806-3826.1:CAS:528:DC%2BC3MXhtFKhs7fO [DOI: https://dx.doi.org/10.1002/adfm.201101241]

6. Raccichini, R; Varzi, A; Passerini, S; Scrosati, B. The role of graphene for electrochemical energy storage. Nat. Mater.; 2015; 14, pp. 271-279.1:CAS:528:DC%2BC2cXitFKltLbP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25532074][DOI: https://dx.doi.org/10.1038/nmat4170]

7. Yu, Q et al. Ion energy-induced nanoclustering structure in a-C:H film for achieving robust superlubricity in vacuum. Friction; 2022; 10, pp. 1967-1984.1:CAS:528:DC%2BB38XjtV2gtbo%3D [DOI: https://dx.doi.org/10.1007/s40544-021-0555-8]

8. Wang, K et al. Unraveling the friction evolution mechanism of diamond‐like carbon film during nanoscale running‐in process toward superlubricity. Small; 2021; 17, 1:CAS:528:DC%2BB3cXisFGmu7%2FP [DOI: https://dx.doi.org/10.1002/smll.202005607] 2005607.

9. Cahangirov, S; Topsakal, M; Akturk, E; Sahin, H; Ciraci, S. Two- and one-dimensional honeycomb structures of silicon and germanium. Phys. Rev. Lett.; 2009; 102, 236804.1:STN:280:DC%2BD1MrjvFCktg%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19658958][DOI: https://dx.doi.org/10.1103/PhysRevLett.102.236804]

10. Erdemir, A. Genesis of superlow friction and wear in diamondlike carbon films. Tribol. Int.; 2004; 37, pp. 1005-1012.1:CAS:528:DC%2BD2cXptlSnt70%3D [DOI: https://dx.doi.org/10.1016/j.triboint.2004.07.018]

11. Shin, D; Kim, S. Effect of hydrogen embrittlement on mechanical characteristics of DLC-coating for hydrogen valves of FCEVs. npj Mater. Degrad.; 2024; 8, 1:CAS:528:DC%2BB2cXhtlShtbvI [DOI: https://dx.doi.org/10.1038/s41529-024-00460-y] 47.

12. Pastewka, L; Moser, S; Gumbsch, P; Moseler, M. Anisotropic mechanical amorphization drives wear in diamond. Nat. Mater.; 2011; 10, pp. 34-38.1:CAS:528:DC%2BC3cXhsFKisrvE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21113152][DOI: https://dx.doi.org/10.1038/nmat2902]

13. Chen, W; Ma, T. Interfacial tribochemical kinetics: a new perspective on superlubricity of diamond-like carbon films. Sci. China Technol. Sci.; 2024; 67, pp. 2050-2052.1:CAS:528:DC%2BB2cXhsFeju7%2FK [DOI: https://dx.doi.org/10.1007/s11431-024-2720-x]

14. Carpenter, C; Ramasubramaniam, A; Maroudas, D. Analysis of vacancy-induced amorphization of single-layer graphene. Appl. Phys. Lett.; 2012; 100, 203105. [DOI: https://dx.doi.org/10.1063/1.4717116]

15. Qiu, K et al. Self-healing of fractured diamond. Nat. Mater.; 2023; 22, pp. 1317-1323.1:CAS:528:DC%2BB3sXhvFyrtb3N [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37735525][DOI: https://dx.doi.org/10.1038/s41563-023-01656-4]

16. Kunze, T et al. Wear, plasticity, and rehybridization in tetrahedral amorphous carbon. Tribol. Lett.; 2014; 53, pp. 119-126.1:CAS:528:DC%2BC3sXhslWjur%2FO [DOI: https://dx.doi.org/10.1007/s11249-013-0250-7]

17. Chu, P; Li, L. Characterization of amorphous and nanocrystalline carbon films. Mater. Chem. Phys.; 2006; 96, pp. 253-277.1:CAS:528:DC%2BD28XhtV2ntLo%3D [DOI: https://dx.doi.org/10.1016/j.matchemphys.2005.07.048]

18. Zhang, S et al. Domino-like stacking order switching in twisted monolayer–multilayer graphene. Nat. Mater.; 2022; 21, pp. 621-626.1:CAS:528:DC%2BB38XhtVKqtLnJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35449221][DOI: https://dx.doi.org/10.1038/s41563-022-01232-2]

19. Dag, S; Ciraci, S. Atomic scale study of superlow friction between hydrogenated diamond surfaces. Phys. Rev. B; 2004; 70, [DOI: https://dx.doi.org/10.1103/PhysRevB.70.241401] 241401.

20. Pei, L et al. Regulating vacuum tribological behavior of a-C:H film by interfacial activity. J. Phys. Chem. Lett.; 2021; 12, pp. 10333-10338.1:CAS:528:DC%2BB3MXit1yhu77O [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34664963][DOI: https://dx.doi.org/10.1021/acs.jpclett.1c02841]

21. Chacham, H; Kleinman, L. Instabilities in diamond under high shear stress. Phys. Rev. Lett.; 2000; 85, 4904.1:CAS:528:DC%2BD3cXosFals7Y%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11102147][DOI: https://dx.doi.org/10.1103/PhysRevLett.85.4904]

22. Bartók, A; Kermode, J; Bernstein, N; Csányi, G. Machine learning a general-purpose interatomic potential for silicon. Phys. Rev. X; 2018; 8, 041048.

23. Pastewka, L; Klemenz, A; Gumbsch, P; Moseler, M. Screened empirical bond-order potentials for Si-C. Phys. Rev. B; 2013; 87, [DOI: https://dx.doi.org/10.1103/PhysRevB.87.205410] 205410.

24. Pastewka, L; Pou, P; Pérez, R; Gumbsch, P; Moseler, M. Describing bond-breaking processes by reactive potentials: Importance of an environment-dependent interaction range. Phys. Rev. B; 2008; 78, [DOI: https://dx.doi.org/10.1103/PhysRevB.78.161402] 161402.

25. Lin, S et al. Machine-learning potentials for nanoscale simulations of tensile deformation and fracture in ceramics. npj Comput Mater.; 2024; 10, 67.1:CAS:528:DC%2BB2cXotVWgsr0%3D [DOI: https://dx.doi.org/10.1038/s41524-024-01252-3]

26. Yu, M; Zhao, Z; Guo, W; Zhang, Z. Fracture toughness of two-dimensional materials dominated by edge energy anisotropy. J. Mech. Phys. Solids; 2024; 186, 105579.1:CAS:528:DC%2BB2cXkvFaitLY%3D [DOI: https://dx.doi.org/10.1016/j.jmps.2024.105579]

27. Zhang, D; Yi, P; Lai, X; Peng, L; Li, H. Active machine learning model for the dynamic simulation and growth mechanisms of carbon on metal surface. Nat. Commun.; 2024; 15, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38184678][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10771457][DOI: https://dx.doi.org/10.1038/s41467-023-44525-z] 344.

28. Hedman, D et al. Dynamics of growing carbon nanotube interfaces probed by machine learning-enabled molecular simulations. Nat. Commun.; 2024; 15, 1:CAS:528:DC%2BB2cXht1ClsrnF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38744824][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11094095][DOI: https://dx.doi.org/10.1038/s41467-024-47999-7] 4076.

29. Marchant, G; Caro, M; Karasulu, B; Pártay, L. Exploring the configuration space of elemental carbon with empirical and machine learned interatomic potentials. npj Comput. Mater.; 2023; 9, 131. [DOI: https://dx.doi.org/10.1038/s41524-023-01081-w]

30. Caro, M; Deringer, V; Koskinen, J; Laurila, T; Csányi, G. Growth mechanism and origin of high sp³ content in tetrahedral amorphous carbon. Phys. Rev. Lett.; 2018; 120, 166101.1:CAS:528:DC%2BC1MXltFSgt78%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29756912][DOI: https://dx.doi.org/10.1103/PhysRevLett.120.166101]

31. Deringer, V; Csányi, G. Machine learning based interatomic potential for amorphous carbon. Phys. Rev. B; 2017; 95, [DOI: https://dx.doi.org/10.1103/PhysRevB.95.094203] 094203.

32. Rowe, P; Deringer, V; Gasparotto, P; Csányi, G; Michaelides, A. An accurate and transferable machine learning potential for carbon. J. Chem. Phys.; 2020; 153, 034702.1:CAS:528:DC%2BB3cXhsVWrtr7M [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32716159][DOI: https://dx.doi.org/10.1063/5.0005084]

33. Faraji, S; Liu, M. Transferable machine learning interatomic potential for carbon hydrogen systems. Phys. Chem. Chem. Phys.; 2024; 26, 22346.1:CAS:528:DC%2BB2cXhvVOnurrL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39140158][DOI: https://dx.doi.org/10.1039/D4CP02300E]

34. Ibragimova, R., Kuklin, M., Zarrouk, T. & Caro, M. Unifying the description of hydrocarbons and hydrogenated carbon materials with a chemically reactive machine learning interatomic potential. Chem. Phys. 08194; arXiv:2409.08194 (2024).

35. Zhang, S et al. Exploring the frontiers of condensed-phase chemistry with a general reactive machine learning potential. Nat. Chem.; 2024; 16, pp. 727-734.1:CAS:528:DC%2BB2cXltlyiu74%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38454071][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11087274][DOI: https://dx.doi.org/10.1038/s41557-023-01427-3]

36. Qi, J; Ko, T; Wood, B; Pham, T; Ong, S. Robust training of machine learning interatomic potentials with dimensionality reduction and stratified sampling. npj Comput. Mater.; 2024; 10, 43. [DOI: https://dx.doi.org/10.1038/s41524-024-01227-4]

37. Hayashi, K et al. Tribochemical reaction dynamics simulation of hydrogen on a diamond-like carbon surface based on tight-binding quantum chemical molecular dynamics. J. Phys. Chem. C.; 2011; 115, pp. 22981-22986.1:CAS:528:DC%2BC3MXht1yqu7fK [DOI: https://dx.doi.org/10.1021/jp207065n]

38. Zhang, L; Lin, D; Wang, H; Car, R; E, W. Active learning of uniformly accurate interatomic potentials for materials simulation. Phys. Rev. Mater.; 2019; 3, 023804.1:CAS:528:DC%2BC1MXhtFGrsrnF [DOI: https://dx.doi.org/10.1103/PhysRevMaterials.3.023804]

39. Thaler, S; Doehner, G; Zavadlav, J. Scalable Bayesian uncertainty quantification for neural network potentials: promise and pitfalls. J. Chem. Theory Comput.; 2023; 19, 4520.1:CAS:528:DC%2BB3sXmslKmsro%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37014758][DOI: https://dx.doi.org/10.1021/acs.jctc.2c01267]

40. Smith, J et al. Automated discovery of a robust interatomic potential for aluminum. Nat. Commun.; 2021; 12, 1:CAS:528:DC%2BB3MXltFOmtrs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33623036][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7902823][DOI: https://dx.doi.org/10.1038/s41467-021-21376-0] 1257.

41. Tan, A; Urata, S; Goldman, S; Dietschreit, J; Gómez-Bombarelli, R. Single-model uncertainty quantification in neural network potentials does not consistently outperform model ensembles. npj Comput. Mater.; 2023; 9, 225.1:CAS:528:DC%2BB3sXis1elsbrK [DOI: https://dx.doi.org/10.1038/s41524-023-01180-8]

42. Vandermause, J et al. On-the-fly active learning of interpretable Bayesian force fields for atomistic rare events. npj Comput. Mater.; 2020; 6, 20. [DOI: https://dx.doi.org/10.1038/s41524-020-0283-z]

43. Gal, Y; Ghahramani, Z. Dropout as a Bayesian approximation: representing model uncertainty in deep learning. Proc. Mach. Learn.; 2016; 48, 1050.

44. Podryabinkin, E; Shapeev, A. Active learning of linearly parametrized interatomic potentials. Comput. Mater. Sci.; 2017; 140, 171.1:CAS:528:DC%2BC2sXhsVOiurrJ [DOI: https://dx.doi.org/10.1016/j.commatsci.2017.08.031]

45. Molina-Taborda, A; Cossio, P; Lopez-Acevedo, O; Gabrié, M. Active learning of Boltzmann samplers and potential energies with quantum mechanical accuracy. J. Chem. Theory Comput.; 2024; 20, 8833.1:CAS:528:DC%2BB2cXitFehur3M [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39370622][DOI: https://dx.doi.org/10.1021/acs.jctc.4c00506]

46. Jung, G; Choi, J; Lee, S. Active learning of neural network potentials for rare events. Digit. Discov.; 2024; 3, 514. [DOI: https://dx.doi.org/10.1039/D3DD00216K]

47. Wang, H; Zhang, L; Han, J; E, W. DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics. Comput. Phys. Commun.; 2018; 228, pp. 178-184.1:CAS:528:DC%2BC1cXmt1ajsrk%3D [DOI: https://dx.doi.org/10.1016/j.cpc.2018.03.016]

48. Jinnouchi, R; Miwa, K; Karsai, F; Kresse, G; Asahi, R. On-the-fly active learning of interatomic potentials for large-scale atomistic simulations. J. Phys. Chem. Lett.; 2020; 11, pp. 6946-6955.1:CAS:528:DC%2BB3cXhsFWgsLfK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32787192][DOI: https://dx.doi.org/10.1021/acs.jpclett.0c01061]

49. Zhang, L; Wang, H; E, W. Reinforced dynamics for enhanced sampling in large atomic and molecular systems. J. Chem. Phys.; 2018; 148, 124113. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29604808][DOI: https://dx.doi.org/10.1063/1.5019675]

50. Chen, W. et al. Ultralow-friction at cryogenic temperature induced by hydrogen correlated quantum effect. Small20, 2400083 (2024).

51. Marks, N; McKenzie, D; Pailthorpe, B; Bernasconi, M; Parrinello, M. Ab initio simulations of tetrahedral amorphous carbon. Phys. Rev. B; 1996; 54, 9703.1:CAS:528:DyaK28Xmt1yku7o%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.54.9703]

52. McCulloch, D; McKenzie, D; Goringe, C. Ab initio simulations of the structure of amorphous carbon. Phys. Rev. B; 2000; 61, 2349.1:CAS:528:DC%2BD3cXlvVGltA%3D%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.61.2349]

53. Salinas Ruiz, V et al. Interplay of mechanics and chemistry governs wear of diamond-like carbon coatings interacting with ZDDP-additivated lubricants. Nat. Commun.; 2021; 12, 1:CAS:528:DC%2BB3MXhslylu77F [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34315887][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8316475][DOI: https://dx.doi.org/10.1038/s41467-021-24766-6] 4550.

54. Bartók, A; Kondor, R; Csányi, G. On representing chemical environments. Phys. Rev. B; 2013; 87, 184115. [DOI: https://dx.doi.org/10.1103/PhysRevB.87.184115]

55. Himanen, L et al. DScribe: Library of descriptors for machine learning in materials science. Comput. Phys. Commun.; 2020; 247, 106949.1:CAS:528:DC%2BC1MXhvV2itrzI [DOI: https://dx.doi.org/10.1016/j.cpc.2019.106949]

56. McInnes, L., Healy, J. & Melville, J. Umap: uniform manifold approximation and projection for dimension reduction. Preprint at https://arXiv.org/ abs/1802.03426 (2018).

57. Wang, H; Guo, X; Zhang, L; Wang, H; Xue, J. Deep learning inter-atomic potential model for accurate irradiation damage simulations. Appl. Phys. Lett.; 2019; 114, 244101. [DOI: https://dx.doi.org/10.1063/1.5098061]

58. Wolloch, M; Levita, G; Restuccia, P; Righi, MC. Interfacial charge density and its connection to adhesion and frictional forces. Phys. Rev. Lett.; 2018; 121, 026804.1:CAS:528:DC%2BC1MXltFSgt74%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30085711][DOI: https://dx.doi.org/10.1103/PhysRevLett.121.026804]

59. Qamar, M; Mrovec, M; Lysogorskiy, Y; Bochkarev, A; Drautz, R. Atomic cluster expansion for quantum-accurate large-scale simulations of carbon. J. Chem. Theory Comput.; 2023; 19, pp. 5151-5167.1:CAS:528:DC%2BB3sXht1Kgt7fI [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37347981][DOI: https://dx.doi.org/10.1021/acs.jctc.2c01149]

60. Khosrownejad, S; Kermodw, J; Pastewka, L. Quantitative prediction of the fracture toughness of amorphous carbon. Phys. Rev. Mater.; 2021; 5, 023602.1:CAS:528:DC%2BB3MXmslWmt7w%3D [DOI: https://dx.doi.org/10.1103/PhysRevMaterials.5.023602]

61. Jonnalagadda, K; Cho, S; Chasiotis, I; Friedmann, T; Sullivan, J. Effect of intrinsic stress gradient on the effective mode-I fracture toughness of amorphous diamond-like carbon films for MEMS. J. Mech. Phys. Solids; 2008; 56, pp. 388-401.1:CAS:528:DC%2BD1cXhvFyhtbs%3D [DOI: https://dx.doi.org/10.1016/j.jmps.2007.05.013]

62. Schaufler, J; Schmid, C; Durst, K; Göken, M. Determination of the interfacial strength and fracture toughness of a-C:H coatings by in-situ microcantilever bending. Thin Solid Films; 2012; 522, pp. 480-484.1:CAS:528:DC%2BC38XhsVamsLzP [DOI: https://dx.doi.org/10.1016/j.tsf.2012.08.031]

63. Shin, D; Jang, D. Crack-tip plasticity and intrinsic toughening in nano-sized brittle amorphous carbon. Int. J. Plast.; 2020; 127, 102642.1:CAS:528:DC%2BB3cXivVSrsQ%3D%3D [DOI: https://dx.doi.org/10.1016/j.ijplas.2019.102642]

64. Tamura, M; Kumagai, T. Hydrogen permeability of diamondlike amorphous carbons. J. Vac. Sci. Technol. A; 2017; 35, [DOI: https://dx.doi.org/10.1116/1.4977106] 04D101.

65. Erhard, L; Rohrer, J; Albe, K; Deringer, V. Modelling atomic and nanoscale structure in the silicon–oxygen system through active machine learning. Nat. Commun.; 2024; 15, 1:CAS:528:DC%2BB2cXlsVKms7w%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38431626][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10908788][DOI: https://dx.doi.org/10.1038/s41467-024-45840-9] 1927.

66. Lu, D et al. 86 PFLOPS deep potential molecular dynamics simulation of 100 million atoms with ab initio accuracy. Comput. Mater. Sci.; 2021; 259, 107624.1:CAS:528:DC%2BB3cXitVKgt7vI

67. Kresse, G; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci.; 1996; 6, pp. 15-50.1:CAS:528:DyaK28XmtFWgsrk%3D [DOI: https://dx.doi.org/10.1016/0927-0256(96)00008-0]

68. Kresse, G; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B; 1996; 54, pp. 11169-11186.1:CAS:528:DyaK28Xms1Whu7Y%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.54.11169]

69. Perdew, J; Burke, K; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett.; 1996; 77, pp. 3865-3868.1:CAS:528:DyaK28XmsVCgsbs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10062328][DOI: https://dx.doi.org/10.1103/PhysRevLett.77.3865]

70. Dion, M; Rydberg, H; Schröder, E; Langreth, D; Lundqvist, B. Van der Waals density functional for general geometries. Phys. Rev. Lett.; 2004; 92, 246401.1:STN:280:DC%2BD2czksVKnsA%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15245113][DOI: https://dx.doi.org/10.1103/PhysRevLett.92.246401]

71. Thompson, A et al. LAMMPS—a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput. Phys. Commun.; 2021; 271, 108171. [DOI: https://dx.doi.org/10.1016/j.cpc.2021.108171]