Atomistic fracture in bcc iron revealed by active learning of Gaussian approximation potential

The prediction of atomistic fracture mechanisms in body-centred cubic (bcc) iron is essential for understanding its semi-brittle nature. Existing atomistic simulations of the crack-tip deformation mechanisms under mode-I loading based on classical interatomic potentials yield contradicting predictions. To enable fracture prediction with quantum accuracy, we develop a Gaussian approximation potential (GAP) using an active learning strategy by extending a density functional theory (DFT) database of ferromagnetic bcc iron. We apply the active learning algorithm and obtain a Fe GAP model with a maximum predicted error of 8 meV/atom over a broad range of stress intensity factors (SIFs) and for four crack systems. The learning efficiency of the approach is analysed, and the predicted critical SIFs are compared with Griffith and Rice theories. The simulations reveal that cleavage along the original crack plane is the crack tip mechanism for {100} and {110} crack planes at T=0K, thus settling a long-standing dispute. Our work also highlights the need for a multiscale approach to predicting fracture and intrinsic ductility, whereby finite temperature, finite loading rate effects and pre-existing defects (e.g. nanovoids, dislocations) should be taken explicitly into account.


Introduction
Brittle fracture is a key failure mechanism of body-centred cubic (bcc) transition metals, which limits their application and can jeopardise the safety of infrastructures.Brittle fracture usually takes place in the form of cleavage, which can be accompanied by dislocation plasticity.The competition between thermally activated dislocation mobility and (atomic-scale) crack tip deformation mechanisms controls the fracture process of bcc iron [1][2][3][4] .Experiments on single-crystal bcc iron reveal that cleavage takes place on {100} planes for pre-existing {100}, and {110} crack planes with 100 and 110 crack fronts 5 .
However, atomic-scale experimental data is unavailable for crack tips in bcc iron, leaving the atomistic crack-tip mechanisms unclear.
Atomistic modelling based on molecular dynamics (MD) approach has been widely used to predict the crack-tip mechanisms [6][7][8][9][10][11][12][13][14][15] (see Supplementary Material S1 for a summary of existing studies).Mode-I atomistic crack tip simulations at T=0K arXiv:2208.05912v2[cond-mat.mtrl-sci]14 Sep 2022 are typically used to assess the intrinsic ductility of metals, which is controlled by the competition between crack propagation (cleavage) and dislocation emission from the crack tip 3,4,16 .Yet, as highlighted in extensive reviews 14,15 , classical interatomic potentials (IAPs) for bcc iron predict contradicting the crack-tip mechanisms at T=0K (i.e.cleavage, crack propagation planes, dislocation emission, and phase transition) for the same crack system.Remarkably, none of them agree with the low-temperature experimental fracture mechanisms as a function of the pre-existing crack system 5 .A indicates amorphous structure forming at the crack tip.The crack propagates on the original crack plane if no symbol is specified.K G and K Ie are predictions according to Griffith 17 and Rice theories 2 for cleavage and dislocation emission, respectively.Fig. 1 shows the critical stress intensity factor (K Ic ) under mode-I loading as computed by performing molecular statics (MS) simulations using seven classical IAPs (see MS/MD simulation setup in Methods) that were not investigated in Ref. 14 nor Ref. 15 , including two embedded atom method (EAM) potentials 18,19 , three Modified EAM (MEAM) potentials ?, 21,22 , one bond order potential (BOP) 23 , and one angular dependent potential (ADP) 24 .Results are compared with the critical stress intensity factors according to Griffith theory 17 where γ s is the surface energy and B is a constant determined by the elastic constants (see Supplementary Material S2.2).
Dislocation emission is predicted according to Rice theory 2 with an anisotropic implementation ?
where o(θ , φ ) is a function of θ and φ .θ is the angle between slip plane and crack plane, φ is the angle between the slip direction and a vector lying on slip plane and perpendicular to the crack-front direction, and γ usf is the unstable stacking fault energy of the slip plane.The computation details of o(θ , φ ) and F 12 (θ ) are given in Supplementary Material S2.3.These potential-dependent properties (elastic constants, γ s and γ usf ) used for computing K G and K Ie are listed in Supplementary Material S3.For (100)[010] (crack plane/crack front) crack system, no IAP predicts pure cleavage on the original crack plane, in stark contrast with both the low-temperature experimental observations 5 and the theoretical predictions (K G and K Ie ).Instead, all IAPs predict either cleavage on {110} planes (rather than {100} planes), or phase transition, or cleavage accompanied by phase transition.With only two exceptions (MEAM- The lack of agreement between atomistic predictions and experiments of fracture modes and critical stress intensity factors motivates the development of a potential that is capable of reproducing the experimental observations.Machine learning (ML) potentials enable MD/MS simulation to describe a complex potential energy surface (PES) with density functional theory (DFT) accuracy, and are order of magnitudes faster than DFT.Among the existing ML potential frameworks, Gaussian Approximation Potential (GAP) has been shown to accurately describe the complex motion of screw dislocations in bcc iron 25 and the cleavage process in silicon 26 .As a Gaussian process regression method, GAP can predict both the mean value and the variance of the atomic energies 27 .The square root of the GAP predicted variance is used as an indicator for error, which is useful to evaluate the quality of the model and practical in iterative/active learning 28,29 .
We first consider an existing GAP for ferromagnetic bcc iron 30 , hereafter named Fe-GAP18.The Fe-GAP18 database has been developed to study thermodynamics and defects in bcc iron 31 .The database includes stretched primitive cells, point defects, interfaces etc.Using the GAP predicted variance, we show (Fig. 2a) that Fe-GAP18 cannot predict the fracture behaviour of bcc iron.In particular, Fig. 2a reports the maximum per-atom energy error of crack-tip atoms during mode-I MS simulations at multiple applied K I s.Only crack propagation along (110) [1 10] is predicted within an accuracy of ∼10 meV/atom.Fig. 2b shows snapshots with typical artifacts produced during the fracture process.We argue that Fe-GAP18 is unable to capture the fracture behaviour due to the lack of crack tip atomic environments in the DFT training dataset rather than due to an intrinsic limitation of GAP.Therefore, here we develop a systematic strategy to improve GAP for fracture predictions, including a preliminary fracture-relevant database extension, which is subsequently supplemented with an active learning algorithm.We show that the new database enables GAP to predict atomistic fracture mechanisms within an accuracy of 8 meV/atom at T=0K and on the order of ∼10 meV/atom at finite temperatures.The paper is organised as follows.We describe the approach and the active learning algorithm used to expand the DFT database in section Results "Preparation of the Reference Database".In section Results "Implementation of active learning", we implement the algorithm and show the convergence of our approach (hereafter indicated as Fe-GAP22).Based on Fe-GAP22, we predict the fracture behaviour of single crystal bcc iron at T=0K and at finite temperatures under high loading rate conditions in section Results "Prediction of fracture mechanisms and critical stress intensity factors (K Ic )".In section Discussion, we summarise our findings and outline further research directions to be explored with the novel, quantum-accurate Fe-GAP22.

Preparation of the Reference Database
Preliminary Fracture-relevant Database Enrichment We first extend the existing DFT database of Fe-GAP18 31 with primitive cells strained according to the large deformations that may occur at crack tips.The strain states near a crack tip under plane strain and plane stress conditions are estimated from classical linear elastic crack-tip fields.Thus, we find that the maximum tensile strain is 0.27 at the point that is 0.1 nm away from the crack-tip under K I =1.5 MPa √ m.By analysing the strain components of the primitive cell from the original database 31 comprising 6, 000 distorted primitive bcc cells, we found that 95% of original data is within the strain of 0.15 and the rest are highly stretched along one crystal direction (2 ∼ 4 times).We therefore create a new set of highly strained primitive cells to enrich the DFT dataset to encompass the classical linear elastic predicted strain states, referred to as DB9.Note that the original database from Dragoni 31 includes 8 subsets, referred to as DB1-DB8.In order to ensure an unbiased sampling of atomic environments, we use the random uniform sampling algorithm SOBOL 33  shows that DB9 expands the original database and encompasses a more extensive set of strain states than those calculated by LEFM (see Supplementary Material S5).Common neighbour analysis shows that the original database and DB9 also include face-cubic centred (fcc) structures.Besides, we include distorted hexagonal close-packed (hcp) primitive cells (referred to as DB10) to ensure that GAP is capable of predicting the hcp phase.
The magnetic contribution to the total energy in iron plays an important role and thus cannot be neglected.Our DFT calculation of the primitive bcc cell shows that the magnetisation degree decreases when the volume decreases.The primitive cell becomes completely non-magnetic when the volume reduction is larger than 37.5%.We found that including magnetic states beyond ferromagnetic deteriorates the performance of GAP, leading to unrealistic predictions, such as negative surface energy.This is because the transition between ferromagnetic and other magnetic states (e.g.non-magnetic) is not unique (i.e., not smooth and not single-valued), which gives rise to a partly spurious GAP predicted PES since magnetism is not explicitly accounted for in the current GAP implementation.Depending on the magnetisation parameter settings in DFT calculations, the self-consistent calculation may be stuck in a ground state of ferromagnetism that has larger energy than the non-magnetic ground state, or vice-versa.The mixture of configurations with the same atomic environment and distinct magnetic states leads to a spurious GAP fitted PES.Here, we make the assumption that the crack tip stays ferromagnetic, since the locally highly strained crack-tip bonds are surrounded by ferromagnetic iron atoms; for this reason, we train Fe-GAP22 on a consistent database of ferromagnetic configurations.Furthermore, DFT calculations used within the active learning scheme (see next Section) show that the crack-tip atoms stay ferromagnetic (see Supplementary Material S6).
Bond rupture and dislocation emission from the crack tip are two fundamental crack-tip mechanisms 34 .The DFT database should be sufficiently informed with configurations relevant to those mechanisms.Maresca et al. showed that Fe-GAP18 predicts the single-hump Peierls potential and the compact core structure of screw dislocations with DFT accuracy 25 .However, Fe-GAP18 does not include the separation process of atomic planes explicitly.Here, we add rigid separation configurations of {100}, {110}, {111} and {112} crystallographic planes for thin and thick slabs (referred to as DB11).The rigid separation includes two sets, i.e., ideal crystal structures and "rattled" structures.The rattled structures are obtained by adding a Gaussian noise N ∼ [0, 0.05a 0 ] (a 0 lattice constant) to the three Cartesian coordinates of all atoms to explore the PES beyond the highly symmetric, minimum energy configurations.In either case, separations are performed from 0 to 4 Å with a step of 0.4 Å.
The detailed information about the preliminary database is listed in Supplementary Material S7.As shown in Fig. 4, from the original database ("Original") to preliminary database ("Iter-0"), adding those relevant configurations considerably reduces the maximum GAP predicted per-atom energy error during the fracture simulations.Yet, GAP based on "Iter-0" database still cannot predict the fracture process within a target predicted energy error (see Supplementary Material S8).

Active Learning Database
Since the crack tip deformation fields are inhomogeneous, while DB9 and DB10 include only homogenous deformations, the "Iter-0" GAP is not well-informed about crack tip local atomic environments.By virtue of the localisation assumption that the atom only interacts with its neighbouring atoms within a finite cut-off radius (5 Å in this study), the crack tip atomic environment can be represented by a computable DFT cell that contains the atomistic crack tip configuration.To learn efficiently from the crack tip, we develop an active learning algorithm by extracting crack tips from extrapolating configurations.The algorithm (see Data Availability) automatically identifies the extrapolated atomic environments with respect to a predefined criterion and constructs a periodic crack cell that is computable with plane-wave DFT while minimising spurious surface effects.

5/17
The GAP predicted variance per-atom is a measure of the extrapolating degree from the existing DFT dataset, which is used to assess the accuracy of the new configurations during the simulation.The square root of GAP predicted variance per-atom is compared with a predefined criterion to determine the extrapolation.The predefined criterion is set to 10 meV/atom in this study.Table 1.Summary of active learning algorithm.Roman numbers correspond to the ones shown in Fig. 3a.
Step Function Code

II
Evaluation of the GAP predicted variance QUIP 36

Construction of the DFT cell Python
V DFT calculation Quantum Espresso 37

VI
Fitting GAP model QUIP Fig. 3a shows the workflow, schematically indicating the main steps of the active learning algorithm.We first run a K-test and evaluate the GAP predicted variance for all frames dumped during the fracture simulation.Then, the GAP predicted per-atom energy error is compared with the predefined extrapolation criterion.A periodic cell is subsequently constructed by extracting the atomic configuration from the extrapolating frame.DFT calculations of the periodic crack-tip cell are performed, and the new DFT data are added to the training database.Next, GAP is retrained with the new DFT database.We repeat this process until convergence is achieved, i.e., the maximum GAP predicted per-atom energy error is less than the predefined criterion during the entire fracture process.A summary of the software/code used in the active learning scheme is listed in Table 1.
One of the main challenges is the construction of an appropriate DFT cell that contains the crack tips in step (IV) of Fig. 3a.
It is necessary to prevent spurious boundaries to avoid learning irrelevant/artificial atomic environments.The active learning has been applied to study the screw dislocation in bcc tungsten, in which the effective extrapolative configuration is constructed by taking advantage of symmetry 38 .The idea here is to symmetrise the crack tip along the crack plane normal to construct an approximate periodic cell.The construction of the DFT crack configuration is illustrated in Fig. 3b.First we identify the atom (ID 0 ) with the largest GAP predicted error at the crack tip in the deformed configuration.Then, a group of atoms (ID g ) in a square centred on ID 0 is selected in the undeformed configuration; once identified, these atoms form a crack tip configuration in the deformed state.Another crack tip configuration is created by duplicating the original one and rotating it by 180 • .We construct the periodic DFT configuration by aligning the rotated replica and the original one by a rigid displacement.
The periodic DFT configuration is surface-free and preserves the local atomic environment of the atom with the largest GAP predicted error.Such construction allows both efficient DFT calculations and learning speed, as illustrated in the following section.

Implementation of active learning
We implemented the active learning algorithm with the original SOAP 39 and an optimised SOAP descriptor 32 (referred to as "Turbo SOAP").The presentation here focuses on active learning results with Turbo SOAP because this leads to higher computing speed and faster convergence.The results based on the original SOAP are discussed in Supplementary Material S9.We applied the algorithm to two crack systems in parallel, i.e., (100)[011] and (110)[001], such that GAP can learn the local environments from both {100} and {110} crack planes at the same time.Two sets of small crack tip regions are used, i.e., configurations extracted from fracture simulations with and without Gaussian noise (N ∼ [0, 0.05a 0 ]) added to the three Cartesian coordinates.Fig. 4 shows the convergence of the maximum GAP predicted per-atom energy error at multiple K I 's with respect to active learning iterations.The GAP predicted energy error is plotted as a function of K I during a mode-I K-test at T=0K.
As shown in Fig. 4a and b, the maximum GAP predicted error converges to 8 meV/atom for crack systems (100)[011] and (110)[001] after two iterations of active learning.The GAP predicted error of (100)[010] converges to 8 meV/atom at the third iteration (Fig. 4c) while (100)[011] and (110)[001] crack systems remain of the same accuracy.As shown in Fig. 4d, the GAP predicted error of (110)[1 10] increases after two iterations (Iter-0 and Iter-1) and then decreases but does not converge to 8 meV/atom after another two iterations (Iter-1 and Iter-2).This is because the original GAP is able to extrapolate the atomic environments that are close to crack tip of (110)[ 1 10].However, the extrapolation is unstable and is easily deteriorated by adding more data, which means that the prediction is not converged yet.We further added a crack tip of (110)[1 10] for one additional iteration of training and converged the GAP predicted error to 8 meV/atom (Fig. 4d).

Prediction of fracture mechanisms and critical stress intensity factor (K Ic )
We predicted the lattice constant, elastic constants and surface energies as a benchmark of the newly developed Fe-GAP22 and found similar accuracy as Fe-GAP18 (see Supplementary Material S10).Subsequently, we performed fracture simulation of single-crystal ferromagnetic iron at T=0K and finite temperatures (T=1, 10, 100, 200, 300K) under high loading rate (10 9 MPa √ ms −1 ).The fracture mechanism at T=0K is found to be cleavage on the original plane in all considered cases, and the GAP predicted per-atom energy error remains below 8 meV/atom during the entire cleavage process (Fig. 5).Neither phase transformation nor {110} planar faults appear according to Fe-GAP22 predictions.
The cleavage is always taking place on the plane that has the maximum normal stress.To evaluate the normal stress, we performed rigid body separation of {100} and {110} planes by using DFT and Fe-GAP22 to calculate the universal binding energy relation (UBER) (Fig. 6a).Fe-GAP22 predicted UBER of {100} plane perfectly overlaps with DFT ones, while the Fe-GAP22 predicted UBER of {110} plane slightly deviates from DFT calculations.We further computed the normal stress during the separation (Fig. 6b) by taking the derivative of the fitted curve in Fig. 6a.The minimum stress required to separate {100} and {110} surfaces predicted by DFT are 32.68 and 34.62 GPa, respectively, and Fe-GAP22 predictions are 0.3% and 9.7% larger than DFT results.For all crack systems considered here, the deviation of the crack from the original plane to {100} or {110} is not likely to occur because it requires a lower cohesive strength than on the original plane.The Fe-GAP22 prediction of cohesive strength confirms that the cleavage will take place on the original crack plane, which is also consistent with anisotropic LEFM calculations, i.e., the largest normal stress exists on the original plane under pure mode-I loading (see Supplementary Material S11).At finite temperatures, cleavage is consistently taking place on the original plane for all crack systems, which shows that cleavage is the controlling mechanism of fracture in iron at high loading rates when considering finite temperatures (T≤300K).
At T=300K, the GAP predicted error raises to ∼30 meV/atom for (110) [1 10] and ∼10 meV/atom for the other three crack 9/17 systems (see Supplementary Material S12).Five independent K-tests are performed at each temperature for the statistic prediction of K Ic at finite temperatures, as indicated by the error bar shown in Fig. 7.Here we propose a criterion to determine K Ic for cleavage from atomistic simulation results based on the traction-separation curve.The first pair of atoms at the crack tip can still interact with each other after the initial debonding process until the entire separation is achieved, i.e. the distance between them reaches the cutoff distance in atomistic simulations.Therefore, K Ic is defined as the point where the crack tip bond is separated to a cutoff distance and will remain open upon further increasing K I .We also considered the GAP trained on the database developed based on the original SOAP descriptor (see Supplementary Material S9), referred to as Fe-S-GAP22 in Fig. 7. Inspection of the atomic snapshots reveals that the fracture mechanism for all crack systems remains cleavage up to T=300K as predicted by Fe-S-GAP22, which demonstrates that the predicted cleavage mechanism is independent of the specific crack-tip converged DFT database.The solid and dashed line indicate critical K predicted by Griffith (K G ) and Rice (K Ie ) theories, respectively.Fe-S-GAP22 is the GAP trained on the active learning database developed based on original SOAP.MD simulations at finite temperatures are performed 5 times to get the average of the predicted K Ic .
Fe-GAP22 predicted K Ic for all crack systems decreases from 0K to 1K and converges from T=1K to 200K for all crack systems except for (110)[001] (Fig. 7).The decrease from 0K to 1K is around 0.05MPa √ m, which is an accuracy limit rather than lattice trapping effects and is due to negligible barriers.Those minute barriers might not be overcome by geometry optimisation/static calculations.Instead, low temperature MD helps to overcome small barriers and move away from saddle points/symmetric configurations.Fe-GAP22 prediction of (110)[001] significantly decreases from T=0K to 10K and converges at T>100K, which is a consequence of the rough GAP predicted PES.Fig. 7b shows the Fe-S-GAP22 predicted K Ic for (110)[001] (see Supplementary Material S13 for prediction of the other 3 crack systems).The predicted K Ic 's of Fe-GAP22 and Fe-S-GAP22 converge to the same value within a 0.004 MPa √ m error bound.All predictions at T>10K are within 10% of Griffith prediction (K G ), indicating mild lattice trapping effects 41 .At high temperatures (T>100K), the surfaces near the crack tip can be closed due to thermal fluctuations, which leads to a slight increase of K Ic at T=200K and 300K for all crack systems.

Discussion
Our analysis unveils that the T=0K fracture mechanism for {100} and {110} crack planes is cleavage on the pre-existing crack planes.There are no fracture experiments performed close to T=0K, thus qualitative comparison is possible only with the work of Ref. 5 , which is performed at T=77K and finite loading rates that are inaccessible to direct MD simulations.The prediction of Fe-GAP22 qualitatively agrees with experiments on {100} crack plane, i.e., cleavage on the original crack plane.This is a considerable step forward with respect to all classical potentials considered by Möller 15 and in the introduction of the present paper, that, in contrast with Fe-GAP22, cannot predict pure cleavage on {100} plane for (100)[010] crack system.Furthermore, Fe-GAP22 is free from commonly observed artifacts of classical potentials, such as the occurrence of {110} planar faults 42 .
However, in the case of {110} crack planes, Fe-GAP22 also predicts cleavage on the original plane while experiments show crack kinking into {100} planes 5 .Our calculation of UBER by Fe-GAP22 and anisotropic LEFM analysis shows that, under K I loading and in the absence of other defects close to the crack tip, cleavage should take place on the original crack plane for both {100} and {110}.The experiments, instead, show a semi-brittle fracture behaviour of bcc iron because cleavage is accompanied by extensive dislocation activity 5 .
Therefore, in order to perform quantitative connection to experiments, the temperature and loading-rate dependent competition between cleavage and dislocation emission must be assessed, e.g. by performing Nudged Elastic Band (NEB) calculations and using transition state theory 4,43 .Also, the interaction between pre-existing defects (e.g.nanovoids, dislocations, etc) and cracks is expected to change the stress field around the crack tip and hence affect the crack-tip response.In the present atomistic fracture simulations however, the specimen is defect-free (except for the sharp crack), and the boundary conditions are those of a half-infinite crystal in a 2D setup.Therefore, quantitative predictions of fracture require a multiscale approach.
Despite the fact that a quantitative comparison of K Ic is not possible, we extrapolate the experimental K Ic linearly to T=0K.The experimental K Ic for {100} and {110} crack planes are ∼4 and 7 ∼12 MPa √ m, respectively 5 .K Ic for {110} crack plane is 2 ∼ 3 times larger than {100}, suggesting pronounced thermally-activated dislocation plasticity at finite temperatures on {110} crack plane.Such a significant energy dissipation is expected to account for the crack deviation from {110} to {100}crack plane.Thus, large-scale simulations in which rate and temperature dependence is considered along with pre-existing defects should be performed to enable prediction of the experimental K Ic .
Griffith 17 (K G ) and Rice 2 (K Ie ) theories are used to predict cleavage and dislocation emission from the crack tip.Since cleavage is the controlling fracture mechanism in bcc iron at T=0K, we compare the Fe-GAP22 predicted K Ic with Griffith theory K G .Because K G only accounts for the energy of newly created surface, it serves as a lower bound of K Ic .Thus, the excess part of Fe-GAP22 predicted K Ic compared with K G could be used as an indicator of lattice trapping effects.For all three crack systems except for (110)[110], the predicted K Ic is around 10% larger than Griffith criterion, which is also found in bcc high entropy alloys 4 .(110)[110] is 5% larger than K G , indicating less lattice trapping effects.
We obtained a Fe GAP that is able to predict atomistic fracture in iron with an accuracy of 8 meV/atom based on an existing database, showing that GAP can be improved for fracture study.Our strategy can be generalised to other materials and other ML frameworks, providing a systematic way to improve machine learning potentials for studying fracture behaviour.The current approach is based on two essential ingredients, i.e., fracture-relevant configurations and an active learning database.
Based on the fracture-relevant database, our active learning approach produced convergence to predicted energy error of 8 We conclude by pointing out that Fe-GAP22 can describe cracks and screw dislocations with quantum accuracy 25 .The Fe-GAP22 is currently being applied to calculate the energy barriers for the fracture-related processes with NEB.Based on transition state theory, the competition of atomistic crack-tip mechanisms (cleavage and dislocation emission) can be connected to microscale behaviour to assess the brittle to ductile transition 4 at finite temperatures and finite loading rates.Mobility laws for microscale simulations can thus be developed by using Fe-GAP22.

DFT calculation
All DFT calculations were performed in collinear spin-polarised plane wave method implemented within QUANTUM ESPRESSO 37 .An ultrasoft GGA PBE pseudopotential from 0.021pslibrary is employed 44 .The kinetic energy cutoff for wavefunctions and charge density are set to 90 Ry and 1080 Ry, respectively.The Brillouin Zone is integrated with a Monkhorst-Pack grid and a Marzari-Vanderbilt smearing scheme at an effective temperature of 0.01 Ry 45 .The starting value of magnetisation is set to 0.34.An adaptive k-mesh approach is used to calculate primitive cells (DB9 and DB10), ensuring that the k-spacing is smaller than 0.015 Å −1 .For other calculations, the simulation cell size is chosen to satisfy the condition that k spacing is smaller than 0.03 Å −1 while only the Γ point is used along the separation direction.A benchmark study compared with Dragoni's database is conducted to ensure all the parameter settings can reach the convergence of 1 meV/atom, 0.01 eV/Å and 0.01 GPa for energies, forces and stresses, respectively.

GAP training
A parallel version ( 1645177290) of the open-source package QUantum mechanics and Interatomic Potentials (QUIP) is used to fit the GAP model 36 .We use three descriptors to encode the local atomic environment, i.e., one distance-based two-body interaction and two many-body turbo-SOAP descriptors.Two turbo-SOAP descriptors are used to describe the inner and outer atomic environment within the cutoff radius of 3.5 Å and 5 Å, respectively.Dot product kernel is used for both turbo-SOAP descriptors, and the CUR matrix decomposition procedure is applied to find the optimal representative local atomic environments.The cutoff smoothing distance is set to 1 Å.The number of radial and angular basis for turbo-SOAP are set to 8.
The training command line is available in Data Availability.

MS/MD simulation setup
We use a cylinder-shaped half crack setup (see schematic plot in Supplementary Material S14), where x, y and z are aligned with the crack propagation direction, crack plane normal and crack front, respectively.The radius of the model in x-y plane is set to 150Å, which satisfies the requirement that the fracture process zone is much smaller than the simulation cell, as proved in Ref. 3,15 .We consider plane strain conditions in a 2D setting by imposing periodic boundary conditions along the crack front direction.
We used a K-test loading scheme that allows controlling the fracture process by directly increasing the stress intensity factor K I .K-test implements a displacement-controlled loading scheme, which makes use of anisotropic LEFM.The displacement field of a half crack in an infinite anisotropic medium can be derived via Lekhnitskii's formalism 46 , as summarised in Supplementary Material S2.1.The asymptotic stress field near the crack tip can be uniquely characterised by the stress intensity factor K I , which enables a single-parameter controlled loading scheme 3 .The K I displacement field can then be applied to the boundary of the simulation cell directly.During K-test, K I is increased with a step of ∆K I = 0.01MPa √ m.The displacement is applied incrementally at each step, and the system is equilibrated under the constraining of fixed boundary atoms.We use a combination of conjugate gradient (CG) and FIRE minimizers 47 with a force tolerance of 10 −9 eV/Å and 10 −3 eV/Å, respectively.For MD, we use the Nosé-Hoover thermostat with a timestep of 0.001 ps and the system is equilibrated for 10 ps at each incremental step.
To preserve an atomistically sharp crack, we neither screen the interaction between free crack surfaces nor delete any atoms (so-called screening and blunting 3 ).Instead, we start with a K init that maintains the current crack tip position.However, it should be noted that finding K init may require trial-and-error.All molecular statics/dynamics calculations in this work are performed by using LAMMPS 35 .The atomic configurations are visualised in OVITO 48  The displacement and stress fields of a half crack in an infinite anisotropic medium can be derived via Lekhnitskii's formulation [15].
The governing equation of the plane strain problems in an anisotropic linear elastic medium is where U is the Airy stress function and b i j are determined by the compliance S i j [16, 17]    for general plane strain problem can be expressed as where z i = x + s i y, (i = 1, 2), U 1 and U 2 are arbitrary functions to be determined by the boundary conditions.Hence, the displacements fields are expressed as where p i and q i are of the form By considering the boundary conditions of a crack in an infinite medium, one could get the stress function φ and ψ.Thus the displacement fields at the crack-tip, expressed in polar coordinates, is obtained.
The the stress fields are Since the stress and energy approaches in linear elastic fracture mechanics (LEFM) are equivalent, there is a unique relation between energy release rate G c and stress intensity factor K Ic [18].The relation can be derived with the aid of crack closure method.First, we consider a Mode-I crack before and after extension distance da.Two corresponding coordinates x − y and x − y centred at the crack tip are established (x = x − da and y = y ).According anisotropic LEFM crack-tip fields, the normal stress ahead of the crack-tip (x = x, y = 0) before extension can be obtained by plugging r = x, cos θ = 1, sin θ = 0 into σ yy of equation ( S8) Assuming crack propagation of a small distance da, the displacement along y direction of new crack surfaces (0 x da) can be calculated by plugging x = x − da into u y of equation ( S7) The strain energy associated with this process can be calculated as the work done by σ yy before extension, which closes up the crack opening after the crack extension where s 1 and s 2 are obtained by solving the governing equation (S3 (S14) q 1 and q 2 are determined by elastic constants (equation ( S6)).Then the energy release rate reads As discussed in section S2.1, B depends on the elastic constants.Since, according to Griffith theory [12], G c = 2γ s , then the Griffith prediction of the critical stress intensity factor K G is which is equation (1) in the main manuscript.

S2.3: Anisotropic Rice criterion for dislocation emission from crack tip
Equation (2) of the main manuscript predicts the critical K Ie for dislocation emission from the crack tip [13, ?, 14] where θ is the angle between slip plane and crack plane, φ is the angle between slip direction and a vector lying on the slip plane and perpendicular to the crack-front direction, and γ usf is the unstable stacking fault energy of the slip plane.F 12 (θ ) is a geometrical factor where σ i j is the near crack-tip stress field given in equation S8 where and Ω is the rotation matrix and Λ is the Stroh tensor [?] A and B can be solved by using the Stroh formalism [?, 15] Equation S25 is an eigenvalue problem where A and B are matrices formed by eigenvectors and ν is a diagonal matrix filled up by eigenvalues N is the fundamental elasticity matrix where The detailed derivation can be found in chapters 5 and 11 of Ref. [15].

S3: Summary of calculated properties for the classical potentials
Table S2 reports the potential-dependent properties that have been used to calculate K G and K Ie for the interatomic potentials considered in Fig. 1 of the main manuscript.

S6: Magnetic moment of crack-tip atoms
We performed DFT calculations of approximately periodic crack-tip cells.The magnetic moment change of each atom with respect to bulk bcc ferromagnetic iron is plotted in Fig. S2.We can see that the atoms at crack-tip remain in the ferromagnetic state, with a maximum local magnetisation change of ∼30%.GAPs with Turbo-SOAP based on the dataset developed by using original SOAP, referred to as Fe-S-GAP22.We found a better convergence in terms of GAP predicted per-atom energy error, as shown in Fig. S5.Except that the error of (100) crack systems shows slight fluctuations, the predicted errors are converged to less than 10 meV/atom for all crack systems .The predicted fracture mechanisms are the same as obtained with Turbo SOAP, suggesting the fracture mechanism is independent of database as long as the GAP predicted error is converged.

Figure 1 .
Figure 1.Summary of fracture simulations for crack system (100)[010], (100)[011], (110)[001], and (110)[110].K Ic is the critical stress intensity factor predicted by MS simulations.The letters above the symbols indicate the fracture mechanism observed in MS simulations with the different IAPs.C0 and C1 indicate cleavage on {100} and {110} planes that are different from the original crack plane.P and D indicate phase transformation and dislocation emission from the crack tip, respectively.A indicates amorphous structure forming at the crack tip.The crack propagates on the original crack plane if no symbol is specified.K G and K Ie are predictions according to Griffith 17 and Rice theories 2 for cleavage and dislocation emission, respectively.

Figure 2 .
Figure 2. Fracture predictions of Fe-GAP18 trained on Dragoni et al. original database 30 .(a) Maximum per-atom GAP predicted energy error during fracture simulations of four crack systems.The dashed line indicates the 10 meV/atom, and the arrows indicate the configuration shown in (b).(b) Simulation snapshots coloured by GAP predicted per-atom energy error.Note that this is a new version of Fe-GAP18 trained on the database by Dragoni et al. using an optimised descriptor 32 (see section Methods: GAP training).

Figure 3 .
Figure 3. (a) Schematic workflow of the active learning algorithm.(I) Fracture simulation by K-test at T=0K. (II) Evaluation of the GAP predicted per atom energy error for each frame generated during K-test .(III) Detection of extrapolation by comparing the maximum GAP predicted error per atom of each configuration with the predefined criterion.(IV) Construction of the crack-tip cell that is computable by DFT.(V) DFT calculation of the crack-tip cell.(VI) Refitting GAP.Steps (I)∼(VI) are repeated until convergence is achieved.(b) Construction of the crack-tip DFT cell.(I) Identification of the atom with the largest error (ID 0 ) in the deformed configuration.(II) Selection of a group of atoms (ID g ) in the square centred on ID 0 in the undeformed state.(III) Identification of the crack tip in group ID g in the deformed state.(IV) Duplication and rotation of the selected region.(V) Generation of DFT cell containing two crack tips by merging the rotated replica and the original ones.

Figure 7 .
Figure 7. Critical K I predicted by Fe-GAP22 at different temperatures (T=0K∼300K).(a)(100) and (b)(110) crack plane.The solid and dashed line indicate critical K predicted by Griffith (K G ) and Rice (K Ie ) theories, respectively.Fe-S-GAP22 is the GAP trained on the active learning database developed based on original SOAP.MD simulations at finite temperatures are performed 5 times to get the average of the predicted K Ic .

1 2 σ
. The work done by σ yy traversing u y equals to the energy release rate G I da under mode-I extension, leading to G c da = 2 da 0 yy u y dx.(S11) Substitution of σ yy (equation S9) and u y (equation S10) into equation (S11) results in the expression

Figure S1 :
Figure S1: (a) Principal component analysis of strain components from DB1 (original database), DB9 and LEFM.(b) Zoom-in of the data that is centred in the origin of (a).

Figure S2 :
Figure S2: Magnetisation degree of a DFT crack-tip cell.Atoms are coloured by their magnetic moment change with respect to the perfect bulk bcc ferromagnetic iron.

Figure S4 :
Figure S4: Convergence analysis of GAP trained on original SOAP.Maximum GAP predicted per-atom energy error as a function of applied K I during the fracture simulation for four crack systems.

Figure S5 :
Figure S5: Convergence analysis of GAP trained on turbo SOAP by using the database developed by original SOAP.Maximum GAP predicted per-atom energy error as a function of applied K I during the fracture simulation for four crack systems.

Fig. S6
Fig.S6shows that the prediction of elastic constants of Fe-GAP18 and Fe-GAP22 are in good agreement with each other.C 44 predicted by Fe-GAP22 is 4% off with respect to DFT.The predicted surface energies of {100} and {110} planes have a maximum difference of 1%.

Fig. S7
Fig. S7 plots the normal stress distribution as a function of the angle with respect to the crack plane for radius from 1 Å to 15 Å.The normal stress σ θ θ is calculated from equation (S16).Anisotropic LEFM (see Section S.2) predicts the maximum stress normal to the crack plane.

Figure S7 :
Figure S7: Normal stress (σ θ θ ) distribution around the crack tip computed using anisotropic LEFM.r is the radius.

S12:
GAP22 predicted error at T=300K during fracture simulation We performed fracture simulation of four crack systems at T=300K under a loading rate of 10 9 MPa √ ms −1 .Fig. S8 plots the maximum GAP prediction energy error per-atom as a function of K I for four crack systems.For all crack system expect for (110)[110], the maximum error increases to 15 meV/atom due to the thermal fluctuation.For (110)[110], the maximum error is around 30 meV/atom.

Figure S9 :
Figure S9: Critical K I predicted by Fe-S-GAP22 at different temperature (T=0K-300K).(a)(100) and (b)(110) crack plane.The solid and dashed line indicate critical K I predicted by Griffith (K G ) and Rice (K Ie ) theories, respectively.

Figure S10 :
Figure S10: (a)Schematic plot of atomistic fracture simulation setup.x, y and z are aligned with the crack propagation direction, crack plane normal and crack front, respectively.(b) K-test iteration algorithm.
instead of the strain states directly obtained by classical elastic crack-tip fields (see Supplementary Material S4).Principal component analysis of the strain components

11/17 meV
/atom in 4 iterations.We found that crack propagates on the original crack plane for (100)[010], (100)[011], (110)[001]and (110)[110] crack system.Based on the convergence of Fe-GAP22 predicted K Ic with respect to temperatures, we conclude that the increase of K Ic at T=0K is caused by minute barriers.The predicted K Ic at T=0K is converged with respect to local symmetry breaking, which we have studied by performing high loading rate mode-I tests from T=10K up to T=300K.Our results confirm the T=0K atomistic fracture mechanism in ferromagnetic iron under mode-I loading, settling the inconsistency of crack-tip mechanisms based on classical potentials.

Table S1 :
. Method summary of existing atomistic studies.

Table S2 :
Summary of calculated properties for the classical potentials used in Fig.1.(a 0 : lattice constant.C ij : elastic constants, γ s : surface energies, γ usf : unstable stacking fault energies)