A new ab initio modeling scheme for ion self-diffusion coefficient applied for {\epsilon}-Cu3Sn phase of Cu-Sn alloy

We present a new modeling scheme for ion self-diffusion coefficient, which broadens the applicable scope of ab initio approach. The essential concepts of the scheme are `domain division' and `coarse graining' of the diffusion network based on the barrier energies predicted by the ab initio calculation. The scheme was applied to evaluate Cu ion self-diffusion coefficient in {\epsilon}-Cu3Sn phase of Cu-Sn alloy, which is a typical system having long-range periodicity. The model constructed with the scheme successfully reproduces the experimental values in a wide temperature range.


INTRODUCTION
Ion diffusion attracts the broad interests of the material researchers because it dominates the various important phenomenas in the solid: corrosion, monotectoid, fracture, degradation and so on. In order to elucidate the microscopic mechanism of the ion diffusion, ab initio calculation is one of the most powerful tools. For example, it can evaluate the barrier energy of the diffusion route with such techniques as nudged elastic band (NEB) method [1]. Even the ion self-diffusion coefficient [2] has been reported theoretically, [3][4][5][6][7][8] combining the ab initio predictions with modeling schemes such as fivefrequency model. [9][10][11] These theoretical works however target just simplest crystals mainly due to the following two difficulties. First, the system appearing in the practical situation often requires a large size of unitcell because of long-range periodicity, which makes the calculation expensive, otherwise infeasible. Such a problem is especially severe, for example, when optimizing the ion path through the diffusion route using NEB method because it requires a lot of force-field calculations for a number of structures. Secondly, it is very hard to count up and evaluate all of the ion diffusion processes in most cases so the application range of ab initio approach is limited within just simplest systems.
To overcome the difficulties, we introduced a couple of novel concepts, which significantly simplifies the diffusion network based on the barrier energies. First, suppose to classify the diffusion routes according to the barrier energies into the three groups I-III as the energies ordering like ∆E I ≪ ∆E II ≪ ∆E III . Since the barrier energy exponentially contributes the diffusion coefficient (i.e. Boltzmann factor), the diffusion routes (III) can be excluded from the diffusion network. Then, the network may be divided into several disjunct domains and the expensive calculation for the large unitcell can be replaced by the several cheaper calculations for the small domains. This is a concept to solve the first difficulty. In addition, the diffusion network in each of the domains may be further simplified by coarse-graining: The vacancy can move almost freely through the diffusion routes (I) compared with moving through the routes (II), so the ion sites connected by the routes (I) can be represented by just a site. Eventually, the diffusion network falls into coarse-grained one with the representative sites, where it is much easier to count up the diffusion processes.
We established a modeling scheme based on the above concepts for the example of the Cu ion self-diffusion coefficient in ε-Cu 3 Sn phase of Cu-Sn alloy. This is a typical system having long-range periodicity and there is abundant experimental data for the Cu self-diffusion coefficient. [12][13][14][15][16][17] From simulation side, classical molecular dynamics (MD) was applied [18] but it overestimated the experimental values by a digit. On the other hand, the model constructed with our scheme successfully reproduced them in a wide temperature range.
This paper consists of the following sections: In the next section (Formalism), we introduce the equations for the physical quantities relevant to the ion diffusion and also define the technical terms and the notation rules. In the 3rd section (Domain division), we explain how to divide the diffusion network into the multiple types of the simple domains, based on the barrier energies. In the 4th section (Details and results of ab initio calculation), we present the predictions for the quantities (ex. barrier energy, vacancy formation energy) calculated with the density functional calculation The details of the calculations are also given here. In the 5th section (Coarse graining), we discuss how to represent the multiple ion sites with just a site and coarse-grain the diffusion network, based on the barrier energies. In the 6th section (Formula of self-diffusion coefficient), We discuss how to model the self-diffusion coefficient when there are two or more types of diffusion routes. In the 7th section (Evaluation of correlation factor), it is deeply discussed how to count up the possible diffusion processes in the simplified diffusion network and evaluate the correlation factor with some approximations. In the 8th section (Results and discussions), the validity of the model constructed with our scheme is discussed, compared with the experimental values for the self-diffusion coefficient or the correlation factor for ideal hexagonal 2-dimensional lattice. In the 9th section (Summary and prospects), we summarize this paper and refer to the advantage of our scheme from ab initio MD.

