Ab initio modelling of intergranular fracture of nickel containing phosphorus: Interfacial excess properties

In the present work, the impact of phosphorus impurities on the grain boundary strength of nickel has been investigated by means of density functional theory (DFT) modelling. Owing to different outcomes and trends previously reported in the literature, it is unclear whether P is strengthening or weakening the Ni grain boundary. To address this issue, we utilize three different DFT based methods: the excess-energy approach, rigid grain separation, and Rice–Wang’s thermodynamic approach. The results show that the commonly used rigid model predicts P to have an increasing effect on the peak stress of Ni of up to 14%, as opposed to a reduction, which is indicated by the excess-energy approach. Employment of the Rice–Wang approach, on the other hand, displays a slight reduction in work of separation. The results show that the discrepancies between previous works can be attributed not so much to the physics of the system, but to the applied model, the partition scheme and the interpretation of the outcomes. This underlines the importance of a proper description of the fracture process, and shows that common simplifications can have a decisive impact on the observed trends.


Introduction
Nickel and its alloys are used in many industrial applications, such as aircraft-and steam turbines, heating equipment, medical applications, thermocouples, and nuclear power systems. The alloys are particularly useful as the properties can be tailored and have good hightemperature characteristics and corrosion resistance. However, if the grain boundaries (GBs) are enriched with impurities, the material may become brittle and susceptible to GB cracking under tensile loading. This is of particular importance in irradiated environments, e.g. nuclear power plants, where Ni alloys are commonly used as structural materials. The radiation field in the reactor core will generate high concentrations of point-defects such as vacancies and self-interstitial atoms in the surrounding materials. An uneven distribution of defects will in turn induce a flux of atoms in the material, which may lead to very high local concentrations of impurities, in spite of insignificant bulk concentrations [1,2]. This in turn can make Ni-based alloys susceptible to creep, swelling, and embrittlement [3,4].
In a recent study, segregation tendencies of Cr, Fe, P, Si, Mn, and Ti in irradiated dilute Ni binary alloys were evaluated in a self-consistent mean field framework [5]. Results showed that Cr, Si and P enrich at sinks following irradiation, whereas Fe, Mn, and Ti are depleted. Within the context of irradiated Ni alloys, it is thus of interest to evaluate the impact of Cr, Si and P on the GB strength. Very few experimental studies regarding the impact of impurities on Ni GBs are available in the literature [6][7][8], instead it is common to evaluate the cohesive properties of interfaces by means first-principles methods. To this end, the Rice-Wang model [9,10] can be used to determine the relative impact of a species on the GB cohesion. The model evaluates the change in work of fracture for GBs containing different amounts of impurities, and based on the outcome it is possible to assess whether an element contributes to either GB strengthening or weakening. However, the model does not give any indication about properties such as critical separation or peak stress associated with cleavage. This information is instead given by energy-and traction-separation curves, which can also be acquired from first-principles methods. But because atomic relaxations during first-principles tensile tests induce an extra level of complexity, such tests are often performed within the approximation of rigid fracture. This approach, known as the rigid grain shift (RGS), neglects the effects of atomic relaxations, which are important for a correct description of surface energies and work of fracture [11][12][13]. However, including atomic relaxations can lead to ambiguity since the peak stress and critical separation depend on the number of atomic layers that are allowed to relax [14]. Specifically, Nguyen and Ortiz derived that the peak stress is inversely proportional to the square root https://doi.org/10.1016/j.nme.2021.101055 Received 5 February 2021; Received in revised form 25 July 2021; Accepted 3 August 2021 of the number of free layers, whereas the critical separation is directly proportional the square root of the number of free layers [15]. Thus, ab initio derived traction-separation data are generally non-unique and appropriate scaling is necessary to ensure accurate interpretation of the results.
To circumvent the issue of non-uniqueness associated with relaxed ab initio traction-separation modelling, Van der Ven and Ceder introduced the notion of excess quantities for cleavage simulations [11]. Within the excess approach, the material comprises interacting atomic layers where interfaces (e.g. GBs) are considered to be localized deviations from the surrounding bulk. Thus, by describing the interaction between the atomic layers by means of parametrized interplanar potentials, the impact of the deviating layer on the mechanical properties can be uniquely determined. Such a description makes it possible to allow for atomic relaxations, while the results are independent of the number of atomic layers in the system. An additional attractive benefit of the excess approach is that together with bulk potentials, the excess properties can be used to generate atomistically informed multi-layer potentials that relates the traction to the separation on a larger scale. Consequently, the relaxed excess properties can be used to construct analytical coarse-grained cohesive zone models of arbitrary thickness, which can be adopted for upscaled continuum failure modelling. Thus, the approach enables a physically accurate bridging between the atomic scale and the larger scale models, which provides a sound basis for a multi-scale description of the fracture response of a brittle material.
GB embrittlement of Ni due Si, Cr and P has been the focus of several theoretical ab initio studies [6,[16][17][18][19][20][21][22][23]. The reported results generally indicate that Si and Cr act as GB cohesion enhancers in Ni [6,18,19], whereas the impact of segregated P impurities in Ni is less clear. Previous studies give contradictory results, and the element has been predicted to either strengthen or embrittle the material in different studies [6,20,21,[24][25][26]. This discrepancy has previously been suggested to be due to a concentration effect, where P at low concentrations increase Ni GB cohesion, whereas at higher concentrations, repulsion between P atoms weakens the GBs and leads to embrittlement [23,24]. In light of the contradictory results in the literature regarding the role of P impurities on the GB strength of Ni, and because P in GBs is notorious for weakening other metal materials (e.g Fe [27][28][29][30], W [31][32][33][34][35], and common austenitic and ferritic alloys [36][37][38][39][40]), further investigation of its impact on the GB strength of Ni is needed. Moreover, the conjecture that the underlying mechanisms for the contradictory behaviour is due to concentration effects was made without detailed consideration of how system parameters and relaxation schemes can impact the results. Indeed, it has been found that the results in the literature indicating that P has a negligible embrittling effect on Ni GBs, have mainly been produced using the Rice-Wang framework [6,20,25], while studies involving the RGS approach indicate that interstitial P can provide a significant increase of the peak stress for the same GB [24,26]. Although impurity-induced changes in peak stress and work of separation are intimately connected and often exhibit similar trends, ultimately the impact on the interfacial work of separation, as evaluated in the Rice-Wang model, may differ from the impact on peak stress, which is assessed in the ab initio stress-test. Thus, the interpretation of the results and the simplifications associated with the two frameworks may produce different conclusions regarding the embrittling potency of impurities, which may be an explanation of the discrepancies found in the literature. A detailed investigation of such effects is for this reason of great interest to the community, in addition to providing new insights on the impact of P in this context.
The goal of the current work is for this reason twofold; the first and most important is to provide a full description of the impact of P on the Ni traction-separation behaviour, and to quantify the deviations associated with the more simplistic models commonly used in the literature. To elucidate the impact of the solute, the stress-response of a Ni GB is evaluated using the relaxed interfacial excess approach for different P impurity contents. The impact of common approximations has also been evaluated by performing the corresponding analysis in both the RGS and the Rice-Wang frameworks and we have investigated the role of different partition schemes. The second objective is to derive a full set of excess energy parameters associated with bulk and the decohering region, to enable the formulation of large-scale models to evaluate brittle fracture of the material. The results of the current study can in this way provide important insight, not only on Ni cohesive properties, but also on how common approximations can impact the results in this context.

Geometry and preferential segregation
In the present work, decohesion was modelled within the framework of density functional theory (DFT). The considered GB geometry is the symmetric tilt Σ5 [100] (012) GB, which was generated using the GB studio software [41] and is illustrated in Fig. 1. The choice of studying the stress response of this particular GB is motivated by the fact that it has previously been modelled in the literature in both the RGS and the Rice-Wang frameworks with different reported outcomes. The clean (impurity free) GB structure was first optimized using full coordinate relaxation along with supercell relaxation in the GB normal direction (the [012]-axis in Fig. 1). For the systems containing impurities, the structure was again optimized along the GB normal direction following the addition of the P atoms. Decohesion simulations were performed by straining the supercell in the GB normal direction whereby the system was allowed to undergo ionic relaxation at constant volume. To reduce interactions between mirror images and minimize the impact of finite-size effects, large bulk and GB supercells were used, both with dimensions 7.0 × 7.9 × 31.5 Å 2 , and comprising a total of 160 Ni atoms. Four possible solute sites on one of the two GBs in the supercell were considered when investigating the impact of P on the GB properties. Three of them (denoted 1 − 3 in Fig. 1(a)) are substitutional configurations in direct vicinity of the GB. The final position is an interstitial site, which is referred to as in Fig. 1(a). The four interstitial sites illustrated in Fig. 1 To evaluate the impact of P on the cohesive properties of Ni, the relative stability of the impurities at various GB sites was first assessed. This was done based on segregation energies according to Eq. (1). Seg = (GB + ) + (Bulk) − ( (GB + ( − 1)) + (Bulk + 1)) where E(GB+ ) is the energy of the GB containing impurities, in which the subscript corresponds to or , representing the substitutional and interstitial configurations, respectively, (Bulk) is the energy of pure bulk, and (Bulk+1) is the energy of the bulk containing a substitutional/interstitial solute. Eq. (1) describes the energies associated with an incremental increase in defect concentration in the GB and, according to this definition, negative energies represent preferential segregation of the solute to the GB.
Since the segregation energy evaluates the energy of the solute in the GB, in reference to its energy in the bulk, the choice of an appropriate bulk reference is of importance. Under equilibrium conditions, P is more stable as a substitutional atom in the bulk [5], and thus segregation energies should under such prevailing conditions be used with this configuration as reference. However, in irradiated environments, the concentration of interstitial atoms (self or foreign) is significant. In such conditions, P can migrate independently of vacancies or other defects [5]. As the interstitial arrives at the GB, if the latter is previously enriched in vacancies, the solute will there take on a substitutional configuration. In the absence of vacancies in the GB, P will there take on an interstitial configuration instead. Since the P atom can originate from either interstitial or substitutional bulk sites in an irradiated environment, both configurations have here been used as reference when calculating the segregation energies of P. As the presence of vacancies is required for a migrating interstitial to take on a substitutional configuration in the GB, calculations were in this case performed under the assumption that the GB contained a vacancy prior to the arrival of the solute.

Cleavage fracture modelling
The response of a material interface subjected to external stresses can be estimated by a cohesive zone (CZ) model, which describes the emerging traction as the material is strained. The purpose of such a model is to predict the behaviour in the decohering region, which for maximum accuracy should be separated from the elastic response of the underlying atomic planes in the bulk. The distinction between the two regions can be made by determining the system's excess quantities, defined as the difference in properties of the fracturing solid and those of the ideal bulk when subjected to the same local stress state [11].
By means of first-principles methods, the energy of locally stable configurations can be evaluated as a material is deformed in controlled steps. Under the assumption that the atomic planes of the crack remain parallel during decohesion, first-principles methods can provide a good basis for assessing the stress-response of a decohering system. This is done by fitting the strain energy versus separation from ab initio calculations to a model curve, such as e.g. the universal binding-energy relation (UBER) derived by Rose et al. [42,43]. From the parameterization of the UBER curve, the potential energy of the interactions between atomic planes is characterized. This gives the traction-separation response of the system, as the first derivative of the energy with respect to the separation. In the literature, fracture processes are often modelled from first-principles by the RGS approach, in which the distance between two rigid grains is increased in discrete steps at a predefined fracture plane. However, to improve the physical basis for the excess properties, atomic relaxations should be incorporated in the model. If such relaxations are included, it has been found that the original UBER may not be sufficient to describe the interactions between atomic planes, but an extension of UBER, the xUBER, given by Eq. (2), has been found to yield an improved description of the excess potential [35,[42][43][44].
where is the separation, is the stiffness and the parameters , and are to be fitted. When describing the fracture process as the decohesion of a cohesive zone embedded in an elastic medium, all deviations from bulk behaviour are incorporated in the excess properties. Results obtained in this framework (henceforth referred to as the excess energy assessment (EEA) method) are independent of the number of atomic planes included in the model [11,44,45]. For this reason, the EEA method enables an unambiguous evaluation of the interplanar potential. However, the approach makes parameterization of energy curves significantly less straightforward compared to rigid calculations. Unlike in the case of the RGS, the ab initio results can in the case of the EEA method not be directly fitted to the UBER model curve. Instead, in the EEA method, the energy obtained through ab initio calculations represent the combined response of the interfacial strength of the material, and the elastic response of the underlying atomic planes. Thus, to parameterize the xUBER curve within the EEA framework, the response of the material must be separated into the two respective contributions. This can be done by describing each separate binding interaction between two atomic planes as a nonlinear xUBER potential [15].
Following the procedure outlined in [35], the total energy, tot , obtained by ab initio calculations, comprises the sum of all individual contributions, and is in the case of a homogeneous bulk material given by Eq. (3) [15].
where Bulk is the energy of an individual xUBER potential, is the number of atomic layers in the model that are allowed to relax, and the individual elongations, 1 and 2 , correspond to the interplanar stretch associated with the surrounding bulk and the decohering interface, respectively.
The objective of the EEA method is to parameterize Bulk∕GB ( 2 ), which is the excess energy that characterizes the stress-response of the decohering region. In the case of bulk material, the excess energy can be extracted by fitting energy versus separation curves obtained from first-principles calculations to Eq. (3), where the energies of the individual UBER potentials are on the form of Eq. (2). Additionally, the condition that the total elongation, tot , consists of the individual elongations must be fulfilled, i.e. tot = 1 + 2 (4) and the individual potentials should adhere to the requirement of mechanical equilibrium of global and local stresses at the CZ/bulk boundary [35], i.e.
For the case of GB decohesion [35], the total energy of the system is given by Eq. (6).
In this case, the elastic response of the surrounding medium, Bulk ( ), should follow bulk behaviour as determined from Eq. (3). Thus, in order to extract the excess properties of a GB, the xUBER analysis of the bulk material must first be performed. Once Bulk ( ) has been extracted, the excess behaviour of the GB can be obtained from a similar approach as for the bulk; the total energy relation is given by Eq. (6), the individual elongations follow the relation of Eq. (4), and the requirement of mechanical equilibrium, Eq. (7), must be fulfilled.
For a given and tot , the relation between 1 and 2 is given by the set { 1 , 2 } that minimizes the total energy and satisfies mechanical equilibrium. For GBs, this relation is always given by Eq. (4), whereas for bulk material and small elongations, the strain will be uniformly distributed in the material, and 1 = 2 . The elongation at which damage is initiated is known as the critical separation, c , and the stress at this point is known as the critical stress, c . When the crystal is strained beyond the critical separation, one of the interlayer spacings is assumed to undergo maximum separation, while the others remain uniformly elongated, such that 1 > 2 and Eq. (4) is satisfied for both bulk and GB [35,45]. If additional stress is applied, the crystal will eventually crack at 2 , and the 1 elongations will return to their equilibrium distance.
By using the EEA approach, both the bulk and GB excess energies can be extracted. The GB excess energy marks the behaviour of the decohering GB, isolated from the elastic response of the bulk. The method provides a framework to assess how the stress-response of a material is affected by the presence of a GB. Additionally, by performing the EEA analysis on a GB containing impurities, it is possible to evaluate the relative impact of a solute on the GB cohesive properties.

Work of separation and embrittling potency
In both the RGS and the EEA methods, the energy curves converge towards the work of separation (WOS) in the asymptotic limit, given by Eqs. (8) and (9), for the bulk and GB, respectively [46,47].
where is the surface energy of the formed surface and WOS GB = S1 + S2 − GB (9) where GB is the GB energy, S1 , and S2 are the energies of the two newly created surfaces, respectively. For the case of symmetric tilt GBs with identically terminated surfaces on both sides of the decohering plane, S1 = S2 = S . The embrittling potency was calculated to estimate the impact of P on the GB strength by means of the Rice-Wang approach [10]. This method describes the change in WOS as impurities are added to the GB. In the case of the symmetric tilt GB, the embrittling potency is given by Eq. (10) [48].

Modelling details
The traction-separation was modelled by both the RGS and the EEA methods. In the EEA analysis, two initial configurations were evaluated from which the system subsequently was allowed to relax; one in which the crystal was homogeneously stretched, and one in which the entire elongation was localized in 2 [11]. For all elongations, the configuration corresponding to the lowest energy, and thus the global minimum, was used for further analysis. As periodic boundary conditions were implemented in this work, the supercells contained two GBs/free surfaces (FSs). For this reason, some atom layers were kept fixed at their equilibrium distance, so that the entire supercell was not subjected to strain, and only one decohering interface underwent deformation. We converged the WOS with respect to the number of atomic planes free to relax for 4 ≤ ≤ 10 and found that = 6 was sufficient to yield well-converged results, thus allowing seven interatomic spaces to be strained. The calculations were performed under the assumption that the interplanar potential was truncated beyond the nearest neighbours.
All DFT calculations were conducted using the Vienna ab initio Simulation Package (VASP) [49,50]. The atomic relaxations were performed using the conjugate gradient algorithm, with a maximum force Table 1 Segregation energies and excess volume for P in Ni, as function of impurity coverage, . and correspond to the interstitial and substitutional sites in Fig. 1 convergence criterion of 10 −2 eV/Å. Standard valence pseudopotentials from the VASP library generated with the projector augmented wave (PAW) method using the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional were used in all calculations [51]. To ensure convergence for the plane wave representation, a Methfessel-Paxton smearing with a broadening of 0.2 eV was used. All calculations were spin-polarized, the energy cut-off was set to 400 eV, and the Brillouin zone was sampled with 3 × 3 × 1 k-points using the Monkhorst-Pack scheme. The equilibrium lattice parameter was set to 3.522 Å following the convergence study of our previous work [5].

Preferential impurity segregation sites and cleavage planes
To investigate the impact of P on the Ni stress-response, the preferential configurations and relative stabilities of the solute in the material should first be determined. This was done based on segregation energies, Eq. (1), which describes the energy change associated with taking a solute from the bulk and place it in the GB. As illustrated in Fig. 1, the considered GB had room for up to four interstitial P atoms. With the exception for I = 2, the placement of the P atoms gave symmetrically equivalent configurations. However, when two interstitials were present in the system, two possible configurations emerged. They correspond to the second interstitial placed either at [010] or (1/2)[021], relative to the position of the first, see Fig. 1(b), such that their separations correspond to ideal second (2nn) and third nearest neighbour (3nn) distances in the fcc structure, respectively. Ultimately the latter configuration was found to be more stable (energy difference 0.3 eV), and was used in the further calculations. The interstitial P atoms were added sequentially to the GB following the order of Fig. 1 As discussed in Section 2.1, we consider both interstitial and substitutional P as reference for the segregation calculations to account for irradiated and unirradiated conditions. As can be seen in Table 1, the interaction between P and the GB is stronger in reference to the interstitial configuration. This was expected, since the solute is more stable as a substitutional compared to an interstitial in the bulk. More important effects, as for instance a shift in stability of the considered configurations between the two environments, were not seen. This is an indication that P can be found as both an interstitial and a substitutional in the GB in both irradiated and equilibrium conditions. For this reason, no further distinction between the impact of the solute in the two environments was considered beyond this point. It should however be noted that the concentration of interstitials is orders of magnitude lower in equilibrium conditions, for which reason the solute will mainly take on the substitutional configuration in such environments. Results in Table 1 show that both the interstitial and substitutional (positions E. Toijer et al. configurations are either unstable or considerably less favourable compared to position 2 . The two positions are for this reason omitted from further consideration. Regarding the interstitial configuration, results in Table 1 indicate an increasing stability with increasing number of P atoms in the GB up to the highest filling degree. Thus, based on the results in Table 1, we investigate the impact P segregation at the 2 site and up to complete filling of the interstitial sites. When multiple interstitial impurities are present in the system, it is of importance to ensure that the traction-separation behaviour is evaluated along the minimum energy path. Thus, the cleavage plane termination that yields the lowest WOS should be identified. In the present case, a number of different partitioning schemes were considered; the P atoms could end up in either collected configurations, with all P atoms on the same FS, or in a dispersed configuration, in which the P atoms are distributed over the two FSs. The WOS of a large number of partitioning schemes were for this reason evaluated, and a selection of the obtained results are presented in Table 2. Both rigid and relaxed calculations have been used to identify the corresponding minimum energy fracture path. Additionally, since traction-separation curves are often smoother and more easy to fit when relaxations are considered in the GB normal direction only (1D relaxations), compared to when relaxations are allowed in all dimensions (3D relaxations), both cases were considered. This allowed for a proper evaluation of how the number of relaxation degrees of freedom could impact the WOS.
The results in Table 2 show that, if atomic relaxations are included in the calculations, it is possible for the system to lower its energy by moving atoms to a more favourable configuration. A consequence of this was that the relative stability of the considered partitioning schemes varied for the different relaxation schemes. The configurations representing the minimum energy partitioning schemes for the different relaxation schemes are presented in Table 3. Since the lowest energy configuration should be used to evaluate the stress-response of a system, the RGS and EEA analyses in the current work were performed along their corresponding minimum energy path. As a consequence, identical terminations were not used in the RGS and the EEA analyses. This issue will be further discussed in Sections 3.4 and 4. In the case of interstitial P, the results in Tables 2 and 3 indicate that 1D relaxations were not sufficient to capture the minimum energy configurations, for which reason 3D relaxations were used in all further calculations related to the EEA model.

Decohesion of bulk Ni and the Ni GB
In Fig. 4, the results from the EEA analysis of the clean GB and bulk Ni are presented, together with total energies obtained from ab initio calculations. The parameters of the fitted excess energy curves obtained in the EEA analysis are presented in Table 4. As can be seen in the figure, the results from DFT are well represented by the model of the current work. This indicates that the energy-separation behaviour of the system is captured by the sum of the individual xUBER potentials, according to Eq. (3).
From the EEA analysis, the excess energy associated with the decohering plane, Bulk∕GB ( 2 ), was obtained. The results are presented in Fig. 5a, together with the corresponding results from the RGS analysis. In Fig. 5b, the traction-separation curves associated with the two methods are presented. The figures indicate that the only important difference between the two models, is the GB WOS, with 3.5 J m −2 obtained from EEA and 3.8 J m −2 from RGS. In the bulk system, however, the difference in this quantity is significantly less important (4.8 J m −2 from EEA, 4.9 J m −2 from RGS). The two methods also display comparable results for the traction-separation response of the material. Indeed, a GB critical stress of 25.0 GPa was obtained by E. Toijer et al.   both methods, whereas the bulk critical stress was calculated to 30.0 GPa and 29.2 in the EEA and the RGS frameworks, respectively. The minor differences indicate that the RGS approach can be considered a reasonable approximation for describing the fracture of bulk Ni, and the clean Ni GB.

Decohesion of the Ni GB containing P
The excess energy and traction-separation data from the EEA and RGS analyses for the Ni GB containing P are presented in Figs. 6 and 7, respectively. The results have been separated into low concentration effects ( , = 1) and higher concentration effects ( = 2-4), for clarity. The corresponding results for bulk Ni and the clean GB are included in both figures for reference. All fitted parameters of the xUBER model curve from the EEA analysis are presented in Table 4.
Figs. 6 and 7 indicate that P only has minor effects on the WOS of Ni, regardless of concentration or applied model. However, the EEA and the RGS methods give significant differences in traction-separation behaviour. In the RGS approach, a clear trend is seen with increasing P concentration. The peak stress increases up to 14% for = 4, a strengthening which corresponds to the peak stress of the perfect bulk. For low concentrations ( ∕ =1), on the other hand, the computed RGS peak stress only deviates marginally from that of a clean GB (difference in peak stress is less than 1%). In the case of the EEA method, the low P concentrations result in a decrease in critical stress, indicating that the P in this case contributes to ∼7% weakening of the material. A weakening effect (∼3%) is also seen for = 4. However, for = 2-3, no significant impact on the GB strength is seen. Thus, a trend corresponding to the findings in the RGS framework, could not be seen when the EEA model was applied. It has previously been shown that the presence of P give rise to magnetically dead layers in the Σ5 (210) Ni GB [7]. This phenomenon could possibly contribute to the observed impact of the solute on the system's stress response. However, a full evaluation of how such effects could impact the peak stress is beyond the scope of this work. The impact of P on the peak stress of the Ni GB is summarized in Fig. 8 for the two methods, respectively. The disparate behaviour displayed in the figure is a clear indication of the importance of inclusion of atomic relaxations in the model. It should be noted that including atomic relaxations in the calculations were shown to impact the relative stability of the considered surface terminations. For this reason, identical terminations were not used in the RGS and the EEA analyses. Consequences of this will be further discussed in Section 3.4.
The effect of P in Ni has also been evaluated based on the embrittling potency, Eq. (10). The results are presented in Table 5, where the corresponding effect of the solute from the EEA and RGS analyses are also presented for comparison. For all cases, > 0, which is an indication that WOS decreases by the introduction of P impurities. This is valid regardless if the EEA or RGS method is employed, or if the impurities occupy interstitial or substitutional sites. However, the reductions can be considered to be very small for both, as 0 < < 0.18 J m −2 and 0 < < 0.25 J m −2 for the EEA and RGS approaches, respectively.

Impact of partitioning scheme on the Ni GB stress-response
As can be seen in Figs. 6-8, the impact of P is significantly different in the EEA and RGS approaches. The fact that different partitioning schemes were used in the two methods may be the possible source of this discrepancy. In order to evaluate if this was indeed the case, the impact of all three = 4 configuration, illustrated in Fig. 3, were assessed in the EEA framework. The resulting excess-energy and traction-separation curves are presented in Fig. 9. The corresponding changes in peak stress are quantified in Fig. 10.
As can be seen in Fig. 9, the stress tolerance depends on the partitioning scheme, and is particularly high when the system is fractured according to the zigzag configuration. Due to the considerable increase in surface energy of this configuration, it is however deemed unlikely to occur. The main interest of the current section is to evaluate whether E. Toijer et al.

Table 4
Parameters of the UBER curve, in Eq. (2), obtained from the EEA analysis. 0 and 1 were set to 1, and 2 was set to 0 to assure consistency with the original UBER format. [Å] C [eV Å −4 ]  the significant difference between the impact of P when evaluated in the RGS and the EEA frameworks respectively (Fig. 8) can be explained by the use of different partitioning schemes. Since the collected configuration was used in the RGS method, the corresponding result obtained with the EEA method should be used to evaluate such effects. As can be seen in Fig. 10, this configuration results in a small increase in peak stress compared to the clean GB. However, the effect is minor compared to the impact of the same configuration when evaluated in the RGS framework. Thus, these results indicate that although a correct choice of partitioning scheme is of considerable importance, the difference between the RGS and the EEA approaches cannot be fully attributed to such effects.

Discussion
Although the stress-response of the clean Ni systems is equivalent when evaluated in the EEA and the RGS frameworks, more significant differences emerge as P is added to the system. Owing to the lack of experimental data regarding the effect of P in Ni, such information cannot be used to estimate the accuracy of the results. However, the current work gives a quantitative description of uncertainties associated with the considered models. The Rice-Wang approach is a computationally efficient method that can be used to estimate the direct effect of solute addition to the GB. However, conclusions in this model are based solely on the sign of the embrittling potency, and the importance of small energy differences risk being exaggerated. This is of particular importance in the Ni-P system, for which it was shown in the current work that the WOS varied very little between the clean and the P covered GBs. Consequently, the small energy variations can explain why previous work in the literature suggested that the impact of P in Ni ranges from cohesion enhancing [20], via no effect, [6] to embrittling [21], although in all these studies, the Rice-Wang method was used.  In order to evaluate the traction-separation behaviour of a decohering system, it is possible to utilize either the RGS or the EEA approach. The two give very similar critical stresses for the clean GB and bulk material, however, Figs. 6 and 7 display significant differences on the impact of P. Indeed, the EEA method indicates that the solute has a weakening effect, although more important at lower concentrations, whereas the RGS method suggests a considerable strengthening effect. Both results are in contrast to the findings in the Rice-Wang framework, where the effect is considered minor to negligible. While results in the literature indicate that P has a negligible effect in Ni, they have mainly been obtained within the Rice-Wang framework [6,20,25]. In studies performed using the RGS approach, interstitial P at high concentration has been shown to significantly strengthen the GB [24,26]. These results are very well in line with the results of the current work, and can be seen as a clear indication that the applied method can significantly impact the observations. These results strongly suggest that, although the RGS method can provide an acceptable description of the tractionseparation response of clean Ni, it cannot fully capture the effect of impurities on the system. One important source of error in this context, is that the minimum energy partitioning scheme could not be found using the rigid approach. However, results in Section 3.4 show that this is not a sufficient explanation as to why such a significant deviation is observed in this framework. Instead, the seemingly disparate results from the reduced WOS and increasing peak stress with increasing impurity content for the RGS approach, can be explained by the fact that even though the peak stress increases, the strain energy, i.e. the overall area under the traction-separation curve, decreases. This is an indication that, despite the peak stress increase, the toughness of the material is generally decreased, such that the amount of energy that can be absorbed by the material before fracture is reduced.
As discussed in the introduction, atomic relaxations are important for a correct description of surface energies, and should be included for a proper description of the system. The results in the current work indeed underline this fact, but also the importance of adopting the lowest energy fracture path. Since the EEA method takes both factors into consideration, in addition to providing a full description of the emerging forces during decohesion, results in this framework are considered more reliable compared to the two others. The fact that P was shown to have a weakening effect further increases the reliability of the results. However, significantly more computational resources and post-processing of data are required in the EEA model, for which reason this approach is perhaps not always a convenient choice. With that said, an important benefit of the EEA model is that it also can be used to obtain a full set of parameters which characterizes the cohesive zone, and can thus be used in larger scale cohesive zone models.

Conclusion
In the present work, the impact of P on the stress response of Ni GBs has been estimated. This has been done by applying three separate ab initio approaches; the Rice-Wang model, the rigid grain shift (RGS)  Fig. 3. Results from obtained in the EEA framework. Bulk Ni and the clean GB are included for reference. Fig. 10. Impact of partitioning scheme (see Fig. 3) in the EEA analysis on peak stress of the Ni GB containing four interstitial P atoms. and the excess energy assessment (EEA). The goal was to improve the current understanding of how P can impact the cohesive properties of Ni, and to investigate how the results can depend on applied model. The Rice-Wang approach is a computationally efficient method, which can be used straight-forwardly to estimate the impact of a solute on the cohesive properties of a material. However, it does not provide any information about the gradual decohesion process or insight on the associated critical stress. This information can instead be obtained in both the RGS and the EEA frameworks. An important difference between the two methods is that the EEA approach allows for atomic relaxations, all the while results are kept independent of the number of free planes in the system. This is not the case for more simplistic models. As a consequence, result obtained in this framework can be used straight forwardly in larger scale models to further investigate the process of brittle fracture of the material.
Results in the current work display a similar stress-response of bulk Ni and the clean Ni GB in the RGS and the EEA frameworks, for which reason the former provides an appropriate approximation for the two systems. However, the relative trends regarding the impact of P were not consistent between the two methods. The RGS model indicated that P has a significant strengthening effect on the material, whereas the EEA method indicated that the solute had a negligible to weakening effect. The EEA approach is in many ways more thorough compared to the two others, and results obtained in this framework are generally considered more reliable. It is for this reason concluded that the RGS approach, which is very common in the literature, is associated with significant uncertainties when evaluating the impact of solutes in the current and similar systems. The method should for this reason not be used without proper consideration of how the many simplifications can impact the results.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The data required to reproduce these findings are available from the corresponding author upon reasonable request.