Grain boundary relaxation and reconstruction: effect on local magnetic moment

We present a detailed numerical study on structure and local magnetic properties of $\langle 100 \rangle$ symmetric tilt grain boundaries in bcc-iron. Particular attention is paid to connection between type of grain boundary relaxation and local magnetic properties. Results from first principles calculation showed that grain boundary reconstruction leads to non-uniform distribution of local magnetic moments in grain boundary plane. This is in contrast with the result obtained in grain boundary plane, where simple relaxation is observed. Well optimized atomic configurations in the vicinity of the interface were achieved by simulated annealing optimization technique improved by combination with genetic algorithm.


Introduction
Grain boundaries (GBs) as interfaces between two grains significantly influence wide range of physical properties of polycrystalline and nano-crystalline materials. From microscopical point of view, the distorted crystal structure along the GBs is mainly responsible for changes in various physical properties of materials [1]. In case of ferromagnetic materials, there is also connection between the distorted structure and local magnetic properties [2,3]. Recently, Ii et al. [4] for the first time made direct measurements of local magnetic moments at grain boundaries in iron. They used electron energy loss spectroscopy on a transmission electron microscope. Since then there have been several attempts to simulate local magnetic moment at iron grain boundaries from first principles [5][6][7][8]. The increase of local magnetic moment on GB plane by 15 to 18% due to magneto-volume effect and oscillatory behaviour in planes parallel to GB explained by Stoner model [9,10] was reported. The results were achieved on symmetric tilt GBs: Σ5(310) [5,7], Σ3(111) [6,8] and Σ5(210) [8]. The experimental study by Ii et al. [4] confirmed the increase of local magnetic moment at grain boundaries. Moreover, the local magnetic moment showed an increasing trend up to misorientation angle 45 • and the cusps on the magnetic moment, i.e., misorientation angle curve at grain boundaries with low coincidence site density, specifically Σ9.
The inevitable part with non-negligible effect on the result is GB type and structure. Consequently, structure optimization and detailed description of optimized GB structure is an essential part of the study of GB properties. Optimization process invokes either relaxation or reconstruction of the ideal geometrical boundary. Relaxation and reconstruction are terms especially connected with free surfaces [11,12]. However, the use of these terms in case of GBs is also reasonable [13][14][15][16]. The most common relaxation observed is oscillatory behaviour of inter-planar distances in the direction perpendicular to GB [7,8,17]. The oscillatory behaviour can be combined with small atomic displacements within planes parallel to GB plane. Whenever the structure of grain boundary plane differs from the structure of the corresponding lattice planes in the grain, one speaks of a reconstructed GB structure. Reconstruction is a result of complicated series of atomic displacements leading to the change of period or coordination number [16], migration of GB plane [15], formation of vacancies or anti-site defects [13], etc.
Due to the high computational demands, first principles calculations use more simple gradient optimization methods. While the relaxation is easily achieved by "simple" optimization techniques, the process seems to be much more complicated in the case of reconstruction. Recently, the application of genetic algorithms (GA) to search stable atomic structures [18][19][20] has identified GA as very effective structure optimization technique.
The aim of our work is to investigate local magnetic behaviour of four bcc-iron 〈100〉 symmetric tilt GBs: Σ5(210), Σ5(310), Σ17(410) and Σ13(510). The first two as representatives with a high coincidence site density and the last two as representatives with a low coincidence site density. The investigated GBs cover the range of misorientation angle from 22.6 • to 53.1 • . To ensure well-optimized structures of the investigated GBs, we improved well-known simulated annealing (SA) [21] technique by addition of some GA features. Inter-atomic interaction is in the process of optimization described by Embeddedatom method (EAM) [22] and final optimized structures act as an input to first principles calculation. The article is arranged as follows: GB supercell construction, new optimization method and first principles simulation details are described in section 2. In section 3 the main results are presented and discussed. Finally, the main contribution of our work is summarized in section 4.