FORMALISM
The self-diffusion coefficient D(w) represents the macroscopic ion flow in a direction w driven by self-diffusion. D(w) is obtained by evaluating the ion jump, which is a minimum process of the ion diffusion. We index the diffusion routes connecting to a site with j and define θ j (w) as the angle of the ion jump thorough route j from direction w. We also define D * j as the angle-independent contribution of the ion jump to the self-diffusion coefficient (explained in the next paragraph) and the self-diffusion coefficient D(w) is evaluated from the following equation: [19] D * j consists of the four quantities appeared in the following equation (index j is omitted): [19] (2) d is the distance of the diffusion route. C v is the vacancy formation rate given as the Boltzmann factor of the vacancy formation energy ∆E vac . η is the jump rate depending on the barrier energy ∆E barrier and the vibration frequency ν of the attentional ion (tracer) in the jump direction: f is called correlation factor and explained in the next paragraph. The evaluation of the vacancy formation energy ∆E vac needs the chemical potential µ Cu of Cu ion. We calculated it for face-centered Cu mono-crystal and got the ∆E vac from the following equation: Here, E perfect is the total energy of perfect crystal and E vac is that of the crystal including a vacancy defect.
The evaluation of the correlation factor f is the most awful part for the modeling of self-diffusion coefficient. The factor reflects the position exchanges sequentially occurring between the vacancy and the tracer: After the tracer jumping via the vacancy, it is at the behind of the tracer. Thus, the vacancy tends to pull back the tracer at the next ion jump around the reverse direction of the ion jump. In addition, such a pullback can occur again and again after the 1st pull-back, always around the reverse direction of the preceding tracer move. The correlation factor f represents how much the sequential pullbacks affect the self-diffusion coefficient in total. Since the 1st pull-back denying the tracer jump has the largest contribution to the correlation factor, it always reduces D * j (i.e. 0 ≤ f ≤ 1). To evaluate correlation factor, we have to count up all of the vacancy tracks from the one to the next position exchange between the tracer and the vacancy. This is the identity of the first difficulty. The contribution from each of the pull-backs to the self-diffusion coefficient is given as average cosine: Here, P k is the sum of the realization probability of the vacancy tracks causing the pull-back from site k. θ k is the angle of the tracer move at the pull-back from the preceding tracer move. We also define nth-order average cosine cos θ (n) , which represents how much the first ion jump is denied by the n times pull-backs in average. The correlation factor f consists of them as the following equation: '(n+1) times of pull-back' is interpreted as one extra pull-back occurring after 'n times of pull-back' so cos θ (n) is given as cos θ (n+1) = cos θ (n) cos θ .
Especially when there in only one diffusion route, [20] the following relationship is hold: [19] cos θ (n) = cos θ n .
This equation is generalized for multiple diffusion routes later in section (formula of correlation factor).

DOMAIN DIVISION
ε-Cu 3 Sn phase has, it is reported, long-range periodic structure (Fig. 1a). It has been recognized that the structure consists of the unitcells of Cu 3 Ti-type structure with shifted by 1/2 in fractional coordinate along a-axis, every M/2 unitcells lining up along b-axis (Fig. 1a). The various periodicities M are reported experimentally as the even numbers from 2 to 12. [21,22] Based on the ab initio predictions, it is newly found that the ion jump though the diffusion route passing between the Sn sites hardly contributes to the self-diffusion because their barrier energies are much higher than the others. After excluding these routes from the diffusion network, it is divided into the two types of simple domains, Cu 3 Ti-type and D0 19 -type, as shown in Fig. 1b. The both domains exist in the ratio of (M-2):2 so the self-diffusion coefficient of the target phase is given as the weighted average of them with the existing ratio:

