## 2 Calculation Methods
### Elastic Constants
 \tag{5}\]  
where \(C_{ij}(0)\) is calculated by using MS according to Equation (3), and its structure is established by periodically stacking the unit cells of \(\alpha\)-Fe along the three axes 52 times each, with a total of 281216 iron atoms. The cell lattice constant is 2.85573A. \(C_{ij}^{q}(c_{H})\), which is calculated using a similar method, but with randomly inserted hydrogen atoms (with concentrations of \(c_{H}\in\) [0.0002, 0.01]) in the sites of the tetrahedral interstitials (Tet-sites) using the Knuth-shuffle algorithm [40] after stacking the super cell when making structures. The changes of volume and its influences on the elasticity of \(\alpha\)-Fe can be connected by temperatures, as the variations of volume with temperature (T) (expressed as \(\Delta V(T)=V(T)-V(0)\)) can be calculated by relaxing the structures under different temperatures in the isothermal-isobaric (NPT, N, P, and T denote the given particle number, pressure, and temperature, respectively) ensemble [41; 42]. Each structure consists of \(\alpha\)-Fe unit cells stacked for 10 cycles along the three axes, and a total of 2000 iron atoms. With the time step \(t_{s}=0.1\,f\)s, the equilibrium structure S\({}_{\rm e}\) and the system volume V(T) were obtained after running 50,000 steps at each temperature T. We use Equation (6) to fit the results:  
\[\Delta V(T)=y_{1}+b_{1}e^{-T/b_{2}} \tag{6}\]  
where \(y_{1}\), \(b_{1}\), and \(b_{2}\) are the parameters to be determined. The elastic constants at finite temperatures were obtained using the fluctuation method. Using S\({}_{\rm e}\) as the initial structure, when \(T\) is not too high, there is [43]:  
\[C_{ijkl}(T)=\langle\varepsilon_{ij}\sigma_{mn}\rangle\langle\varepsilon_{mn} \varepsilon_{kl}\rangle^{-1} \tag{7}\]where \(i,j,k,l\), \(m\), \(n=1\), \(2\), \(3\). Given strain \(\varepsilon_{j}\), Equation (7) is simplified as follows for the body-centered cubic (bcc) structures:  
\[C_{ij}(T)=\langle\alpha_{i}\rangle/\varepsilon_{j} \tag{8}\]  
where \(i,j=1\), \(2\), \(4\). Accordingly, \(C_{ij}(T)\) can be calculated in the \(N\hbar T\) ensemble (\(N\), \(h\), and \(T\) denote the given particle number, scale matrix, and temperature, respectively) [44]. In this work, let \(t_{s}\)= 0.1 \(f_{s}\), \(T\)\(\in\)[\(3\), \(500\)] K, \(S_{\rm e}\) is relaxed under the \(N\hbar T\) ensemble [41, 42] for 10,000 steps, and strain \(\varepsilon_{j}=0.01\) is applied in each direction. After each strain was applied, 3000 steps were run under the \(N\hbar T\) ensemble. The variations of \(C_{ij}\) with temperature \(T\) was fitted by Equation (9):  
\[\Delta C_{ij}^{temp}(T)=C_{ij}(T)-C_{ij}(0)=y_{2}^{ij}+b_{3}^{ij}e^{-T/b_{4}^ {ij}} \tag{9}\]  
where \(y_{2}^{ij}\), \(b_{3}^{ij}\), and \(b_{4}^{ij}\) are the undetermined parameters. Substituting them to Equation (6), we obtain:  
\[\Delta C_{ij}^{vol}(\Delta V)=y_{2}^{ij}+b_{3}^{ij}\bigg{(}\frac{\Delta V-y_{ 1}}{b_{1}}\bigg{)}^{b_{2}/b_{4}^{ij}} \tag{10}\]  
Averaging the results of many structures under the same \(c_{H}\) eliminates \(\Delta C_{ij}^{dis,q}(c_{H})\); thus, the \(\Delta C_{ij}^{coh}(c_{H})\) can be calculated from Equation (11):  
\[\Delta C_{ij}^{coh}(c_{H})=\frac{1}{N}\sum_{p=1}^{N}\bigg{(}C_{ij}^{p}(c_{H}) -C_{ij}(0)-\Delta C_{ij}^{vol}(\Delta V^{p})\bigg{)} \tag{11}\]  
where \(N\) is the total number of structures and \(\Delta V^{p}\) is the equilibrium volume of structure p. First-principle calculations show that, for maintaining structural stability, the space between hydrogen atoms in the equilibrium structure of iron should be greater than 2.1 A [45], or greater than the fourth-nearest neighbor (4 NN) distance between the Tet-sites that can be occupied by hydrogen atoms [27]. Therefore, we generated hundreds of structures under each \(c_{H}\) for relaxation and screening, so that the number of calculated structures at each \(c_{H}\) was no less than 50. At \(c_{H}\) = 0.8% we obtained 101 structures, which was the most of all \(c_{H}\).