Grain boundary model
Grain boundary supercells were constructed by rotating the two grains in the opposite direction by the same angle around the [001] rotation axis and then matching the two grains together. The result of this is that in the direction of boundary plane the GB shows periodic structure. Atomic positions in the computational supercell are generated according to coincidence site lattice theory. This geometry in combination with periodical boundary conditions creates computational cell which contains two GBs. The first of them is positioned in the middle of the computational cell (GB1) while the second (GB2) is from geometrical reasons split into two parts positioned on the top and bottom edge of the cell. Model of computational supercells is shown in figure 1. GB structure was optimized by our improved optimization technique, which we will refer to as OPSA (optimization by parallel simulated annealing), in combination with EAM. First principles calculation was applied only for magnetic moment computation. Note, two different types of supercells were used for GB structure analysis and for magnetic moment calculations. Due to high computational demands of first principles simulations, the dimensions of supercells designed

43603-2
for magnetic calculation were cropped to provide maximum 100 atoms per supercell. Supercell details are summarized in table 1. Original and cropped GBs were optimized separately, so in both cases both grain boundaries GB1 and GB2 were optimized and, therefore, are equivalent.

Optimization technique
Structures of GBs were optimized by our newly improved and parallelized SA optimization algorithm OPSA. The basic strategy of optimization is based on the technique of SA. It was proved that by carefully controlling the rate of cooling the temperature, SA can find the global optimum. However, this requires infinitely long computational time. In real simulation, the result strongly depends on experimental details, like random generator seed. In the case of GBs, which are in discussion, there are many metastable configurations separated by relatively high energetic barriers. The existence of barriers makes the problem complicated and it is almost impossible to solve it in real time.
To avoid this problem, we applied the technique of parallelization which is based on the idea of many simultaneous runs of identical problem on parallel CPUs. Parallelization was included via certain GA features. OPSA works with a group of individuals (supercells with different GB structure) which are individually annealed in several generations (annealing cascades) ending with crossover operation. The whole optimization process was carried out in detail as follows: At the first step, we generate the geometric position of atoms positioned in the computing cell. The process by which the positions of atoms were generated is described in section 2.1. The data are multiplied N -times where N is the number of processors on which the parallel job is running. Subsequently, atoms positioned in the vicinity of GBs are slightly shifted using the random numbers generator so that displacement of individual atom does not exceed 0.5 Å in each direction. Next step is SA which runs parallel on CPUs for different individuals. The process runs from 350 K temperature to 7.6 K using a stepwise exponential decrease of temperature involving a total of 250 steps. At the end, temperature of 0 K is reached by an acceptance of the position of atoms with lower energy. This step corresponds to zero generation.
At the second step, we create a next generation. The next generation is created by the SA process which runs from temperature 100 K to temperature 7.7 K in 100 steps. The process is finalized with 0 K temperature like in the previous step. The result is a new generation of N individuals.
The third step was motivated from the experiences stating that in the case of complex interface, a group of atoms can occur into local minima that cannot be overcome by the method of SA. The group of these atoms increases the overall energy of GB. To solve this problem, we proceed as follows: after the end of SA in the second step, the energy of each GB is determined and compared with the previous GB energy. If the energy of GB after annealing does not reduce by more than a threshold value (in our case it is 0.005 J m −2 ), this individual is labelled as stuck, and the group of atoms that has caused it is identified as the critical group. Subsequently there is selected the same group of atoms consisting of a critical group located in a successful individual. Successful individual is the individual, whose critical group of atoms has the lowest overall energy. This group of atoms replaces the stuck group of atoms in an unsuccessful individual. Identification is based on the position of the central atom of the critical group of atoms. This is the procedure of mating which is a part of GA. Having carried out this, the algorithm loops back to the second step and repeats this procedure until the desired result is reached.

43603-3
The energy of GB, which is a very important parameter for the assessment of the optimization, is defined as follows: where E GB is the energy of a supercell with two optimized GBs, N is the number of atoms in a supercell, E C is the cohesive energy and s x s y is GB area. The energies were calculated using EAM potential with parametrization provided by Mendelev et al. [23] According to this, parametrization was a lattice parameter for bcc-iron set up to 2.855 Å.
The described method was tested on four GBs, which are considered in our study. Details about dimensions, geometry and the number of atoms in supercells which were used in this calculation are summarized in table 1. We refer to these supercells as SMALL. As a comparison, we performed the same simulation on extended supercells with four times more atoms in GB plane. We will refer to these supercells as LARGE. The energies of GBs obtained by standard SA and OPSA algorithm are summarized in table 2. Note that standard SA corresponds to the first step in the OPSA algorithm. The number of individuals was set to 15 for SMALL and to 23 for LARGE supercell. The number of generations was limited to 50. The results show that OPSA algorithm yields better results for all cases compared to SA algorithm. However, in the case of LARGE supercell, the efficiency of OPSA algorithm increases. This can be explained by a higher incidence of local minima which incurred as a result of the extension of GB. We must emphasize that the result for SA is the best one out of 15 (SMALL), 23 (LARGE) individuals. The results in table 2 also show that while relaxation is easily achieved by standard optimization techniques like SA, reconstruction is not. LARGE GBs Σ17(410) and Σ13(510) remained after SA optimization in an unreconstructed form with unacceptably high energy around 2 J m −2 . The reconstruction was not fully accomplished also for SMALL GB Σ13(510). Therefore, the optimization technique is an essential part of GB study with a non-negligible impact on the result.