DETAILS AND RESULTS OF AB INITIO CALCULATION
We employed 2×2×2 supercell for both Cu 3 Ti-type and D0 19 -type domains to reduce the spurious interaction between the vacancies. The lattice parameters are fixed at the experimental values. [21] We used VASP [24] for the density functional calculation to evaluate the physical quantities. The barrier energies are evaluated with climbing-nudged elastic band (c-NEB) method [1] implemented in VASP. [24] c-NEB needs the interpolated structures between the both edges of the diffusion path as inputs and optimizes them with connected by the virtual springs. c-NEB guarantees that one of the structures is positioned at the exact saddle point of the given potential surface, in which c-NEB is superior to the original NEB method. [25] We prepared 15 interpolated structures and set the spring coefficient as 5 eV/Å 2 . We used PBE functional [26] and PAW pseudopotentials provided in VASP [27] for all of the calculations. We determined the cutoff energy and k-mesh for the total energy to be converged within 0.5 kJ/unitcell for the both Cu 3 Ti-type and D0 19 -type domains. We calculated the vibrational frequency of the tracer in the jump direction with the harmonic approximation of the potential surface.
The diffusion routes in D0 19 -type (Cu 3 Ti-type) domain is shown in Fig. 2 (Fig. 3). The diffusion barriers for D0 19 -type domain are written in the figure. Those for the Cu 3 Ti-type domain are written in Tab. I with the vibrational frequencies ν in the jump direction and jump rate η at 300 K and 423 K. It can be observed from the table that the route 3 of D0 19 -type domain and the route 6, 6' of Cu 3 Ti-type have the highest barrier energies, which is the ground of the domain division of the diffusion network. D0 19 -type domain has only one type of Cu site and the vacancy formation energy was calculated as 16.60 kJ/mol. Cu 3 Ti-type domain has two types of Cu site denoted by 'top' and 'base' in Fig. 3. The vacancy formation energies are calculated as 25.94 kJ/mol for 'top' and 21.71 kJ/mol for 'base'.

