**2. Method**

In this work, molecular dynamics (MD) simulations are performed by Large scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [24] which is an open source MD software to simulate atomic interactions for a given system. The output data are analyzed by OVITO software (Version 3.5.4, OVITO GmbH, Darmstadt, Germany) [25]. Meanwhile, as stated above, 7 different grain boundaries including ∑3(112), ∑3(111), ∑5(012), ∑5(013), ∑9(221), ∑11(113) and ∑17(410) are built in bcc Fe based on coincidence site lattice (CSL) theory [26]. The GB angle of these GBs is from around 36◦ to 109◦. Periodic boundary condition (PBC) is applied along 3 directions. In order to avoid the interaction between two GBs induced by PBC, the distance along the normal direction of GB is set to at least 14.8 nm, as listed in Table 1. A schematic of the symmetrical GB model is shown in Figure S1, provided in the Supplementary Materials. The other information for each simulation box is also included in Table 1. Atomic structures of these 7 GBs before full relaxation are shown in Figure 1. The Fe-Fe interaction is described by Fe potential developed by Mendelev et al., [27]. This potential has been well applied for grain boundary simulations [5,13,28] and is suitable for the present purpose.

**Table 1.** Parameters of symmetrical GBs used in the present work.


**Figure 1.** GB atomic structures before and after MD relaxations shown in the top and bottom panels, respectively: (**a**) ∑3(112) GB, (**b**) ∑3(111), (**c**) ∑5(013), (**d**) ∑5(012), (**e**) ∑9(221), (**f**) ∑11(113) and (**g**) ∑17(410). Structure types are analyzed by the common neighbor analysis method. Atoms near the GB region are colored white, while the atoms in the grain are colored blue.

Based on the above simulation models, conjugate-gradient (CG) method is used for relaxation at 0 K, and subsequently followed by MD relaxation at 300 K and 600 K. After full MD relaxation, the CG method is applied again to obtain the atomic structure and energy at 0 K. In this way, the local structure may overcome the energy barrier and reach a lower energy state after the CG-MD-CG relaxation process. In CG relaxation, the specified energy tolerance is 10−<sup>10</sup> and the specified force tolerance is 10−<sup>10</sup> eV/Å3. During the MD relaxation, the timestep is 1 fs and at least 20 ps relaxation is applied for each process to ensure the system is fully relaxed at the given temperature. It should also be noted that during relaxation, the constant number of atoms, pressure and temperature (NPT) ensemble is applied in order to relax both the atomic position and simulation volume. Furthermore, for volume relaxation, each direction of the box is allowed to relax independently to fully release the internal elastic stress and obtain a more stable state for further calculations. The GB formation energy, *EGB*, is then calculated according to the following equation as applied in previous studies [13,29,30], which is defined as the difference between the potential energy *Etotal* of *n* atoms in the supercell containing GBs and the potential energy of a computational cell with the same number of atoms in a perfect crystal, divided by the cross-sectional area, *S*, of two GB planes.

$$E\_{GB} = \frac{E\_{total} - NE\_p}{2S} \tag{1}$$

where *Etotal* is the total energy of a system containing 2 GBs. *N* is the number of atoms in this system. *Ep* is the cohesive energy of one atom in a perfect bcc Fe crystal, which is −4.12 eV according to Mendelev potential [27].

To investigate the mechanical property of the grain boundary, after full relaxation, external strain is applied to the system along the normal direction of GB plane with a strain rate of 5 × <sup>10</sup>9/s. The simulation temperature is set at 300 K. It should be noted that during the strain application along the normal direction of GB plane, the box length along the other two directions of the computational box is allowed to change automatically to ensure zero pressure along these two directions. The stress related to the increases in strain is also calculated and thus, the stress–strain curve is obtained. For comparison, this process is also performed for a perfect structure with strain direction along the same normal direction of GB plane. During the strain application along GB normal direction, the length of the box along the directions perpendicular to GB normal direction is also allowed to vary, which is similar to the real experimental condition and ensures the system is at a single stress–strain state [31].

To investigate the underlying factors affecting the mechanical property of GB, different methods are used to analyze GB structures and properties including the common neighbor analysis (CNA) [32], free volume (FV), single-atom potential energy distribution, local stress field and the single-atom displacement magnitude. In this work, CNA analysis, atomic energy distribution and displacement are performed by OVITO software and details can be found in [25]. The free volume is obtained by introducing enough grids along three directions in computation box and calculating the maximum distance (*Dmax*) of each grid point to surrounding atoms. The free volume is then calculated as the sphere volume with radius of *Dmax*-*r*, where *r* is the atomic radius in perfect lattice. This method has been well used in previous studies [33]. The FV calculation under the application of external stress was performed for polymer and amorphous materials [34], which have an FV distribution through the whole system. In this work, this method is also applied firstly for GB investigation since the FV in GB is expected to change in the GB failure process. Furthermore, as stated in the Introduction, the results from the present work may be used for comparison with neutron diffraction measurements under the effect of external stress

field in future. The stress of each atom is calculated according to the following equation and can be viewed via the Ovito software:

$$
\sigma\_{ij}^V = \frac{1}{V} \sum\_{a} \left[ \frac{1}{2} \sum\_{\beta=1}^{N} \left( R\_i^{\beta} - R\_i^a \right) F\_j^{a\beta} - m^a \upsilon\_i^a \upsilon\_j^{\beta} \right] \tag{2}
$$

where (*i,j*) take values of x, y and z (directions). *β* takes values 1 to *N* neighbors of atom *α*. *Rα <sup>i</sup>* is the position of atom *<sup>α</sup>* along direction *<sup>i</sup>*. *<sup>F</sup>αβ <sup>j</sup>* is the force (along direction *j*) applied on atom *β* from atom α. V is the total volume of the system. *m<sup>α</sup>* is the mass of atom *α* and *v<sup>α</sup> <sup>i</sup>* is the thermal excitation velocity of atom *α* along *i* direction.
