Hybrid quantum/classical study of hydrogen-decorated screw dislocations in tungsten: Ultrafast pipe diffusion, core reconstruction, and effects on glide mechanism

The interaction of hydrogen (H) with dislocations in tungsten (W) must be understood in order to model the mechanical response of future plasma-facing materials for fusion applications. Here, hybrid quantum mechanics/molecular mechanics (QM/MM) simulations are employed to study the 〈111〉 screw dislocation glide in W in the presence of H, using the virtual work principle to obtain energy barriers for dislocation glide, H segregation, and pipe diffusion. We provide a convincing validation of the QM/MM approach against full DFT energy-based methods. This is possible because the compact core and relatively weak elastic fields of 〈111〉 screw dislocations allow them to be contained in periodic DFT supercells. We also show that H segregation stabilizes the split-core structure while leaving the Peierls barrier almost unchanged. Furthermore, we find an energy barrier of less than 0.05 eV for pipe diffusion of H along dislocation cores. Our quantum-accurate calculations provide important reference data for the construction of larger-scale material models.

atoms show a strong attraction to screw dislocations resulting in hardening by a pinning mechanism, while Hf, Ta, and Re reduce glide barriers and thus facilitate plastic slip [11]. Re is known to affect the stable core structure as well as the glide plane of screw dislocations in W [12,13], while interstitial carbon stabilizes the hard-core configuration in Fe, Mo, and W [14,15].
In a fusion environment, a PFM will be continuously populated with hydrogen isotopes. A full picture of the interactions of H with dislocations is essential for understanding PFM ageing, in particular the H isotope retention phenomenon [16]. H and He behavior in W has been investigated by means of atomistic simulations; thorough reviews of recent modeling activities can be found in Refs. [17][18][19]. According to DFT studies, H atoms occupy tetrahedral interstitial positions in bulk W with a migration barrier of around 0.25 eV [20][21][22][23].
Since the dislocation glide is governed by a very localized rearrangement of atomic bonds at the dislocation core, accurate atomistic methods are required to capture the migration mechanism [8]. While empirical interatomic potentials are constantly improving, the accuracy and transferability of existing interatomic potentials in application to dislocations remains limited [24,25], especially when accounting for chemical effects induced by impurities, meaning that ab initio simulations are essential. However, dislocations have a net elastic disregistry which generates long-range elastic fields whose accommodation in principal demands large model systems beyond the computational size limits of ab initio methods. Combining a quantum mechanical description of FIG. 1. Visualizations of a fully periodic 135-atom quadrupole supercell (a) and a 237-atom cluster supercell (d), with atoms colored by displacement along 111 direction. The purple dashed circles show the radius of the unconstrained region for the cluster cell (a) and the circle with a diameter equal to the distance between two dislocations with opposite Burgers vectors for the quadrupole cell (d). Two cells were considered equivalent if the sizes of these two regions were equal (see Sec. III A). The red triangles in (a), (d), and (e) show two equivalent easy-core positions. Solid and empty triangles show initial and final core positions used for NEB calculation of the dislocation glide MEP with a red dashed line corresponding to the initial linear guess for the MEP. (e) Differential displacement map corresponding to the initial easy-core configuration. A red star and red cross show hard-core and split-core positions, respectively. The structure of a perfect bcc crystal in the (111) plane and (b) in the (110) plane (c). Left part of (f) shows the same bcc structure as (c), with colored arrows showing the 111 displacement for the easy core. Middle part of (f) shows resulting easy-core structure. Right part of (f) shows hard-core structure, which can be obtained by reversing the displacement. Red, green, and blue colors of W atoms in (b), (c), (e), and (f) show three nonequivalent atomic planes along 111 direction in bcc structure. Note change of ordering of W atoms along 111 direction from red-green-blue for ideal bcc structure to green-red-blue for easy core and the single horizontal red-green-blue plane for the hard core (f). the dislocation core (along with any nearby impurities) where bond rearrangement takes place with a classical description of the long-range elastic relaxation is thus an appealing prospect.
A popular solution is a cluster arrangement, where a dislocation core is contained in a cylindrical disk of atoms governed by ab initio methods oriented perpendicular to the 111 Burgers vector direction, which is then coupled to a linear elastic medium controlled by the lattice Green's function [26][27][28] able to accurately capture the far-field deformation. Another option is to embed the ab initio cluster in an atomistic region governed by empirical force fields. The embedding region is then made large enough in order to accommodate far-field deformation. This latter technique, used here, is an example of a hybrid QM/MM simulation [29,30]. In all cluster methods, an additional layer of redundant "buffer" atoms is required in the ab initio cluster simulation to mitigate electronic effects induced by the artificial free surface of the cluster. It is well known that, due to the delocalization of the electronic energy density, cluster-based simulation methods possess well-defined forces but no energy function. Recently, a partial solution to this problem was developed, allowing the evaluation of energy barriers from hybrid QM/MM simulations using the principle of virtual work to integrate the projected force along pathways in configuration space [31].
To account for the long range elastic fields, the empirical atomistic region for dislocation simulations can either be an axially periodic cylinder with the outer boundaries of the disk fixed according to linear elastic solution for an isolated dislocation [32] (so-called cluster approach), shown in Fig. 1(d), or a fully periodic supercell with an additional dislocation of opposing Burgers vector, creating an infinite dipole array with zero net disregistry [27,33]. The dislocation dipoles can be arranged in either a honeycomb or square lattice; a square lattice was found to be the most effective for reducing elastic interactions between dislocations [34], referred to as a quadrupole arrangement. An example quadrupole supercell is shown in Fig. 1(a). While in principle it is possible to create a dislocation dipole in a supercells of a few hundred atoms small enough to be treated in fully periodic density functional theory (DFT) calculations, where a total energy is defined, the strong elastic interactions between dislocations in the dipole and their periodic images is typically too large to give useful results. However, 111 screw dislocations in bcc metals are one of the small number of cases where the interactions are sufficiently weak that such simulations are meaningful [8,27,33,34], although errors can be introduced when a defect or an impurity is present near one of the dislocations due to the broken symmetry.
The three nonequivalent atomic planes along the 111 direction in bcc materials [shown with red, green, and blue circles in Figs Burgers vector it is possible to obtain two types of dislocation cores: "easy" and "hard" cores [35][36][37][38].
In our calculations we applied displacement fields corresponding to positive Burgers vector centered on an upwardpointing triangle in order to obtain an easy-core configuration, shown by a red solid triangle in Fig. 1(e). This results in reversed chirality of the triangle containing the dislocation core. This effect is illustrated in the middle part of Fig. 1(f). The hard-core configuration can be obtained by application of the same displacement to an downward-pointing triangle, shown by a red star in Fig. 1(e). In this case the applied displacements cancels out the natural chirality of the three nonequivalent planes in 111 direction resulting in the arrangement of core atoms in straight lines perpendicular to 111 direction. The resulting hard-core structure is shown in the right part of Fig. 1(f).
Previous DFT calculations of 111 screw dislocation dipoles in bcc metals have found that the easy core is symmetric, nondegenerate, and the most stable [10,26,32,39]. Glide migration thus occurs between two easy-core positions [34,37,40,41] shown with (filled and open) red upwardpointing triangles in Fig. 1(e). The minimum energy path (MEP) of the dislocation core passes close to a third type of core configuration-the "split" configuration when the dislocation core is centered close to an atomic column [red cross in Fig. 1(e)] resulting in a single hump-shaped glide MEP profile [37,[41][42][43][44][45] presented with dashed lines in Fig. 2. However, most interatomic potentials for tungsten incorrectly predict a degenerate easy core with a double-hump glide profile [24,25]. Recent bond order potentials (BOPs) [46] as well as embedded atom method (EAM) potentials [42] are in better agreement with DFT in terms of core stability and the shape of the glide MEP (cf. brown solid line with filled squares in Fig. 2). However, the necessity to fit a new potential for every new impurity element in the system remains a limiting factor for studying the effect of impurity-dislocation interactions.
Interaction of H with dislocations was mainly studied in the context of H isotope retention in W after plasma exposure. DFT studies of H and He interaction with screw dislocations in tungsten using quadrupole geometries showed attractive interaction with possible acceleration of pipe diffusion of H along screw dislocations [16,47]. However, the effect of H or He on the dislocation core structure and glide properties have not been investigated by DFT calculations and as such definitive conclusions is precluded by the known unreliability of empirical potentials in this regime. An energy-based QM/MM scheme [48] was applied to study H and He interaction with screw and edge dislocation in Fe, where it was demonstrated that both impurities are attracted to both types of dislocations with low diffusion barrier along the screw dislocation core [49]. H was found to affect the structure of the screw dislocation core together with lowering of the glide barrier. However, neither the dislocation glide path in the presence of H nor how the effect of H on the glide barrier depends on the energy correction scheme used during QM/MM coupling were reported.
The first aim of this paper is to provide a valuable validation of the QM/MM virtual work principle for dislocation systems by comparing calculations of the Peierls barrier [7] for 111 screw dislocations in W using both the QM/MM virtual work method and small periodic quadrupole cells where a total energy function is available (Sec. III B). After validation, the QM/MM approach is applied to study interaction of H with dislocation cores (Sec. III C) together with the effect of H on the glide mechanism (Sec. III E).

A. Description of the QM/MM scheme
A QM/MM simulation always defines three principal regions: a QM region, here containing a dislocation core, a buffer region surrounding the QM region, and the remaining bulk region. On each force call a DFT calculation is performed on a cell containing atoms from QM and buffer regions surrounded by vacuum. The buffer region is essential to reduce the effect of artificial free surfaces of the cell on the accuracy of the forces in the QM region. Forces on the buffer atoms obtained by DFT are ignored during minimization, with forces on buffer and bulk atoms given by a classical interatomic potential, here the "EAM3" embedded atom model from Ref. [42], resulting in the "abrupt buffered force mixing" QM/MM scheme described in Ref. [50]. The advantage of the scheme is accurate forces throughout the system including near the QM/MM interface. However, energy per atom in the QM region is not accessible due to the nonlocal nature of electronic energy and thus the value of the total energy of the system is not available. Therefore, the virtual work principle described in Ref. [31] was employed in order to extract the energy barriers of interest.
To briefly summarize the virtual work approach, in a system of N atoms we define a (possibly unconverged) pathway X(λ) ∈ R 3N , where λ ∈ [0, 1] is an affine parameter along the path. For any given configuration X(λ) we also have access to the QM/MM force vector F(λ) ∈ R 3N . The virtual work principle states that the energy difference along the pathway where d dλ X(λ ) is the pathway tangent at λ . Using the chain rule it is simple to show that any nonlinear mappingλ(λ) leaves the energy barrier max λ E (λ) unchanged.
The scheme described above was implemented within the framework of the atomic simulation environment [51] Fig. 1(d). The buffer size was chosen to be 12 Å, resulting in a DFT cluster containing a total 144 atoms for pure DFT simulation. Convergence tests for the glide barrier using buffer sizes up to 14 Å and QM regions containing up to 23 QM atoms showed a relative difference of less than 10% (≈7 meV), comparable to the error bars associated with NEB convergence (cf. Sec. II B below) and in agreement with previous tests of our QM/MM approach for fcc Al and bcc Mo [31]. The PBE generalized gradient approximation [54] was used to describe effects of electron exchange and correlation together with a projector augmented wave (PAW) basis set with a cut-off energy of 550 eV. The Brillouin zone was sampled with a 1×1×12 Monkhorst-Pack k-point grid, reflecting the fact that the cluster is periodic along the dislocation line. Occupancies were smeared with a Methfessel-Paxton scheme of order one with a 0.1 eV smearing width. The values of these parameters were chosen after a series of convergence tests on forces with a tolerance of few meV/Å. The lattice parameter a 0 = 3.17 Å and elastic constants C 11 = 540 GPa, C 12 = 204 GPa, and C 44 = 142 GPa obtained from VASP calculations with the above mentioned parameters were used to generate an initial displacement field using Stroh method [55] as implemented within the atomman [56] package. The EAM3 potential from Ref. [42] was used for force calculation in the MM part of the hybrid scheme. The potential was rescaled to match the DFT lattice parameter and bulk modulus as described in Ref. [31]. Endpoints for each MEP were relaxed with a maximum force tolerance of 0.01 eV/Å, while a force tolerance of 0.05 eV/Å was applied for nudged elastic band (NEB) [57] MEP calculations. The FIRE minimization method was used in all cases [58].
The periodic dislocation quadrupole configurations [ Fig. 1(a)] for full DFT simulations were generated by our own implementation of the method described in Refs. [34,37,59], which is available as part of the matscipy [60] package. For the quadrupole case, the NEB glide calculation was carried out by moving both dislocation cores simultaneously and dividing the resulting barrier by two.
No elastic correction for the resulting energy barrier was applied, though previous studies have shown this correction is minimal [41,44,45], as we confirm below.

B. Error propagation for the virtual work approach
DFT codes such as VASP use a self-consistent field (SCF) procedure to determine the electronic ground state, leading to DFT atomic forces with some residual error, in addition to that induced by choice of exchange-correlation functional and other parameter choices. As a result, the atomic force vector F(λ) used in the virtual work calculation of Eq. (1) will not be the "true" vector of Hellman-Feynman forces F 0 (λ) for a given configuration X(λ). To estimate the resulting error on energy barrier calculations, we model the QM/MM force vector F(λ) for a configuration X(λ) as where η(λ) is a white noise vector inside the QM region and zero otherwise and σ is a RMSD force error for the QM region. From Eq. (1) it is simple to show that the variance on the energy difference is given under these assumptions by where I QM is a diagonal matrix which is unity for atom coordinates inside the QM region and zero outside. To evaluate σ , we monitor the variation in the force vector across the QM region during the last stages of minimization, for all configurations whose maximum force component is less that 0.05 eV/Å. We have also evaluated the error in the energy barrier by sampling new forces according to the per-coordinate force distribution from the same configurations then integrating these sampled forces directly to sample the energy difference. Both cases gave essentially identical results, as represented by the error bars and filled ribbons in Figs. 2, 4, and 5.

A. Quadrupole vs cluster size tests using an empirical potential
In order to estimate the effect of the finite size of the overall cell, a number of NEB calculations for the glide barrier of a 1 2 111 screw dislocation using cluster and quadrupole approaches were performed. In both cases we considered an easy-core configuration (Fig. 1). The "EAM4" potential from Ref. [42] was used. A quadrupole cell containing 135 atoms and a cluster cell containing 237 are shown in Figs. 1(a) and 1(d), respectively. The dislocation core positions were extracted by fitting to an analytical solution for the displacement field [41]. The purple dashed circle shown on the quadrupole cell has a diameter equal to the distance between two dislocation cores. The size of the cluster cell was chosen to have the same diameter of unconstrained atoms around the dislocation core. Larger cells were considered to be of equivalent sizes when the same condition was satisfied. The values of the glide barrier obtained with different sizes were compared to reference configurations with core distance (cluster radius) of 818 Å containing 377 394 atoms in cluster configuration and 235 791 atoms in quadrupole configuration. The glide barriers for both reference configurations is 0.0621 eV with FIG. 3. Size dependence of relative error in the glide barrier for quadrupole and cluster configurations computed with EAM4 interatomic potential from Ref. [42]. The reference value of 0.0621 eV was obtained from a cluster with radius 818 Å containing 377 394 atoms. a relative difference of ≈10 −4 (≈10 −5 eV) between cluster and quadrupole configurations.
As can be seen from Fig. 3 the smallest quadrupole configuration containing 135 atoms provides a relative error of a few percent. At the same time, the error for the corresponding cluster configuration containing 237 atoms is of order of 40%. It is important to note that the 135-atom cell is the most widely used for DFT calculations and configurations with larger sizes are used only for convergence tests due to their high computational cost. In order to achieve the same accuracy of a few percent with a cluster configuration, one has to use a cell containing at least 1200 atoms, which is impractical with a full DFT calculation. Moreover, this estimation does not take into account the effect of free surfaces on DFT forces in cluster configuration which can increase the total error for a given size. We note that, since the accuracy of a quadrupole cell relies on error cancellation resulting from the symmetry of two dislocations of opposite signs, adding an impurity to one of the cores reduces this symmetry and thus deteriorates the accuracy. Adding an impurity to both cores is possible in order to keep the symmetry, however it adds extra complexity and affects the overall stability of the system. This result demonstrates that a large cluster cell is thus a more suitable choice for studying interactions of dislocation core with impurities and/or defects or clusters of defects. In this work we use a cluster configuration containing 1437 atoms for calculations using a hybrid QM/MM scheme.

B. Screw dislocation core structure and glide in pure W
The results of the hybrid QM/MM NEB cluster calculation of the screw dislocation glide barrier in pure W are shown with a solid blue line and circles in Fig. 2. The brown line with solid squares shows the results of the same calculation, but where the system was fully treated with EAM3 potential from Ref. [42] used for the MM part in QM/MM scheme. Dashed lines with hollow circles give previously reported DFT results with periodic quadrupole configurations containing 135 atoms [42,44,45]. Our results demonstrate that treating only 23 atoms around the dislocation core at the DFT level allows us to obtain an MEP and absolute barrier height in much better agreement with reference DFT results, within the range of the reported DFT data (≈10 meV). It is possible to reduce the spread of the found migration barriers by using improved NEB minimization algorithms [61].
The structure of the dislocation core corresponding to the initial, saddle, and final points of the MEP are shown as insets in Fig. 2, with the dislocation core migration path shown with a red dashed line in each inset. The optimized dislocation core path significantly differs from the initial linear interpolation guess shown in Fig. 1(e) and shows a curved path between two easy-core configurations (left and right blues boxes). The saddle point is located close to the "split" core configuration in agreement with previous DFT calculations using quadrupole configurations [37,42,44,45].

C. H segregation to the screw dislocation core
A number of QM/MM relaxations were performed in order to investigate the effect of H atom on stable dislocation core structure and position. The starting configurations for each relaxation were obtained by adding an H atom to a relaxed cell containing a dislocation, with initial position given by adding the screw dislocation displacement field to a tetrahedral interstitial position as shown with blue diamonds in Fig. 4(a). Since the possible effects of H interstitials are most interesting in the context of dislocation glide, two sets of simulations were performed: with a dislocation in the initial glide position [brown triangle in Fig. 4(a)] and with a dislocation at the final position [orange triangle in Fig. 4(a)]. The relaxed configurations were then used as initial and final configuration for NEB calculations of H migration as well as the dislocation glide in the presence of an H atom.
The resulting stable H positions at the dislocation core are shown with colored diamonds in Figs. 4(f) and 4(g). The color scale represents clusters of atomic positions with distance lower than b/2. Orange diamonds show the compact cluster of positions inside the triangle of W atomic columns forming the dislocation core. Purple, cyan, and brown diamonds show three groups of positions situated outside of the dislocation core. Due to rotational symmetry these positions can be considered as one group. These findings are in agreement with pure DFT studies using a periodic quadrupole cell [16,47]. Figure 4(b) shows the effect of an H atom situated inside the dislocation core on the core structure and position. Neither the core position (red triangle) nor the differential displacement map are affected by the H atom; the dislocation core remains a symmetric nondegenerate easy-core structure situated in the center of the same atomic triangle. In contrast, Figs. 4(c) and 4(d) show the effect of an H atom situated outside the dislocation core: Fig. 4(c) shows the core has moved to the left towards the H atom, and at the same time the corresponding differential displacement map has lost symmetry and is close to the split-core structure. A similar effect is seen in Fig. 4(d), where the dislocation core has moved significantly towards the H atom and now is close to the saddle point for the pure W core glide MEP (red dashed line). Overall, we find that H stabilizes dislocation cores in configurations close to the split core that corresponds to the saddle point for dislocation glide in pure W, suggesting that the presence of H may affect the dislocation glide mechanism.

D. H migration around the screw dislocation core
The H atoms represented with orange diamonds in Due to the threefold symmetry of the triangular tessellation formed by the lattice in the (111) plane, positions marked by cyan, purple, and brown diamonds in Fig. 4(f)  The first class of barriers represents pipe diffusion along the dislocation line, also reported in [47]. The second class barrier represents the energy barrier for H to leave the dislocation core towards bulk W. As can be seen from the figure, the in-core (orange) migration barriers are almost three times lower than the barriers for leaving the dislocation core (green, purple). Moreover, the value for the migration barrier of 0.05 eV is significantly lower than H bulk migration barrier of 0.25 eV [20][21][22][23] suggesting pipe diffusion along the dislocation line is favorable once an H atom becomes attached to a dislocation core. A similar but less pronounced effect was found in Ref. [47] using pure DFT simulation with a periodic quadrupole cell.

E. Screw dislocation glide in the presence of H
Our observation from Figs. 4(c) and 4(d) shows that adding an H atom shifts the stable dislocation core towards the splitcore configuration. This core configuration corresponds to a saddle point for the dislocation glide in pure W (see Fig. 5), thus it is interesting to investigate a possible effect on dislocation glide mechanism. We carried out NEB calculations using relaxed configurations with H attached to the dislocation core as initial and final images, using a number of possible combinations of initial and final positions. First we looked at a migration barrier when the initial and final structure is not affected by the H atom shown in Fig. 4(b). Here the NEB calculations failed to converge to a reasonable reaction path resulting in unrealistic core MEPs with energy differences along the path larger than 0.3 eV. A possible explanation is that the configuration when both dislocation core and H atom are situated in the downward-pointing triangle of W atoms close to the hard-core configuration [see Fig. 1(e)] is highly unstable. Indeed, the configurations shown in Figs. 4(c) and 4(d) demonstrate that the dislocation core remains within the initial upward-pointing triangle of W atoms after relaxation. All calculations with reaction paths passing close to the hard core failed to converge resulting in barriers larger than 1.0 eV.
The results of NEB calculations for the glide in the presence of H are shown in Fig. 5 with an orange line and are compared to the results for the glide in pure W shown with the blue line (reproduced from Fig. 2). The orange-and blue-edged insets in Fig. 5 show the differential displacement maps and core position together with H positions corresponding to initial, final, and saddle points of the glide barriers. These snapshots demonstrate that the reaction path does not involve H atoms positioned in the adjacent downwardpointing triangle of W atoms.
We find that considering the relevant dislocation core configurations helps to find a physical interpretation for the MEP trajectory. As a further example, the lower orangeedged insets in Fig. 5 demonstrate that in the presence of H, glide occurs between two configurations with core structures close to the split-core configuration. These configurations are similar to the saddle point configuration for the glide in pure W shown on the central blue-edged inset. At the same time, in the presence of H the saddle point shifts to an easy-core configuration (central orange-edged inset), identical to the initial and final configurations for the glide in pure tungsten (first and last blue-edged insets). The height of the barrier is slightly higher when H is present in the material, however, this difference is not significant when compared with our estimated error bars. The presence of H also affects the shape of the MEP, which flattens near the saddle point.

A. Comparison to elasticity theory
We have calculated the point defect tensor P H ∈ R 3×3 for the H defect using the method described in Ref. [62]. To a high degree of accuracy we found that P H = P H I 3 , where P H = 5.67 eV. Using the strain field given by anisotropic elasticity theory (r) at a position r away from the dislocation core, the elastic estimation of the interaction energy is obtained as E el (r) = P H Tr (r).
We note that in isotropic elasticity theory Tr (r) = 0 for screw dislocations, thus predicting E el (r) = 0; even when accounting for elastic anisotropy, the predicted interaction energy is still very weak outside of the core region, where elasticity theory is expected to be valid. In Fig. 6 we calculate the predicted change in binding energy for all H positions considered and compare these to the results from force integration. Matching the energy differences furthest from the core, we see reasonable agreement up to the first interpolation point then a significant deviation. This is to be expected as elasticity theory is not able to capture core strain fields or migration barriers.

B. Comparison with experiments
Contamination of tungsten with hydrogen has been studied experimentally. However, the main focus of those studies are the reaction of the material to the heat loads and thermal shock [3][4][5] and H retention [16], as well as interrelations between these phenomena. To fully elucidate the implications of the energy barriers found in this work with experimental data would require mesoscale simulations such as cluster dynamics [63,64], reaction rate theory [65,66], or kinetic Monte Carlo [67,68]. However, these simulations must also account for all other known processes involving H in W, including segregation to void and other interstitial clusters, which is beyond the scope of the current contribution. We therefore leave rigorous parametrization of a full mesoscale model for future work.

V. CONCLUSION
We have carried out QM/MM calculations of the structure and energetics of screw dislocations in tungsten, both with and without the presence of hydrogen. In pure tungsten, our work closely agrees with previously reported results using a full DFT description within a periodic quadrupole unit cell. Given the very small ∼100 meV energy barriers involved, this result provides compelling validation of our QM/MM virtual work approach. A further advantage of our approach is that an error bar on the output quantity of interest can be estimated; here these are found to be of the order of 10 meV, comparable to the accuracy of the DFT calculation itself.
When attached to a dislocation line, an H atom shifts the most stable dislocation core structure from easy to split-core configuration leading to a corresponding change in the dislocation glide MEP without affecting the height of the energy barrier. We anticipate that these structural rearrangements will have a significant effect on double-kink nucleation, which must be understood with chemical accuracy to treat urgent challenges in fusion materials science.
We would like to note that in principle it is possible to achieve similar results with a full DFT calculation by deco-rating both dislocation cores with H in a quadrupole cell in a similar way it has been done for carbon impurities in iron [14,15] or oxygen in tungsten [69]. However, in this case it is required to estimate subtle effects of all the uncertainties related to the geometry of the cell on the final results. The problem becomes even more pronounced if we think about extending the study towards larger defects such as impurity clusters and/or vacancy clusters. Our QM/MM approach significantly simplifies the simulation setup.
Finally, we emphasise that the QM/MM framework demonstrated here is able to treat dislocations with significant prismatic character (e.g., edge dislocations) that have much stronger elastic binding with interstitial defects such as H and He but cannot be contained in the small periodic supercells available to fully DFT simulations. Moreover, it is possible to model realistic 3D dislocation kink structures by extending the MM region in the direction of the Burgers vector; however, it requires a reliable empirical force field. Our methodology will thus allow investigations with ab initio accuracy ( 10 meV) into important solute-dislocation interaction phenomena to provide quantitative energetics. Important fusion-relevant examples include the experimentally observed impurity pinning of self-interstitial atom palettes [70] and the novel pipe diffusion behavior in M111 mixed character dislocations [71]. However, we believe that these applications are beyond the scope of the current paper. The atomic configurations created during this research are openly available from the University of Warwick Research Portal at http://wrap.warwick.ac.uk/132290/. ACKNOWLEDGMENTS P.G. and J.R.K. acknowledge funding from the EPSRC under Grants No. EP/P002188/1, No. EP/R012474/1, and No. EP/R043612/1. Additional support was provided by the Leverhulme Trust under Grant RPG-2017-191. We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC under Grant EP/P020194/1. We are grateful for computational support from the UK national high performance computing service, ARCHER, for which access was obtained via the UKCP consortium and funded by EPSRC grant reference EP/P022065/1. Additional computing facilities were provided by the Scientific Computing Research Technology Platform of the University of Warwick.