COARSE GRAINING
In the structure of D0 19 -type domain, the hexagonal tubes consisting of Cu sites surrounded in the lines connecting Sn sites. Each of the tubes consists of the laminated Cu ions arranged in the triangle shape located in the different abplanes by 1/2 in fractional coordination (hatched and unhatched sites). Since the barrier energy of route 3 connecting the two tubes is much higher than the others, the ions diffuse in each of the independent tubes. In addition, route 1 has much lower barrier energy than route 2 so the vacancy can move inside of the triangle through route 1 much more easily than moving across the triangle through route 2. Therefore, it  Fig. 3. ∆E i is the energy barrier, ν i is the vibration frequency towards the jump direction, and eta i is the ion jump frequencies at 300 K and 423 K. 103.66 --- The diffusion routes of Cu ion in Cu 3 Ti-type domain. The blue balls represent Cu ions and gray balls do Sn ions. The hatched balls are located on c=1/2 plane and the unhatched ones are on c=0, 1 plane in fractional coordinate. The diffusion routes are shown as the red arrows (1,1,',2) and blue ones (5) and they are named routec and a respectively. The both-side arrow represents that the ion jumps in normal and reverse directions are equivalent. The ion diffusion occurs in the 2-dimensional layers partitioned by the thick black lines, since the energy barriers of the route 6,6' across the lines are much higher than the others. This picture is drawn with VESTA [23].
would be possible to represent the triangle with just a site, and the diffusion network is coarse-grained into the 1-dimensional one long c-axis composed of the representative sites, named rep-site (a schematic picture given in Fig. 4). The correlation factor of 1-dimensional diffusion is known to be zero in the textbook [19] and hence The diffusion network of Cu 3 Ti-type domain can be also coarse-grained with the same analogy yet into the 2dimensional one: The barrier energy of route 6,6f is much higher than the others so the ion diffusion occurs in the layers partitioned by the thick lines written in Fig. 3. In addition, the barrier energies of route 3,3f are much lower than the others, so the three sites connected by route 3,3f can be represented by a rep-site. The vacancy formation rate of the rep-site would reasonably be evaluated as After the coarse-graining with using a rep-site, the diffusion network becomes the 2-dimensional one consisting of two types of diffusion routes a,c shown in Fig. 5. (Noted that route a is parallel to a-axis but routec deviates from c-axis.) Route a consists of route 5 and routec consists of route 1,1',2 so the jump rates through routes a,c can be given as: We define γ (top) , γ (top) as: and substitute γ (top) , γ (base) into eq. (12), (13). Then, we can get the effective jump rates η a , ηc as: a-axis Noted that γ (top) and γ (base) correspond to the probabilities of that the vacancy found in a top (base) site after the vacancy wondering in a rep-site for a sufficiently long time.
Denoting the self-diffusion coefficients in direction a, c as D(a), D(c), the angle averaged self-diffusion coefficient of ε-Cu 3 Ti-type domain can be evaluated as Once the D * a and D * c are obtained from eq. (2), they are projected into D(a) and D(c) based on eq. (1): Here, θ ac (θ cc ) is the angle between the direction a (c)-axis and routec.

FORMULA OF CORRELATION FACTOR
In order to obtain D * a,c from eq. (2), we have to evaluate the correlation factors for both of the route a and routec. We established the equations alternating eq. (8) for multiple types of diffusion routes.
Before explaining about it, we will confirm the notation rules used in this paper. Suppose that the tracer (open cross) has just moved through route l = a orc shown as the arrows towards the tracer in Fig. 6. The six sites surrounding the tracer are indexed with k = 1 − 6 in clockwise order beginning from the site where the tracer was positioned before moving. θ k in eq. (5) is the angle between the directions of the sequential tracer moves. Hence, for example, θ 2 = θ ac for ion move l = a and θ 2 = 2θ cc for ion move l =c. For convenience, the l dependency is noted as θ (l) k and a new function L (l) (k) is also introduced, which returns the type of the diffusion route (a orc) connecting site k and the tracer position. The return values are listed specifically as:  The nth-order average cosign cos θ (n) l=a,c may be built as the below recurrence formula similar to the idea for eq. (8): This equation can be solved analytically and, of course, numerically also. The correlation factors for the ion jumps through route a and routec is given as same as eq. (6) f a,c = 1 + 2 which is same as (6). Here, n max should be large enough for f a,c to be converged.

EVALUATION OF CORRELATION FACTOR
Before solving eq. (19) to obtain the average cosines cos θ (n) l=a,c , P (l) k should be prepared first, which represents the probability that the tracer gets pulled back by the vacancy from site k, after the tracer moved through route l. P (l) k is given as the sum of the realization probabilities for the vacancy tracks pulling back the tracer from site k. This value can be reasonably approximated with considering only the vacancy tracks having highest realization probabilities. This is an example to obtain P (l) k : The tracer just moved to the central site (open cross) through the route a (double shafted arrow) in Fig. 7. Suppose that the vacancy moves through a route a,c with the selection probability p a,c (given later) and P (l) k is given as the product of the ones of the routes included in the vacancy track. Thus, the realization probability becomes lower when the track consists of larger number of the diffusion routes in general. Now, we consider the two vacancy tracks shown as the red and blue arrows in Fig. 7. The former !" !" !" !"  consists of two routesc and the latter does of one route so the realization probabilities of the both tracks is given as p 2 c and p a respectively. Hence, P (a) 2 =p 2 c + p a . p a,c is in proportion to η a,c (T ) and given as p a,c (T ) = η a,c (T ) 2η a (T ) + 4ηc(T ) .
The denominator reflects that a rep-site connects to 2 rep-sites through route a and 4 rep-sites through routec. p a,c depends on temperature due to the Boltzmann factors in η a,c , so the different vacancy tracks should be taken into account for different temperature. We considered the high and low temperature cases separately.
First, for the high temperature case, the relationship, p a ∼ pc ∼ 1/6 (e.g., p a /pc = 1.03 at T =1000 K), is given from ab initio calculations. Hence, p l = 17%, p 2 l = 3.8%, and p 3 l = 0.46% so the vacancy tracks including 3 moves or more would be ignored Then, the vacancy cannot pulls back the tracer from site k = 3 ∼ 5 for the preceding ion move through route l = a and hence P (a) k=3∼5 = 0. The other P (l=a) k are also straightforwardly given as P (a) k=1 = p a , P (a) k=2 = pc · pc, and P (a) k=6 = pc · pc. Applying the same analogy for l =c, P (c) k=3∼5 are zero and the others are P (c) k=1 = pc, P (c) k=2 = p a · pc, and P (c) k=6 = pc · p a . Secondly, for the low temperature case, p a is much greater than pc (e.g., p a /pc = 8.93 at T =300K). Thus, we suppose that, while the number of the vacancy moves through routec is limited in 0 or 1, that through route a is not limited. Then, we have to consider the infinite patterns of the vacancy tracks and it is apparently impossible. We introduced the following three tools to count up the primary vacancy tracks which have the high realization probabilities: • Suppose to consider the vacancy tracks as the vacancy being finally at the starting site after the arbitrary n times of the vacancy moves through route a. The sum of the probabilities of the vacancy tracks is denoted by π even .
• Suppose to consider the vacancy tracks as the vacancy being finally at the left (right) site of the stating site after the arbitrary n times of the vacancy moves. The sum of the probabilities of the vacancy tracks is denoted by π odd .
• Suppose to consider the vacancy tracks as the vacancy being finally at the left (right) site of the stating site after the arbitrary 4n times of vacancy moves, imposing that the vacancy never gets into the sites positioned in the right (left) side from the starting site. The sum of the probabilities of the vacancy tracks is denoted as π oneside .
Here, π even and π oneside must be greater than 1 because the vacancy has to be the starting site for n = 0. It is explained in the appendix how to deduce π even , π odd , and π oneside . P (l) k is deduced here for every (k, l) using the three tools: (1) l = a, k = 1 case; We may consider the vacancy tracks that the vacancy roams in the left side from site k = 1 (including site k = 1) through route a before the pull-back. The probability sum [28] of the vacancy tracks before the pull-back from site k = 1 is exactly π oneside . Hence, P (a) k=1 = π oneside · p a . (2) l = a, k = 2 ∼ 6 case; Two vacancy moves through routẽ c or more are required to pull back the tracer from these site. Hence, P (a) k=2∼6 = 0. (3) l =c, k = 1 case; We may consider the vacancy tracks that the vacancy roams through route a before the pull-back. The probability sum of the vacancy tracks before the pull-back from site k = 1 is exactly π even . Hence, P (c) k=1 = π even p c . (4) l =c, k = 2 case; This case is quite similar to case (3). The difference is just that the vacancy pulls back the tracer from the site k = 2 not k = 1 after roaming through route a. The probability sum just before the pull-back is exactly π odd . Hence, P (c) k=2 = π odd p c . (5) l =c, k = 3 case; We may only consider such vacancy tracks that the vacancy moves to site k = 2 from k = 3. The probability sum before the vacancy positioned at site k = 3 is given as π odd p c just like case (4), and that of the left vacancy tracks is given as π oneside p a just like case (1). P (c) k=2 is eventually given as the product of them: P (c) k=2 = π odd p c · π oneside p a . (6) l =c, k = 3 ∼ 5 case; Any of the vacancy tracks have to include two moves through routẽ c or more. Hence, P (c) k=3∼5 = 0. (7) l =c, k = 6 case; We may consider the two patterns of vacancy tracks before the vacancy positioned at the site k = 6: The vacancy moves to k = 1 from k = 6 or site Y (shown in Fig. 8). The probability sum before reaching site k = 6 is given as (π even + π odd ) pc. and that of the left vacancy tracks is equivalent to the one in case (1). P (c) k=6 is eventually given as the product of them: P (c) k=6 = (π even + π odd ) pc · π oneside p a . We discussed the high temperature and low temperature cases separately above. If considering the sum group of the vacancy tracks taken in the both cases to calculate the probability sums P (l) k , such P (l) k may work well in the entire range of temperature. The vacancy tracks for the high temperature case are actually included in the ones of the low temperature case except for l =ã, k = 2, 6. Eventually, P (l) k is given as: π even · pc (l =c, k = 1) π odd π oneside · p a pc (l =c, k = 3) (π even + π odd ) π oneside · p a pc (l =c, k = 6) 0 (otherwise) (22) and eq. (19) becomes from which the nth-order average cosigns cos θ (n) a,c are given for any order n. Fig. 9 shows the convergence of f a,c for n max . f a,c apparently vibrates for n max because every pull-back denies the preceding tracer move. Further on, we consider only the converged values, which is actually used to solve the model.

RESULTS AND DISCUSSION
If route a andc are equivalent, the coarse grained structure shown in Fig. 5 is regarded as ideal 2-dimensional hexagonal lattice. Its correlation factor is analytically known as 0.56 [19] (shown as the dotted line in Fig. 9), so it is interesting to compare f a,c with the ideal value. First, in the case of low temperature (T =300 K), it is observed that f a is lower and fc is higher than 0.56 respectively. This fact comes from that the vacancy moves through routes l = a much more easily than through route l =c: The vacancy easily goes away from the tracer through route a after the tracer move l =c, so the pullback seldomly happens and fc is becomes higher. On the other hand, for the tracer move l = a, the vacancy moves through route a on the same line for many times, so the pull-back often happen and f a becomes lower. In another way of thinking, the latter situation is similar to the 1-dim diffusion so f a gets much lower. Secondly, in the case of high temperature (T =1000 K), it is expected that both f a and fc converge around 0.56 when temperature increasing, because route a and c are almost equivalent at high temperature and the diffusion network is similar to that on ideal 2-dim. hexagonal lattice.
The converged values are actually higher than 0.56 by ∼0.08 probably because of that the vacancy tracks including 3 routes and more are ignored, which can cause the overesimation of the self-diffusion coefficient. Nevertheless, the bias is quite small even at high temperature and presumably it should get less at low temperature because the number of the route a is not limited.
Once getting the correlation factors f a,c , D Cu 3 Ti can be calculated from eq. (16)∼ (18). Since ε-Cu 3 Sn phase consists of Cu 3 Ti-type and D0 19 -type domains with the ratio of (M −2):2, the self-diffusion coefficient in the phase is given as The temperature dependency of D ε−Cu 3 Sn is shown in Fig. 10 for both M = 4 and M = ∞.  FIG. 9. The convergence of correlation factor f a,c for 300, 500, and 800 K evaluated from eq. (20). The black dotted line represents the correlation factor of 2-dimensional hexagonal lattice and the value is 0.56. [19].  [12], e: Expt. [13], f: Expt. [14], g: Expt. [15], h: Expt. [16], i: Expt. [17], This graph compares our predictions for the diffusion coefficients (blue points) with the ones reported from the experiments (black points) and classical molecular dynamics (red points). The blue open (close) points correpond to the analysis with M = ∞, (4).
MD [18] overestimates our results by a digit. It is possible that the MD underestimated the potential top of diffusion route because of too large time step, since the time step is not shown in the previous paper. [18] However, useally the time step is femto-second order and it is sufficiently small compared with the time scale of the vibration shown in Tab. I. Hence, we tentatively conclude that the overestimation comes from the practical force-field, which is presumely less reliable especialy for the description of the saddle point since three ions or more interact with the tracer.

SUMMARY AND PROSPECTS
We established a new modeling scheme for ion selfdiffusion coefficient by introducing some novel concepts and formalisms. The scheme is applied to obtain the Cu selfidffusion coefficient in ε-Cu 3 Sn phase of Cu-Sn alloy. The highlights are 'domain division' and 'coarse grainining' of the diffusion network based on the barrier energies predicted by the ab initio calculation. In the former concept, the diffusion network is divided into the disjunct domains. Then, the calculation for the original large unitcell is replaced by the ones for the different types of simple domains, which significantly reduces the calculation cost. The diffusion network of our target system, ε-Cu 3 Sn, is divided into the two dypes of domains, which have 1-and 2-dimenstional networks respectively. It is concluded based on the textbook [19] that the former domain never contributes to the self-diffusion. In the latter concept, the diffusion network of each domain is coarse-grained by representing the ion sites with just a site, which are connected by the routes having the lowest barrier energies. This is on ground of that the vacancy moves through such routes much more easily than through the others. Eventually, it is significantly simplified to count up the vacnacy tracks. After introducing the concepts, the model for ε-Cu 3 Sn phase was successfully established and much better preditions by a digit than the previous classcial MD study [18] in the comparison with the experimental values.
If the overesimation of self-diffusion coefficient in the previous classical MD work [18] comes from the usage of the practical potential, ab initio MD (AIMD) can be a strong competitor of our modeling scheme. Nevetherless, our scheme is superior to AIMD at the following points: First, the quantities included in the model can be evaluated from more reliable methods such as diffusion Monte Carlo method. [29] The merit is remarkable especially for such systems as transition metal dioxide [30,31], where the severe accuracy is required to evaluate the electronic correlation effect. Secondly, it is challenging to evaluate the self-diffusion coefficient with AIMD in terms of the conputational cost, because ion jump is a rare event. The recent attempt using the color-diffusion algorithm with non-equiblium MD realizes more efficient simulation by increasing the rate of ion jumps. [32] Nevertheless, 1100 steps were used to obtain the self-diffusion coefficient in body centered molybdenum crystal at 1600 K. To make matters worse, more steps are required at lower temperature because ion jump occurs less frequently. The application of AIMD at low temperature is still challenging. On the other hand, while our scheme requires to construct the model by hand, the self-diffusion coefficient can be calculated with the reasonable cost. The merit becomes remakable especially when performing a lot of similar calculations: For example, surveying how the addicted ions affect the self-diffusion coefficient with modeling the inclusion of the addictinves using virtual crystal approximation. [33].