Article pubs.acs.org/JPCA
Molecular Dynamics Simulation of Self-Diffusion Processes in Titanium in Bulk Material, on Grain Junctions and on Surface Gennady B. Sushko,† Alexey V. Verkhovtsev,†,∥ Alexander V. Yakubovich,† Stefan Schramm,†,‡ and Andrey V. Solov’yov*,§,‡,∥ †
Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt am Main, Ruth-Moufang-Str. 1, 60438 Frankfurt am Main, Germany ‡ Department of Physics, Goethe-Universität Frankfurt am Main, Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany § MBN Research Center, Altenhöferallee 3, 60438 Frankfurt am Main, Germany ABSTRACT: The process of self-diffusion of titanium atoms in a bulk material, on grain junctions and on surface is explored numerically in a broad temperature range by means of classical molecular dynamics simulation. The analysis is carried out for a nanoscale cylindrical sample consisting of three adjacent sectors and various junctions between nanocrystals. The calculated diffusion coefficient varies by several orders of magnitude for different regions of the sample. The calculated values of the bulk diffusion coefficient correspond reasonably well to the experimental data obtained for solid and molten states of titanium. Investigation of diffusion in the nanocrystalline titanium is of a significant importance because of its numerous technological applications. This paper aims to reduce the lack of data on diffusion in titanium and describe the processes occurring in bulk, at different interfaces and on surface of the crystalline titanium.
■
INTRODUCTION Grain boundaries, which are interfaces between differently oriented crystals of the same material, play a crucial role for many properties of materials.1 For instance, it has been discussed that the phenomenon of grain boundary sliding may be considered as one of the dominating deformation mechanisms of nanostructured metallic materials,2 which are currently in the focus of attention in view of their application for medical implants.3,4 The diffusion of atoms in grain boundaries is several orders of magnitude faster than that in a lattice. Thus, grain boundary diffusion can affect many processes in materials, including microstructure development, phase transformations, plastic deformation, fracture, and so on.5 Despite some experimental data, which have been accumulated over past years, many aspects of the grain boundary diffusion remain not well understood, including the atomistic level description of the diffusion mechanisms. We investigate the process of self-diffusion of titanium atoms in a nanostructured crystalline sample by means of full-atom molecular dynamics (MD) simulations. Because of the chosen geometry of the sample, we specify five characteristic regions, namely, the sample’s crystalline interior, double and triple junctions, surface of the sample, and its edge. We evaluate the temperature dependence of the self-diffusion coefficient of titanium atoms with respect to their location in the sample. The analysis of the performed MD simulations demonstrates that the diffusion coefficient can significantly vary from one chosen © XXXX American Chemical Society
region to another. These changes can reach several orders of magnitude at a given temperature. From the experimental perspective, the number of studies devoted to diffusion in titanium is rather limited. During past decades, a few experiments were performed to study diffusion in both low- (hcp α-Ti)6 and high-temperature (bcc β-Ti)7 phases of crystalline titanium as well as in the temperature region above its melting point.8,9 The diffusion coefficient on a clean surface of titanium has not yet been measured experimentally due to its chemical activity at the processing temperatures and related challenges in maintaining cleanness of the surface.10 Up to now, the atomistic-level studies of diffusion in nanostructured materials have received much less attention as compared with widely developed continuum models.11 Because the amount of available data on diffusion along the grain boundaries and triple junctions in titanium is also very limited, the diffusion process in nanocrystalline materials deserves further study. In this work, we carry out such an analysis and describe the self-diffusion processes occurring in the bulk, at different interfaces and on a surface of titanium nanocrystals. Our studies demonstrate and quantify the significant increase in Special Issue: Franco Gianturco Festschrift Received: April 17, 2014 Revised: June 23, 2014
A
dx.doi.org/10.1021/jp503777q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
such a potential can be applied not only for monatomic systems but also for bimetallic alloys.20,21 In the latter case, the aforementioned parameters, d ≡ dαβ, and so on depend on the type of an atom, α or β, chosen within the summation. To describe the interaction between titanium atoms, we utilize in this study the parametrization proposed by Lai and Liu,22 which reproduces main mechanical and geometrical properties (such as the cohesive energy, lattice parameters, and elastic constants) of an ideal titanium crystal and a nickel− titanium alloy derived from zero-temperature ab initio calculations. This force field has also been applied recently to study surface diffusion at the interface of nickel and titanium crystals23 as well as to evaluate mechanical properties of ideal crystalline and nanostructured titanium and nickel−titanium samples by means of MD simulation of nanoindentation.20,21,24 The performed theoretical analysis and comparison with experimental data have demonstrated that this force field is capable to reasonably describe mechanical properties of these materials. The utilized values of the parameters are as follows: d = 2.95 Å, A = 0.153 eV, p = 9.253, ξ = 1.879 eV, and q = 2.513.22 Most of the many-body potentials approach zero at large distances, and thus a cutoff radius rcut is frequently introduced to reduce the computation time. In this case, the interatomic potentials and, subsequently, the forces are neglected for atoms positioned at distances larger than rcut from each other. The utilized parameter set22 was constructed with a fixed cutoff radius of 4.2 Å that slightly exceeds the second-neighbor distance in the hcp-Ti crystal, d2 = 4.17 Å. Thus, interaction of a given atom with its second neighbors is not truncated in the case when the simulation is carried out at very low temperature. Because of thermal expansion of the crystalline structure at elevated temperature, the second-neighboring atoms can leave the interaction sphere that may cause discontinuity of the potential. To avoid this effect, in this study, we have chosen a significantly larger cutoff radius, namely, 7 Å. Computational Approach. Numerical simulations of complex atomic and molecular systems require substantial computational time. To reduce it, a GPU version of the MBN Explorer25 software was utilized in this study. MBN Explorer is a software package for the advanced multiscale simulations of complex molecular structure and dynamics.25 It is suitable for structure optimization,26−28 simulation of various dynamical processes,23,29,30 and self-organization31−33 in various molecular and bio/nanosystems. Both CPU and GPU versions of MBN Explorer utilize the so-called linked cell algorithm25 to reduce the computational time. This algorithm splits the simulation box into cells of a fixed size and distributes all particles in the system among these cells. In the CPU implementation, the calculation of forces acting on each particle and the particle’s energy can be performed by summing up the interactions with particles in the same cell as well as in all neighboring cells. In a slightly different manner, the multithreaded GPU version of the code defines the size of the box cells to optimize thread block sizes, as defined by the GPU architecture. Each particle sums up its forces and energy in a separate thread and all particles in a single cell define a thread block. Because the box size might not be congruent with the range of the interaction, this can lead to a neighborhood of cells that might span a larger range than just near-neighborhood cells. In the CPU version of MBN Explorer, the list of cells is distributed between all available cores of CPU to speed up the
the diffusion coefficient at the grain boundaries and at the triple junctions, which should be the very useful input information for the construction of the continuous models of the diffusiondriven processes in materials.
■
THEORETICAL AND COMPUTATIONAL METHODS Interaction Potential for Titanium Atoms. The diffusion process is studied in this work by means of classical MD simulations. Utilizing this approach, one can describe the time evolution of a system of interacting particles by numerically solving classical (Newton or Langevin) equations of motion. This technique allows one to ultimately describe the dynamics of systems consisting of up to several million atoms at substantially large (up to hundreds of nanoseconds) time scales. Within this framework, interactions between constituent atoms are parametrized by various empirical potentials of the force fields. The interaction between titanium atoms is described in this study using the many-body Finnis−Sinclair potential.12 As most of interatomic many-body potentials,13−17 it is constructed as a sum of an attractive long-range density-dependent term and a repulsive short-range term, which results from the repulsion between core electrons of neighboring atoms. In the Finnis−Sinclair representation, the total energy of an N-atom system is postulated in the following form U=
1 2
N
∑ ∑ V (rij) − c ∑ i=1 j≠i
ρi (1)
i
where ρi =
∑ ϕ(rij) (2)
j≠i
The first term on the right-hand side of eq 1 is the sum of pair interactions, while the second term represents the many-body interactions responsible for a significant part of bonding in metals. This term represents the embedding energy of an atom i in the host electron density ρi induced at site i by all other atoms.5,12 The function ρi is given in the Finnis−Sinclair representation by a sum of additional pairwise functions; see eq 2. The form of the embedding function is postulated in this approach to be attractive and proportional to the square root of ρi,5 with c being a positive constant. The square root form is chosen to mimic the result of tight-binding theory, where ϕ(r) is interpreted as a sum of squares of overlap integrals.12 In this approach, the energy of the d electron band in metals is proportional to the square root of the second moment of the density of states.17−19 Similar to the original second-moment approximation of the tight-binding (TB-SMA) scheme,17−19 the functions V(rij) and ϕ(rij) are introduced in exponential forms. Thus, the total potential energy of a system of N titanium atoms, located at positions ri, in the Finnis−Sinclair representation reads as N
U=
N
∑[∑ i=1
N
A e−p(rij / d − 1) −
∑
j=1
j=1
(i ≠ j)
(i ≠ j)
ξ 2e−2q(rij / d − 1) ] (3)
where rij is the distance between atoms i and j, and the rest are free parameters of the potential: d is the first-neighbor distance, parameters ξ and q are related to the hopping integral, and p is related to the compressibility of the bulk metal.17 Note that B
dx.doi.org/10.1021/jp503777q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Exclusion of the Pure Oscillatory Motion. To exclude oscillations with larger periods, one needs to exclude the contribution from particles, whose position in the crystal does not change. We have analyzed the trajectory of each particle to determine the maximal displacement of a particle and a characteristic distribution volume of its positions, as illustrated in Figure 1. This value was then compared with the volume of a
calculation. In the GPU version, two levels of parallelism are used: cells are distributed between cores and then particles in each cell are distributed between SIMD threads. For both versions, this approach leads to a linear growth of the computational time with respect to the number of particles in the system with increasing the system size. For the system used in the present calculations (CPU: AMD Opteron 6172, 2.1 GHz; GPU: Radeon HD 5870), a speed-up factor of 22 for GPU execution was achieved as compared with a single CPU-core. Although the acceleration effect of GPUs is the most pronounced in systems with long-range (e.g., coulomb) interactions, there is a still a sizable gain in computational speed in favor of GPUs as compared with CPU-based approaches. Calculation of the Diffusion Coefficient. The diffusion coefficient of a particle is defined as follows34 D=
⟨Δr 2⟩ 2zΔt
(4)
Figure 1. 2D projection of simulated trajectories of titanium atoms (red); ten uppermost atomic layers of a crystal were considered. The dashed line represents an edge of the constructed cylindrical sample. (See the description below.) Particle 1 (blue area) oscillates around its equilibrium position and does not leave the characteristic distribution volume during the simulation. Particle 2 (green area) oscillates and also diffuses to a neighboring position. The inset illustrates the trajectory of this particle during the simulation.
where ⟨Δr2⟩ is a mean-square displacement of a particle per time Δt averaged along its trajectory ⟨Δr 2⟩ =
1 nt
nt
∑ [R(t0j + Δt ) − R(t0j)]2 j=1
(5)
where R(t) is a position vector of a particle at time t, and nt is the number of time slices, t0j, considered along the trajectory. The parameter z is defined by the dimensionality of space.34,35 In the case of particle diffusion along or over a surface z = 2, while it equals 3 in the case of bulk diffusion. This definition of the diffusion coefficient includes displacement of atoms caused by different types of motion, namely, by diffusion itself as well as thermal vibrations and collective motion of the sample. To study the diffusion process directly, all other types of motion should be excluded from consideration. In this study, exclusion of collective motion is done by subtraction of motion of the center of mass of the system. Exclusion of the oscillatory motion requires additional efforts and is done as follows. Because of a finite simulation time, the value of the calculated diffusion coefficient has a lower limit, which is defined as D=
⟨r02⟩ 2zτ
unit cell. If the particle stayed in the unit cell, its motion was considered to be an oscillatory one (see Figure 1), and this particle was excluded from the diffusion coefficient calculation.
■
NUMERICAL SIMULATIONS To study the diffusion of titanium atoms at various interfaces, the following sample was constructed. As illustrated in Figure 2, the sample has a cylindrical form and is split into three equal sectors. The cylinder has a radius of 7.7 nm and a height of 10 nm, consisting of about 105 000 atoms. Each sector of the cylinder is constructed of titanium atoms packed in the hcp crystalline lattice. The lattices in each sector are rotated with respect to each other by 120°. High symmetry
(6)
where r0 is a characteristic size of the unit cell, z is the dimensionality of space, and τ is the total simulation time. This estimation means that oscillations with the period larger than the simulation time cannot be separated from the diffusive motion. Diffusion Coefficient Extrapolation. To distinguish the diffusion process from thermal vibrations of smaller periods, the diffusion coefficient calculated in the course of simulation can be fitted with the following function a D(Δt ) = + D0 (7) Δt where the first term on the right-hand side corresponds to the oscillatory motion of atoms and D0 is the desired diffusion coefficient. Calculation of the diffusion coefficient in the course of simulation results in a series of D(Δt) values. This series can be fitted by the aforementioned dependence and extrapolated to larger values of t to obtain the desired value of the diffusion coefficient D0.
Figure 2. Structure of the investigated sample. The sample consists of three equal sectors of crystalline hcp titanium rotated with respect to each other by 120°. C
dx.doi.org/10.1021/jp503777q | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Figure 3. Five regions of the titanium sample where different types of diffusion were studied. We considered diffusion in bulk material (blue area/1), on double and triple junctions (red/2 and dark-yellow/3 regions, respectively) as well as diffusion on the surface (dark-green area/4) and border (light-green area/5) regions of the sample.
of the sample leads to stability of the junctions in the course of simulations. Prior to the simulation of diffusion, the initial structure of the sample was optimized. In the course of the optimization procedure, the displacement of atoms was
