Search

Collections

Journals

About

Contact us

My IOPscience

The crystallization process of liquid vanadium studied by ab initio molecular dynamics

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 J. Phys.: Condens. Matter 26 155101 (http://iopscience.iop.org/0953-8984/26/15/155101) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 129.10.52.112 This content was downloaded on 07/04/2014 at 08:19

Please note that terms and conditions apply.

Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 26 (2014) 155101 (7pp)

doi:10.1088/0953-8984/26/15/155101

The crystallization process of liquid vanadium studied by ab initio molecular dynamics T T Debela1, X D Wang1, Q P Cao1, D X Zhang2 and J Z Jiang1 International Center for New-Structured Materials (ICNSM), Laboratory of New-Structured Materials, State Key Laboratory of Silicon Materials, and Department of Materials Science and Engineering, Zhejiang University, Hangzhou, 310027, People’s Republic of China 2 State Key Laboratory of Modern Optical Instrumentation, Zhejiang University, Hangzhou, 310027, People’s Republic of China 1

E-mail: [email protected] Received 21 January 2014, revised 19 February 2014 Accepted for publication 24 February 2014 Published 27 March 2014 Abstract

We present a study of the crystallization process in liquid vanadium over a temperature range from 3000 K down to 1500 K using ab initio molecular dynamics simulations. Short-range order evolution during solidification is studied using various structural analysis methods. We show that the icosahedral-like short-range order is detected in the stable liquid phase and grows upon supercooling. The system undergoes a first-order phase transition (from a liquid to a solid state) at a temperature of about 1600 K. The crystal nucleation process is further studied using the time–temperature transformation mechanism by annealing the system at 1650 K. The nucleation is examined using bond-orientational order and density fluctuation analysis. Our finding is that various precursors appear in the region of high bond-orientational order with the majority having body-centered cubic (bcc)-like symmetry. This bcc-like region grows on annealing via thermal fluctuations. Our results reveal that the bond-orientational order precedes the density fluctuation, and is the main driving factor for nucleation. Keywords: liquid vanadium, atomic structure, ab initio molecular dynamics, nucleation, bond-orientational order (Some figures may appear in colour only in the online journal)

1. Introduction

order (ISRO). Over the last few years, considerable progress has been made, both experimentally [4–13] and theoretically [14–21], in understanding the atomic structure and crystallization of metallic melts, and the interface between the liquid and solid state. However, our understanding of atomic-scale crystal nucleation and growth is very limited. It is difficult to study this process experimentally, owing to the short lengths and time scales associated with nucleation events. To overcome this difficulty colloidal systems (in which particles are large and slow enough to access a nucleation process with particle-level resolution) have been used to investigate the elementary process of crystallization [22–25]. In addition, ab initio molecular dynamics (AIMD) simulations are used to shed light on the crystallization process of realistic systems with reasonable accuracy [16, 26].

An adequate understanding of microstructure evolution during the solidification of metallic melts is a key topic in both condensed matter and materials science [1]. To fully understand this phenomenon, one must start with an accurate description of the liquid structure and track the process of crystallization as a function of temperature. However, our understanding of this mechanism is far from complete owing to the lack of longrange order in liquid and supercooled liquid states, resulting in a complex atomic structure. Following the work of Turnbull in 1950 [2], it has been widely accepted that metallic melts can be cooled to well below their melting temperatures. Frank [3] first proposed that the local short-range order (SRO) in most supercooled metallic melts is icosahedral-like short-range 0953-8984/14/155101+7$33.00

1

© 2014 IOP Publishing Ltd Printed in the UK

T T Debela et al

J. Phys.: Condens. Matter 26 (2014) 155101

(VASP) [27, 28]. The electronic exchange and correlation are described in the generalized gradient approximation (GGA), using the PW91 functional due to Wang and Perdew [29]. The ultrasoft pseudopotentials (USPP) of Vanderbilt type were used to describe the electron–ion interaction [30, 31]. The system consists of 128 V (vanadium) atoms in a cubic supercell with standard periodic boundary conditions. In the process of sample preparation, it was first melted and equilibrated at 3000 K, which is well above the experimental melting point of V (i.e. Tm = 2163 K) in order to avoid the memory effect from the initial random configuration. Then, the system was cooled down to 2700, 2500, 2400, 2300, 2200, 2000, 1800, 1700, 1600 and 1500 K respectively with a cooling rate of 3.33 × 1013 K s−1, in an NVT ensemble by means of a Nosé thermostat to control the temperature [32], during which the cell size was adjusted to give nearly zero pressure at each temperature. The equation of motion is solved via the velocity Verlet algorithm with a time step of 3 fs. Only the Γ point is used to sample the Brillouin zone of the supercell. After equilibrating for 4000 MD steps (12 ps) at each temperature, the final 2000 MD steps (6 ps) were used to collect the trajectories for statistical analyses. The TTT calculation is performed at the temperature between the supercooled and stable liquid state. By first checking the reliability and suitability from a computational time (CPU time) perspective, we chose the temperature 1650 K for our isothermal annealing. After equilibrating for 6 ps (in order to adjust the lattice parameters to make pressure close to zero), we continued for another 18 ps for the analysis of crystal nucleation process. 3. Results and discussion Figure 1. (a) Pair-correlation function g(r) as a function of temperature during solidification of liquid vanadium. (b) The evolution of the calculated coordination number CN with temperature.

3.1 General overview of solidification

Structural evolution upon cooling from the melt can be observed through temperature dependence of the pair-correlation function g(r) analysis. Figure 1(a) shows the calculated g(r) at various temperatures. The position of the first peak is at about 2.60 Å and moves towards the large r value by about 0.015 Å as temperature drops. This value is in good agreement with the experimental interatomic distance of vanadium at ambient condition. During solidification, it is observed that the first peak gets narrower and its height increases indicating the development of the SRO in the first coordination shell of the supercooled or liquid state. The radius of the first coordination shell is defined by a cutoff distance rmin, which is taken to be the first minimum of g(r), is 3.60 Å. This value will be kept throughout the structural analysis. A rapid change occurs at the temperature of 1600 K, where the second and third peaks split into two sub-peaks, and the fourth and fifth peaks become clear. The first peak of g(r) also shows a possibility to decompose into two sub-peaks. As the temperature further drops, these features become amplified. These results indicate that crystallization occurs at 1600 K in supercooled liquid vanadium and the crystalline structure becomes more dominant at lower temperatures. The nearest neighbor coordination number (CN), by integrating the radial distribution functions 4πnr 2g ( r ) (where n is the atomic number density) up to rmin are plotted in figure 1(b). The CN linearly increases from

In this work, we present a detailed crystallization process in metallic melts, taking vanadium as an example. Using AIMD simulations, we show the structural evolution of stable liquid and supercooled liquid states during solidification. The liquid– solid phase transition is a first-order phase transition in which the supercooled liquid state suddenly transforms to the crystalline state. However, the kinetics that create the final stable crystalline phase have so far remained unclear. One possible way of obtaining a full picture of nucleation is by performing the so-called time–temperature transformation (TTT) of the state between the supercooled liquid and the final stable crystalline phase. We studied in detail the atomic evolution during the nucleation process at 1650 K as a function of time, in addition to general solidification at a constant cooling rate. In section 2, we give details of our model and simulation methodology. Section 3 contains the results and discussion of general solidification, crystal nucleation and growth processes. A brief summary and conclusions are given in section 4. 2. Computational details The AIMD simulations were performed based on density functional theory (DFT) using Vienna Ab-initio Simulation Package 2

T T Debela et al

J. Phys.: Condens. Matter 26 (2014) 155101

Figure 3. Fractions of HA indices in liquid V at different temperatures.

Figure 2. Structure factors, S(q), of liquid V at different temperatures.

is for the number of bonds between these common neighbors. HA analysis is able to distinguish between various local structures such as bcc, close-packed structures (fcc and hcp) and icosahedral environments (ICOS). For instance, ICOS is characterized by 155 pairs, while 154 and 143 pairs represent distorted ICOS. The fcc and hcp crystalline structures are composed of the basic components of 142-type pairs, whereas 166 and 144 are typical pairs for the bcc crystal. In figure 3, we report fractions of the most abundant bond pairs as a function of temperature during cooling. A fraction of 155 pairs already exists in the liquid state and grows upon supercooling, while that of distorted ICOS shows a slight decrease. This result indicates that as temperature decreases, 2 sin( qr ) 4πr S(q) = 1 + n [ g ( r ) − 1 ] dr , (1) the distorted ICOS presumably become more perfect ICOS. qr We also observed a small fraction of 132 pairs, indicating that where n is the averaged number density, r is the radial the local structure in liquid and supercooled liquid V is based distance and q is the wave number. As the temperature on 13-atom icosahedra. Conversely, at 1600 K, the population decreases, the main peaks become narrower and their ampli- of both ICOS and distorted ICOS indices drops suddenly, tude increases, indicating an increase in short-range order- and tends to be zero with a further decrease of temperature. ing with supercooling. A typical mark of the ISRO is a The fraction of 142 pairs, characteristic of fcc and hcp local shoulder on the high q side of the second peak [33, 34]. It is atomic structures, slightly decreases from 4.6% at 3000 K, observed that the shoulder of the second peak (see the inset of 1.6% at 1700 K, to almost zero below 1700 K. The fractions figure 2) becomes more pronounced for S(q) calculated in the of 144+166 bond pairs slightly increase from 16.6% to 22% temperature range of 1700–3000 K. This phenomenon even with decreasing temperature from 3000 to 1700 K. Suddenly, becomes clearer below the Tm. Taking the position of the their fraction dramatically increases to about 97% at the temmain peak in S(q) as q1 and that of the second peak and the perature of 1600 K. Consequently, one can conclude that the shoulder as q2 and qs , respectively, for the liquid V in super- system has crystallized into a bcc-like structure. cooled state at 2000 K, 1800 K and 1700 K we have found Further investigation of the evolution of atomic local packthat the peaks are related by q2 /q1~ 1.745, 1.727, 1.718 and ing was performed by using Voronoi tessellation analysis qs /q1~ 1.968, 1.963, 1.976, respectively, which is in excellent [36, 37], which is a more three-dimensional (3D) approach than the HA perspective as it provides a more complete agreement with ideal ISRO [35]. To get more detailed information on the local order in geometrical construction of a central atom to its neighthe system upon cooling from the melt, we analyzed the boring atoms. In this approach, the perpendicular bisectemperature dependence of various bond pairs obtained tors between the central atom and neighboring atoms will by Honeycutt–Andersen (HA) analysis, which is shown in form a polyhedron about the central atom; this is termed a figure 3. In the simplified HA analysis [26], the structure Voronoi cell (VC). This cell is pretty similar to the famous is analyzed by pairs of atoms on which three indices are Wigner–Seitz cell in crystallography. The distribution of assigned: the first digit is 1 if the pair of atom is within the the Voronoi polyhedron has been expressed using Schläfli’s cutoff distance rmin; the middle digit stands for the number of notation whereby a polyhedron is denoted by a set of indices common neighbors owned by the pair of atoms; the last digit , where ni is the number of faces with i vertices.

13.27 ± 0.03 to 13.64 ± 0.05 from a stable liquid at 3000 K to the supercooled state at 1700 K. Then, the coordination number rapidly jumps to 13.95 ± 0.02 at the temperature 1600 K and maintains almost a constant value of about 13.96 ± 0.02 at the temperature 1500 K which reveals a sudden structure change at 1600 K from liquid to a body-centered cubic (bcc)like crystalline state. Figure 2 illustrates the calculated structure factor as a function of temperature. The structure factor, S(q), can be calculated from the Fourier transformation of g(r) by

∫

3

T T Debela et al

J. Phys.: Condens. Matter 26 (2014) 155101

and translation of the clusters with respect to each other. Major SROs in the system at various temperatures can be identified. The detailed procedure can be found in [38]. The collective alignment results for vanadium at five representative temperatures are plotted in figure 5. In the liquid and supercooled state, it demonstrates icosahedrallike five-fold symmetry. This is strong evidence that ISRO dominates in the liquid state and develops on supercooling. However, at 1600 K, bcc-like order is clearly detected. This result further confirms the agreement with different structural analyses given above. 3.2 Crystal nucleation

In the previous section, we have shown that liquid vanadium can be undercooled from 3000 to 1700 K, where the ISRO is detected in the liquid state and develops upon supercooling, before it suddenly solidifies at 1600 K. The interesting question is: what has happened in the range of rapid jump (between 1700 and 1600 K)? On the other hand, the possible pathways through which a supercooled liquid state of metallic melts can evolve from the liquid to the solid state are still a mystery. Here, we focus on the crystal nucleation process by isothermally annealing the system at the temperature of 1650 K to uncover the possible nucleation pathways. Figure 6 shows the evolution of the energy as a function of annealing time. It can be clearly seen that there is a significant jump at about t = 8 ps. The inset snapshots give the clear picture of how the disordered atomic structure becomes more ordered upon annealing, consistent with a general picture of crystallization. So far the crystal nucleation process has been explained by using density as a key order parameter. For instance, the Alexander–McTague theory [39] is based on the expansion of the free energy in terms of density fluctuations. Contrary to this densityparameter scenario, the recently proposed order parameter, a bond-orientational parameter, plays a critical role in the crystallization process [40–43]. To identify the crystal structures, we use the coarse-grained local bond orientation order method, which is the most accurate way of identifying crystal structures [44]. In this method, first, a (2l + 1) dimensional complex vector qlm of atom i is defined as

Figure 4. Evolution of the distribution of major VCs for representative temperatures.

The total number of faces of VC is equivalent to the coordination number (∑ n i = CN) for a chosen central atom. In figure 4, the most abundant VC types in the supercooled liquid state at six representative temperatures are presented. In order to improve the statistics, the average over 2000 atomic configurations was taken. In the high temperature range, more VC types are observed. However, we only present the most abundant ones. The Z13 Frank–Kasper polyhedral (FK) with Voronoi indices , and the Z12 of , are the most abundant. The fraction of VCs with a high coordination number increases with decreasing temperature, while those with a low coordination number (, , and ) decrease upon solidification. It is also interesting to note that the fraction of perfect icosahedral clusters increases with decreasing temperature, although it is not the most abundant. The Z14 FK (, and ) which are related to the 144 and 166 bond pairs in the HA analysis also increase with decreasing temperature, indicating agreement between the two methods. At the temperature 1600 K, the fraction of Voronoi index clusters (not shown here), representing bcc clusters, dramatically increases and reaches 94.6% or so, while other VC types vanish, indicating that the supercooled liquid V has crystallized into a bcc-like crystalline structure. The recently developed cluster alignment method [38] is also used to further explore the SRO of liquid V during solidification. In this method, the local environment of each atom is characterized in the sample by an atomic cluster which includes the atom itself and its nearest neighbors. All the atoms within the first shell, which is the first minimum in the pair-correlation function, are included. We then randomly select one cluster from each 2000 configurations (i.e. 2000 clusters from each temperature). Each cluster is made up of 16 atoms. These clusters are overlapped and aligned by the rigid rotation

1 qlm(i ) = Nb(i )

Nb(i )

∑Y

lm

(riĵ ),

(2)

j=1 where Nb ( i ) is the number of the nearest neighbors of atom i, rij is the vector joining two neighboring atoms (two atoms are considered as neighbors if they are separated by a distance smaller than a cutoff distance, analogous to the first minimum of g(r), i.e. 3.60 Å), and Ylm is a spherical harmonic. The local order parameter which is rotationally invariant can be defined as

4

⎡ ⎢ 4π ql (i ) = ⎢ 2l + 1 ⎣

l

∑

m =−l

⎤2 2⎥ qlm(i ) ⎥ . ⎦ 1

(3)

T T Debela et al

J. Phys.: Condens. Matter 26 (2014) 155101

Figure 5. The atomic density contour plot of the final configuration of collective alignment of V at different temperatures. Five-fold icosahedral-like symmetry dominates the local atomic structures of stable and supercooled liquid V above 1600 K while the bcc-like symmetry is observed at 1600 K.

⎛ N (i ) ⎞ 1 ⎜ ⎟ ( ) Qlm(i ) = q i lm ⎟. N (i ) ⎜ ⎝ j=0 ⎠ Note that N(i) is the number of all the neighbors of atom i (including the ith atom itself). Figure 7 shows the coarse-grained bond-orientational order (Q6) and the local density as a function of time. Q6 provides a measure of the amount of crystalline order in the system, although other parameters were also considered. The value of Q6 is close to zero for a liquid state, and gives a nonzero value for the common crystalline structures. Hence, it is often used as an order parameter to detect any kind of crystal structures. The local atomic density of each atom can be calculated by constructing the Voronoi diagram and measuring the volume linked with each atom. Since our calculation is performed in NVT ensemble, we recalculate the local density using the calculated pressure by assuming Boyle’s law. As is clearly shown in figure 7(a), the average bond-orientational order remains almost constant before 8 ps in the range 0.15–0.2, increases clearly from about 0.2 to 0.48 in a time period of 8–8.5 ps, and then remains almost constant after 8.5 ps. These results strongly correlated with the crystallization process in figure 6, indicating the nucleation is mainly driven by bond-orientational ordering of atoms. Although the average density shows slightly increasing behavior up on crystallization, the temporal change in density fluctuations is not systematic, as shown in figure 7(a). On the other hand, the spatial distributions of local atomic density and the structural change upon crystallization do not seem to correlate well. Further insight can be gained by looking at the spatial distribution of atoms having high Q6. One can define the threshold Q6 (≥0.2) since the critical nuclei are formed from atoms with high bond-orientational order. Figure 8(a) shows the spatial distribution of atoms with high Q6 at representative times taken before and after crystallization. Clear development in Q6 is revealed. The spatial local density distribution of high ordered atoms is displayed in figure 8(b). It can be clearly observed that before crystallization takes place, atoms with high Q6 (≥0.2) still have low densities. This shows that the spatial distribution of high bond-orientational ordered atoms and the local density is not symmetric. The local density order parameter lags the bond-orientational order, although the crystal density is greater than the supercooled liquid density shown in figure 8.

∑

Figure 6. Temporal evolution of internal energy of supercooled liquid V annealed at 1650 K. The insets are snapshots of atom configurations at selected annealing times.

The evolution of the local structure during crystal nucleation can be effectively identified by combining this parameter with other types of rotationally invariant parameter, defined as wl =

⎡l ⎢

∑

l l ⎤ m2 m3 ⎥⎦qlm1 (i ) qlm2 (i ) qlm3(i )

m1+ m2 + m3= 0 ⎣ m1

⎛ ⎜ ⎝

⎞ 2 qlm(i ) ⎟ ⎠

3

∑

l m =−l

2

.

(4)

The parameters in the square bracket are the Wingers 3-j symbols. The sum runs from – l to l with the combination of m1 + m2 + m3 = 0. We present our results by using the coarse-grained form of local bond-orientational order as described in [44] and written as ⎡ ⎢ 4π Ql (i ) = ⎢ 2l + 1 ⎣

Wl =

∑

where

⎡l ⎢

∑

m =−l

⎤2 ⎥ Qlm(i ) ⎥ , ⎦ 1

l l ⎤ m2 m3 ⎥⎦Qlm1(i )Qlm2(i )Qlm3(i )

m1+ m2 + m3= 0 ⎣ m1

⎛ ⎜ ⎝

l

∑

l m =−l

⎞ Qlm(i ) 2 ⎟ ⎠

3

2

(5)

,

(6)

5

T T Debela et al

J. Phys.: Condens. Matter 26 (2014) 155101

Figure 7. (a) Temporal change in bond-orientational order parameter, Q6, of each atom (all atoms are plotted). The average value is plotted in a red line. (b) Temporal change in local atomic density of each atom (all atoms are plotted). The average value at each time is plotted in red line.

Figure 9. Temporal evolution of number of bcc- and hcp-like atoms with Q6 > 0.2, defined by W6.

hcp-like (W4>0) from fcc-like (W4 0.2) at representative times before and after crystallization. (b) The special distribution of local atomic density of atoms plotted in (a).

We further extend our investigation to the mechanism of the polymorph selection during the crystallization process. The order parameters W4, W6 and Q4 together with Q6 were used. Atoms with high bond-orientational order Q6 (≥0.2) were selected, since it is believed that polymorphs will be nucleated from them. W6 is a good order parameter to distinguish bcc-like (W6>0) from closed packed (hcp/fcc) (W6