First principles calculations
First principles calculations were carried out using the density functional theory (DFT) formalism as implemented in the ABINIT code [24]. The electron-ion interaction was described by the norm-conserving Troullier-Martins pseudopotential [25]. The exchange correlation energy was treated in the local density approximation (LDA) using Teter-Pade parametrization [26]. Before magnetic moment computation, a series of convergence tests were performed and basic bcc-iron parameters were computed. Results are listed in table 3. For comparison, results obtained by our EAM simulation and experimental values are listed in table 3. K -point space was sampled by 2 × 2 × 2 Monkhorst-Pack set, which corresponds to 4 kpoints. The effect of k-point sampling on the quality of local magnetic moment was intensively tested on GB Σ5(210). Results which are shown in figure 2 (a) demonstrate that the number of 4 k-points used provides a sufficient quality and does not need to be increased. The cut-off energy of the plane wave set was set to 1632 eV. Thermal broadening was defined as Fermi-Dirac smearing with temperature of 0.27 eV. It should be noted that the lattice parameter was changed from EAM equilibrium value 2.855 Å to DFT equilibrium value 2.823 Å and no additional relaxation was applied within first principles calculation. The effect of additional first principles optimization was tested on GB Σ5(210), as the most unstable one, concerning the highest GB energy connected with this GB. The comparison of local magnetic moment behaviour without and with additional BFGS quasi-Newton method optimization performed within first principles calculation is in figure 2 (b). The result in figure 2 (b) shows a very good coincidence in the qualitative magnetic moment behaviour and only a small (maximum 6% related to bulk value) quantitative change. Another argument supporting good performance of combination of EAM + OPSA optimization and first principles magnetic moment calculation is a very good agreement between our calculations and local magnetic behaviour on GB Σ5(310) presented by Čák et al. [7]. Comparison is in figure 2 (c). The advantage of our approach lies in a very good performance to computational demand ratio. While first principles optimization takes several days, OPSA + EAM optimization can be achieved within one day.
K -point sampling has a significant impact on the results obtained by DFT calculation. In order to compute GB energy by DFT, we prepared a special supercell. The supercell contains 100 atoms arranged into ideal bcc-crystal. The dimensions of the supercell are 2.823 Å × 14.1150 Å × 28.23 Å, which ensures the shape consistent with GB-supercells described in section 2.1. K -point mesh for a special supercell was adopted from simulation with GB-supercell. In this instance, GB energy is computed as follows: where E GB is energy of supercell containing two GBs, E 100 is energy of special 100-atom supercell, N is the number of atoms within GB-supercell and s x s y is GB area.

Results and discussion
The energy of individual GBs in the process of optimization can be reduced in several ways. In particular, by the change of inter-planar distances between planes parallel to GB, by a rigid shift of one grain with respect to another, or by small atomic displacements within planes parallel to GB. In this case, we define the result of optimization process as relaxation, where various atomic positions are separated only by a low energetic barrier and can be described in terms of the changed inter-atomic distances or bond angles. The result is mostly independent of the choice of the optimization method. The relaxation mechanism was observed in the case of GBs Σ5(210) and Σ5(310). In both cases there was identified a shift of inter-planar distances between the planes parallel to the interface as the main process of relaxation.
Moreover, for the Σ5(210), the mechanism was complemented by a rigid shift of one grain with respect to another in the [001] direction. Note that the shift was 16.8% of the lattice parameter.
The meaning of the optimization process is manifested in the reconstruction of the interface where individual positions of atoms are separated by a high energy barrier, and the efficiency of optimization method plays an important role in GB structure prediction. This is the case of the remaining two interfaces Σ17(410) and Σ13(510). The first feature which indicates the interface reconstruction is the comparison of LARGE GB energies obtained by SA and OPSA techniques (see table 2). Application of an improved structure optimization technique OPSA leads to a significant decrease of GB energy. In particular, by 0.689 J m −2 for GB Σ17(410) and 0.831 J m −2 for GB Σ13(510).
The second feature is the changes in the crystalline structure of the interface plane after optimization. Figure 3 shows the position of atoms in the plane of GB before/after relaxation as a two dimensional (2D) crystalline structure of GB planes of all investigated GBs. It could be seen that 2D structure of GB planes Σ5(210) and Σ5(310) remained unchanged while the structure of GB planes Σ17(410) and Σ13(510) has been significantly changed. GB plane Σ5(210) is characterized by a rectangular Bravais lattice while Σ5(310) GB plane is characterized by oblique Bravais lattice. Note that both lattices contain one Fe atom in the basis. Lattice of Σ17(410) GB plane transfers from rectangular Bravais lattice with one Fe atom in the basis to oblique Bravais lattice with two Fe atoms in the basis. Lattice of Σ13(510) GB plane remained unchanged (oblique Bravais lattice). However, its basis changed from one Fe atom to two Fe atoms. This change could be explained by inward-relaxation of planes adjacent to GB plane, but a better performance of an improved optimization technique OPSA indicates much more complicated optimization process. It should be emphasised that for all the investigated GBs, the mirror symmetry is kept, which is not straightforward from figure 3.
The third and decisive feature was discovered after a detailed analysis of the optimization process. An increased atomic density at GB planes Σ17(410) and Σ13(510) was accomplished by concurrent migration of individual atoms between planes (410) resp. (510) and (001). Since this process cannot be described in terms of the changed inter-atomic distances or bond angles, it is defined as reconstruction. In case the 43603-6 0 5 10 15 x(Å)   In order to compare our data with experimental data published by Ii et al. [4], we computed the average local magnetic moment in the investigated GB planes. The local magnetic moment was averaged over all atoms positioned in GB plane. The dependence of the average local magnetic moment in GB plane as a function of misorientation angle for both simulated and experimental [4] data, is in figure 5. Note that experimental data correspond to various types of GBs while in our case there were considered only 〈100〉 symmetric tilt GBs. The averaged local magnetic moment at GB plane increases up to 2.47 µ B and shows a maximum at misorientation angle 36.87 • . The increasing tendency is in good agreement with experimentally observed data. We have demonstrated that the enhancement of the total magnetic moment on GB plane decreases not only as a result of the volume effect but also due to non-uniform distribution of local magnetic moment on GB plane. Non-uniform distribution of local magnetic moment was identified on low coincidence site density GBs Σ17(410) and Σ13(510) which undergo reconstruction. The same effect might be also responsible for low magnetic moment enhancement measured at GB Σ9, which belongs to high Σ GBs as GBs Σ17(410) and Σ13(510).

Conclusions
It was shown that a prominent feature of certain GBs is the ability to reconstruct the structure during optimization process which in fact reduces the final energy. The effect of reconstruction was observed on low coincidence site density GBs Σ17(410) and Σ13(510) and was accompanied by a significant reduction in energy. Reconstruction was allowed by migration of atoms from the plane (001) to (002). The migration of atoms consequently invokes an increase of atomic density in GB interface and the change of 2D crystalline structure at interface.
Response of local magnetic moment on the structure shows a uniform distribution in the relaxed GB planes and non-uniform distribution in reconstructed GB planes. The non-uniform distribution might be the reason for cusps on the magnetic moment -misorientation angle curve measured by Ii et al. [4]. The same phenomenon in addition to magneto-volume effect could be also the reason for a lower local magnetic moment in low-angle GBs.
We also proposed and applied an improved SA optimization technique, which in fact is a combination of parallel SA and GA algorithms. The main advantage of this technique is its ability to overcome the energy barriers that occur in the process of structure optimization. It was shown that GBs Σ17(410) and Σ13(510) are good candidates to test optimization algorithm qualities, because standard techniques like 43603-8 conventional SA algorithm are not capable of overcoming the energy barriers related to the process of reconstruction. The second advantage is low dependence of the result of optimization process on initial adjustment of random number generator.