- Email: [email protected]
Chain structure of liquid Se at high temperature and pressure investigated by ab initio molecular dynamics simulations
Journal of Non-Crystalline Solids 358 (2012) 873–879
Contents lists available at SciVerse ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/ locate/ jnoncrysol
Chain structure of liquid Se at high temperature and pressure investigated by ab initio molecular dynamics simulations Y.B. Wang a,⁎, W.S. Dong b, G. Zhao c, J.X. Ding d, S.H. Li a, Y.J. Ge a a
Institute of Intelligent Machines, Chinese Academy of Science, Post Office 1130, Hefei 230031, China School for the Gifted Young, University of Science and Technology of China, Hefei 230026, China Department of Physics, Ludong University, Hongqi Road, No. 186, Yantai 264205, China d New Star Research Institute of Applied Technology in Hefei City, Hefei 230031, China b c
a r t i c l e
i n f o
Article history: Received 25 July 2011 Received in revised form 30 December 2011 Available online 24 January 2012 Keywords: Liquid selenium; Molecular dynamics simulation; Chain structure
a b s t r a c t Using ab initio molecular dynamics simulations, the local atomic structural, dynamic and electronic properties of liquid selenium were studied under different temperatures and pressures. Compared with experimental data, the calculated structure factors and viscosities are acceptable on the whole. Our results indicate that the chain structure of crystalline selenium still exists in liquid state even at high temperature and pressure. The fraction of twofold-coordinated atoms decreases obviously under high pressure while it remains nearly invariable with the increase of temperature. The Peierls-type distorted structures in trigonal Se still reserve in liquid state even under high temperature and pressure. The calculated DOS displays an obvious dip at EF, and the dip becomes shallower with rising temperature. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The properties of element selenium have been extensively investigated for the past several years mainly due to the following two reasons [1–6]. Firstly, crystalline selenium is one of the elements with the most allotropes next to sulfur and phosphorus, and many allotropes of selenium are of great fascination from the viewpoints of the structure and properties [7]. There are five crystalline configurations at ambient conditions: t-Se, α-, β- and γ-monoclinic Se, and rhombohedral Se, all of which are twofold coordinated that bonded into infinite chains, only with a varying nearest-neighbor-distance from 2.32 to 2.37 Å [8–10]. Due to the valence-electronic configuration s2p4, the structure of liquid Se in the thermodynamic equilibrium state at room temperature and atmospheric pressure is the semiconductor trigonal phase, in which the helical chains of atoms run along the axes of the hexagonal cell [7,11]. This structure may be interpreted by the Peierls-type distortion of a simple cubic lattice: two s electrons are at deep level and do not participate in the covalent bonding; the four p electrons occupy two-thirds of the p-band, so that the simple cubic lattice is unstable against the Peierls-type distortion and the six neatest-neighbor bonds are divided into two short intrachain bonds and four long interchain bonds [12,13]. Secondly, selenium is a good candidate to form glass and amorphous structure with an accessible glass transition temperature (Tg ~ 310 K). As shown in his work previously, Phillips has concluded that covalent materials with an average coordination number N ≤ 2.4
⁎ Corresponding author. Tel.: + 86 1 551 5591195. E-mail address: [email protected] (Y.B. Wang). 0022-3093/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2012.01.001
are potential glass formers [14]. Meanwhile, amorphous selenium exhibits wide applications in photocopying areas in virtue of its photoconducting and semiconducting properties [15,16]. The explicit elaboration of liquid structure is also of great application importance in the control the glass-forming process of Se. Upon heating the infinite chain structure with twofold coordination is retained in the liquid state near the melting point (494 K), which is consistent with experimental-observed high viscosity and low diffusion in liquid Se [17]. However, the magnitude of viscosity decreases more than one order while the diffusion coefficient shows a marked increase as temperature is increased to above 1000 K, suggesting that the chain structure may be strongly destroyed [18]. This point of view has also been advanced by several computer simulations [19–20]. Moreover, it was reported that the chain structure was disrupted in liquid Se under high pressure and the average chain length was reduced sharply to about only ten atoms near the critical point (TC = 1860 K, PC = 380 bar) [21]. These temperature/pressure induced structural variations were also witnessed by a semiconductor–metal (SC-M) transition in the measurements of conductivity [22]. Although a large number of experiments and simulations have been performed on liquid Se, the related structural-change mechanisms are still not fully understood. In present work, we have performed AIMD simulations based on density functional theory (DFT) and the pseudopotential method to study the atomic structural and electronic properties of liquid Se under different temperatures and pressures, in order to illustrate the deeply-inside structural-change behaviors. The paper is organized as follows: in Section 2, we describe the methods of present simulations; the results of our simulations and the corresponding discussions
874
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
are presented in Sections 3 and 4 respectively; and finally a short summary is given in Section 5. 2. Computational methods Our simulations were performed with the framework of the density functional theory (DFT) [23], with the generalized gradient approximation (GGA) formulated by Perdew and Wang (PW91) to the exchange-correlation energy [24], and the valence electron-ion interaction was modeled by the ultrasoft pseudopotential (USPP) of the Vanderbilt type [25,26], as implemented in the Vienna Ab initio Simulation Package (VASP) [27,28]. The 90-atom system was put in a simple cubic box with periodic boundary conditions. Only the Γ point was used to sample the Brillouin zone of supercell. The electronic wave functions were expanded in the plane-wave basis set, with an energy cutoff of 160 eV. Our canonical ensemble simulations were performed at 570, 870 and 1370 K, with a Nosé thermostat to control temperature [29]. The densities as measured in Refs. [30] and [31] were adopted in our simulations. The Verlet algorithm was used to integrate Newton's equations of motion and the time step of ion motion was 4 fs [32]. The Kohn–Sham energy function was minimized by the preconditioned conjugated-gradient method. The initial atomic configuration adopted was a random distribution of 90 atoms on the grid, which was constructed by dividing the supercell into 5 × 5 × 5 square segments. Then the system was heated up to 2000 K by rescaling the ionic velocities. After equilibration for 8 ps at this temperature, we gradually reduced the temperature to 1370 K. For other temperatures, we repeated this procedure and only changed the final temperature to 870 and 570 K, respectively. At each temperature, the physical properties of interest were obtained by averaging over 8 ps after an equilibration taking 4 ps. As is well-known, the increase of pressure corresponds to the increase of density in a normal liquid. We built a system with the density of 3.57 g/cm3 at 870 K from experiments data; after well equilibrating, we tentatively chose other densities 3.8, 4 and 4.5 g/cm3 at 870 K from low pressure to high pressure gradually. The physical parameters were also calculated by averaging over 8 ps after an equilibration of 4 ps.
Fig. 1. (a) The calculated pair correlation functions g(r) of liquid Se at three different temperatures; (b) The calculated pair correlation functions g(r) of liquid Se at 870 K with four different densities.
3. Results The pair correlation function plays a key role in the physics of liquids since it can be measurable directly, from which various properties can be estimated when coupled with appropriate theories. This quantity is used to describe three-dimensional spatial distribution of particles from one-dimensional point of view, and is proportional to the probability of finding a particle at a distance r from a reference particle. Using the atomic coordinates from simulations, the pair correlation functions under different temperatures and densities are displayed in Fig. 1. From which we can easily find two peaks, located at around 0.236 and 0.384 nm, in which the first peak presents a slight shift towards the right side with the rise of temperature and pressure, and the second peak moves to the larger r value at higher temperature but displays a different change tendency with increasing pressure. As expected, the heights of the two peaks decrease obviously under higher temperature and pressure, indicating that the correlations between the atoms become weaker. To gain a better understanding of the local structural properties around Se atoms, we have also calculated the contributions of the first, second…, sixth nearest-neighbor (NN) distances around a given atom as described in Ref. [33], and the results are shown in Fig. 2. It can be observed that the positions of the first and second NNs shift to the right, and the height of them decrease both with the increase of temperature and pressure. The positions and the heights of the third, fourth, fifth and sixth NNs present quite different temperature-dependence and pressure-dependence, and we will give a detailed discussion combined with g(r) in the following section.
In the physics of liquids, the structure factor S(Q) is also a fundamental quantity by which we can investigate the temperaturedependence and density-dependence of the atomic structure in the reciprocal space. S(Q) serves as a connection with experimental results and it can be obtained by two methods: one is calculated by Fourier transformation of pair correlation function g(r); the other is obtained directly from atomic coordinates considering Fouriertransformation calculation on a small simulation box may bring significant errors, SðQ Þ ¼
1 N
〈
→
→
∑ ∑ expð−i Q ⋅ r Þ i
j≠i
ij
〉
−NδQ ;0 ;
ð1Þ
where N is total number of atoms, rij is interatomic distance between atoms i and j. Calculated by these two methods, the total structure factors under different temperatures and densities are given in Fig. 3, in which the experimental results at three temperatures are also presented for a direct comparison. It can be found that S(Q) by two methods are in good agreement with each other on the whole; moreover, they correctly reproduce the peaks and temperaturedependence in experimental results. With rising temperature, one can observe that all peaks move towards the left side and grow broadened; particularly, the intensity of the first peak reduces significantly and it almost turns to be a shoulder at 1370 K. In Fig. 3, S(Q)s at different temperatures and pressures present relatively dramatically changes around the region from 17–27 nm − 1, corresponding to the intrachain bonds and the first two interchain bonds in liquid Se.
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
Fig. 2. Detailed contributions of the six neatest neighbors around a Se atom to g(r) under different temperatures and pressures.
875
876
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
coordination numbers present no significant variation and almost remain a constant of 2. We have also calculated the distributions of the number of neighbors in the intrachain bond when Rcutoff is chosen to be 0.286 nm, and the results are also presented in Fig. 4. We can find that the distribution of the number of neighbors is almost insensitive to the variation of temperature; however, two-neighbors-probability decreases apparently with the increase of pressure while threeneighbors-probability increases steadily. The bond-angle distribution function, g3(θ), is one of threebody distribution function and can be obtained from the atomic configuration in present simulations. The angle noted in is formed by a pair vectors drawn from a reference atom to any other atoms within a sphere of cutoff distance radius Rcutoff. The calculated bond-angle distributions of liquid Se under different conditions when Rcutoff = 0.286 nm are shown in Fig. 5. One can easily find a strong peak at around 103°, which consists with the previously-observed chain structure of crystalline modification [2]. With the increase of temperature and pressure, the height of 103°-peak decreases obviously while the position seems unmoved. As reported previously, the long-bond and short-bond induced by Peierls-type distortion in crystalline Se are 0.232 and 0.347 nm, respectively [13]. So we have calculated the bond-angle distributions with a longer cutoff distance Rcutoff = 0.350 nm, as displayed in Fig. 6. Compared to the bond-angle distributions with Rcutoff = 0.286 nm, the other two peaks can be found: the one peak, located at about 50°, is related to the close-packed configuration of atoms; the other peak, located at around 160°, is closed to the characteristics angles of Peierls-type distortion in trigonal phase Se, which
Fig. 3. (a) The calculated structure factors S(Q) of liquid Se at three different temperatures; (b) The calculated structure factors S(Q) of liquid Se at 870 K with four different densities, in which S(Q) calculated by Fourier transformation is plotted in solid line, calculated directly from coordinates is plotted in gray circle and line, compared with experimental data plotted in dash line.
Given the pair correlation function, the average coordination number can also be obtained by the integration of g(r) to its first minimum, R
2
N ¼ ∫o cutoff 4πr ρg ðrÞdr;
ð2Þ
where ρ0 is the number density, the cutoff distance Rcutoff is chosen to be the position of the first minimum of pair correlation functions, here is 0.286 nm. The temperature-dependence and density-dependence of N are shown in Table 1, from which one can find that the calculated Table 1 The calculated coordination numbers and viscosities under different temperatures and pressures. Temperature (K)
Density (g/cm3)
Coordination numbers
Calculated viscosities
570 870 870 870 870 1370
3.91 3.57 3.8 4.0 4.5 3.08
2.234 2.216 2.232 2.239 2.416 1.99
3.4754 2.7633 3.3274 3.9148 5.4 2.4208
Fig. 4. The distributions of the number of neighbors under different temperatures and pressures with Rcutoff = 0.286 nm.
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
877
Fig. 5. (a) The calculated bond angle distribution g3(θ) of liquid Se at three different temperatures with Rcutoff = 0.286 nm; (b) The calculated bond angle distribution g3(θ) of liquid Se at 870 K with four different densities with Rcutoff = 0.286 nm.
Fig. 6. (a) The calculated bond angle distribution g3(θ) of liquid Se at three different temperatures with Rcutoff = 0.350 nm; (b) The calculated bond angle distribution g3(θ) of liquid Se at 870 K with four different densities with Rcutoff = 0.350 nm.
may drop a hint of the Peierls-type distorted atomic configuration in the liquid state at high temperatures and pressure. As direct evidence of Peierls-type distorted local structures, the calculated angular limited bond–bond correlation functions are presented in Fig. 6. This function, P(r1, r2) is defined as the probability of finding an atom C at a distance r2 from an atom B, which is at a distance r1 far from the reference atom A. A constraint is placed → on the position of atom C. Namely, the angle θ formed by vectors BA and → BC is constraint in a small range. Here, θ is chosen between 150° and 170° in order to investigate the correlation of two bonds that form the angle around the third peak in g3(θ) with Rcutoff = 0.350 nm. For an undistorted structure, the maximum of the distribution should be on the diagonal (i.e. a maximum at r1 = r2) in the angular limited triplet correlation function; however, for a distorted structure, two maxima corresponding to the short-long and long-short correlations exist. In Fig. 6, the maxima are mainly located at two large regions centered at (0.240 nm, 0.335 nm) and (0.335 nm and 0.240 nm) at 570 K. As the temperature increases, the regions of short-long and long-short bonds correlations are reduced obviously. The locations of the maxima show no evident shift with the variations of temperature and pressure, which are very close to the characteristics bonds length of Peierls distortion in trigonal Se. Thermodynamic properties of liquid Se have also been studied in present work. As proposed by Iida et al., the viscosity of liquid Se can be estimated in terms of the pair distribution and average interatomic frequency,
the atom will stay in a oscillation state around a fixed coordination position, and a is the distance over which the transfer of momentum occurs [34]. It can be assumed that momentum interactions exist between the neighboring atoms. In Ref. [34], the nearest-neighbor distance is chosen as the minimum between the first and second peaks in g(r), which is used as the upper limit of integral in Eq. (3). According to the well-known Engkog theory of atomic transport, the collipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 sion frequency Г is given by Γ ¼ 4σ g ðσ Þρ πκ B T=m, in which m is the atomic mass and σ is the position of the first peak in g(r). Assuming that P(T) is inversely proportional to Г due to the fact that the atom will move to position once collision takes place, we can deduce Eq.(3) approximately as
2
a
4
η ¼ ð8π=9Þυ0 P ðT Þmρ ∫0 g ðr Þr dr;
ð3Þ
in which υ0 is a constant, corresponding to the frequency of oscillation when there is not net displacement, P(T) is the probability that
η∼
1 a 4 pffiffiffi ρ∫0 g ðrÞr dr: σ 2 g ðσ Þ T
ð4Þ
Calculated by Eq. (4), the relative viscosities under different temperatures and pressures are listed in Table 1. It can be observed that the calculated viscosities decrease significantly with rising temperature as reported in experiments, and they present an obvious increase under higher pressure as a normal liquid. The microscopic atomic structure is closely related to the electronic structure. Here, we have calculated the electronic density of state (DOS), i.e., the DOS for each atom species is decomposed into angular-momentum-resolved contributions. By projecting all wave functions in a sphere of radius R around atom i onto the spherical harmonic (l, m), we obtained the (l, m) angular momentum component of the atoms. The calculated DOS of liquid Se under different temperatures and densities are shown in Fig. 8. The most remarked characteristic of DOS is the appearance of a deep dip at the Fermi level EF, but the dip becomes shallower at high temperature and pressure.
878
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
Fig. 7. The angular limited bond-bond correlation functions under different temperatures and pressures.
Our simulations reproduce the semiconductor properties in the experiments, and the conductivity increases steadily with rising temperature. 4. Discussions Based on the experimental data, the atomic structural, dynamic and electronic properties have been investigated in present simulations. A good consistency can be found in these calculated results. The calculated structure factors in Fig. 3 correctly reproduce the positions and temperature-induced variations of peaks in experimental results, which indicate that our simulation could give an accurate picture of the local order and distance correlations in liquid Se under different circumstances. Likewise, although the calculated viscosities as displayed in Table 1 are not accurate results as obtained in experiments, they coincide well with the variation tendency as the temperature increases. Our simulations of DOS also reproduce the semiconductor properties in the experiments, and the conductivity increases steadily with rising temperature. The local atomic structure can be interpreted directly by the calculated pair correlation functions and the contributions of the first, second…, sixth nearest-neighbor (NN) distances, as displayed in Figs. 1 and 2, respectively. It can be found that the first and the second NNs contribute to the first peak of g(r), suggesting that they are the intrachain covalent bonds and their variations with rising temperature and pressure lead to the right shift of the first peak in g(r).
The third NN is the interchain bond, which locates at the position of the first trough of g(r) and could be considered as a result of the Peierls-type distortion, which will be discussed subsequently. The fourth, fifth and sixth NNs give rise to the second peak of g(r), which decrease obviously with rising temperature and shift towards the left side with increasing pressure. These variation tendencies fit well with those of the second peak of g(r), and suggest that the interchain bond is weakened significantly at high temperature and the bond length is reduced under high pressure. In other words, the interchain correlation is more sensitive to temperature and pressure compared with the intrachain correlation. As shown in Table 1, the calculated coordination numbers present no significant variation and almost remain a constant of 2, which clearly show that the chain structure with two neighbors reported previously still exists at high temperature and pressure. The distribution of the number of neighbors when Rcutoff = 0.286 nm has been presented in Fig. 4, from which we can observe that twofoldcoordinated atoms still prevail in liquid Se at high temperature but they could be disrupted under high density; meanwhile, threefoldcoordinated atoms increase gradually. In Ref. [13], the Peierls-type distorted local structure has been reported to exist in crystalline Se, with the short distance and the long distance being 0.232 and 0.347 nm, respectively, and this distortion could also be found in our calculated results. Compared to the bond-angle distributions with Rcutoff = 0.286 nm, the third peak of the bond-angle distributions with Rcutoff = 0.350 nm is located at
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
879
dependence are in good agreement with available experimental data. Our results suggest that the chain structure in crystalline modification still exists in liquid state under high temperature and pressure. The fraction of twofold-coordinated atoms decreases obviously under high pressure while it remains nearly invariable with rising temperature. The Peierls-type distorted structures in trigonal Se still reserve in liquid state even under high temperature and pressure. The calculated DOS displays an obvious dip at EF, and the dip becomes shallower with rising temperature.
Acknowledgements This work is supported by the National Science Foundation of China (Grant no. 61072032 and Grant no. 60910005), the China Postdoctoral Science Foundation (Grant no. 2011M501072) and by the Center for Computational Science, Hefei Institutes of Physics Sciences.
References
Fig. 8. (a) The calculated electronic density of state DOS of liquid Se at three different temperatures; (b) The calculated electronic density of state DOS of liquid Se at 870 K with four different densities.
around 160°, which is closed to the characteristics angles of the Peierls-type distortion in trigonal phase Se. The calculated angular limited bond-bond correlation functions as shown in Fig. 7 can provide us direct signature of the Peierls-type distorted structures: two maxima corresponding to the short-long and long-short correlations are found under different temperatures and pressures. The lengths of the bond are also very close to the characteristic bond lengths of Peierls-type distortion in trigonal Se, which indicate that the distortions still retain in the liquid state even at high temperature and pressure. 5. Conclusions In summary, the local atomic structural, thermodynamic and electronic properties of liquid Se under different temperatures and pressures were studied using ab initio molecular dynamics simulations. The structure factor, pair-correlation function, average coordination number, bond-angle distribution function, angular limited bondbond distribution function, viscosity and DOS were calculated. The calculated structure factors and their temperature/pressure
[1] C. Bichara, A. Pelegatti, Phys. Rev. B 49 (1964) 6581. [2] D. Hohl, R.O. Jones, Phys. Rev. B 43 (1991) 3856. [3] F. Kirchhoff, M.J. Gillan, J.M. Holender, G. Kresse, J. Hafner, J. Phys. Condens. Matter 8 (1996) 9353. [4] N.G. Almarza, E. Enciso, F.J. Bermejo, J. Chem. Phys. 99 (1993) 6876. [5] F. Shimojo, K. Hoshino, M. Watabe, Y. Zempo, J. Phys. Condens. Matter 10 (1998) 1199. [6] Y. Katayama, T. Mizutani, W. Utsumi, O. Shimomura, K. Tsuji, Phys. Status Solidi (b) 223 (2001) 401. [7] J. Donohue, The Structure of the Elements, Wiley, New York, 1974. [8] Y. Miyamoto, Jpn J. Appl. Phys. 19 (1980) 1813. [9] K. Tamura, S. Hosokawa, H. Endo, S. Yamasaki, H. Oyanagi, J. Phys. Soc. Jpn 55 (1956) 528. [10] M. Inui, K. Tamura, M. Yao, H. Endo, S. Hosokawa, H. Hoshino, J. Non-Cryst. Solids 117&118 (1990) 112. [11] J.D. Joannopoulos, in: E. Gerlach, P. Grosse (Eds.), The Physics of Selenium and Tellurium, Springer, Berlin, 1979. [12] G. Kresse, J. Furthmüller, J. Hafner, Phys. Rev. B 50 (1994) 13181. [13] J.P. Gaspard, A. Pellagatti, F. Marinelli, C. Bichara, Philos. Mag. B 77 (1998) 727. [14] J.C. Phillips, J. Non-Cryst. Solids 34 (1979) 153. [15] R. Strudel, J. Non-Cryst. Solids 83 (1986) 63. [16] J. Robertson, Philos. Mag. B 51 (1985) 183. [17] R. Bellissent, G. Tourand, J. Non-Cryst. Solids 35&36 (1980) 1221. [18] A. Axmann, W. Gissler, A. Kollmar, T. Springer, Discuss. Faraday Soc. 50 (1970) 74. [19] C. Bichara, A. Pellegatti, J.P. Gaspard, Phys. Rev. B 49 (1994) 6581. [20] S. Balasubramanian, K.V. Damodaran, K.J. Rao, Chem. Phys. 166 (1992) 131. [21] W.W. Warren Jr., R. Dupree, Phys. Rev. B 22 (1980) 2257. [22] H. Hoshino, et al., Phil. Mag. B 33 (1976) 255. [23] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [24] Y. Wang, J.P. Perdew, Phys. Rev. B 44 (1991) 13298. [25] K. Laasonen, A. Pasquarello, A.R. Car, R. Lee, D. Vanderbilt, Phys. Rev. B 47 (1993) 10142. [26] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [27] G. Kresse, J. Furthmüller, Comput. Mater. Sci. 6 (1996) 15. [28] G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169. [29] S. Nosé, J. Chem. Phys. 81 (1984) 511. [30] J. Thurn, J. Ruska, J. Non-Cryst. Solids 22 (1976) 331. [31] S. Hosokawa, K. Tamura, J. Non-Cryst. Solids 117&118 (1990) 52. [32] L. Verlet, Phys. Rev. 159 (1967) 98. [33] G. Zhao, Yu'e Zhao, Wang Yubing, Ji Changjian, Phys. Scr. 82 (2010) 035603. [34] T. Iida, R.I.L. Guthrie, The Physical Properties of Liquid Metals, Clarendon, Oxford, 1988.
Contents lists available at SciVerse ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/ locate/ jnoncrysol
Chain structure of liquid Se at high temperature and pressure investigated by ab initio molecular dynamics simulations Y.B. Wang a,⁎, W.S. Dong b, G. Zhao c, J.X. Ding d, S.H. Li a, Y.J. Ge a a
Institute of Intelligent Machines, Chinese Academy of Science, Post Office 1130, Hefei 230031, China School for the Gifted Young, University of Science and Technology of China, Hefei 230026, China Department of Physics, Ludong University, Hongqi Road, No. 186, Yantai 264205, China d New Star Research Institute of Applied Technology in Hefei City, Hefei 230031, China b c
a r t i c l e
i n f o
Article history: Received 25 July 2011 Received in revised form 30 December 2011 Available online 24 January 2012 Keywords: Liquid selenium; Molecular dynamics simulation; Chain structure
a b s t r a c t Using ab initio molecular dynamics simulations, the local atomic structural, dynamic and electronic properties of liquid selenium were studied under different temperatures and pressures. Compared with experimental data, the calculated structure factors and viscosities are acceptable on the whole. Our results indicate that the chain structure of crystalline selenium still exists in liquid state even at high temperature and pressure. The fraction of twofold-coordinated atoms decreases obviously under high pressure while it remains nearly invariable with the increase of temperature. The Peierls-type distorted structures in trigonal Se still reserve in liquid state even under high temperature and pressure. The calculated DOS displays an obvious dip at EF, and the dip becomes shallower with rising temperature. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The properties of element selenium have been extensively investigated for the past several years mainly due to the following two reasons [1–6]. Firstly, crystalline selenium is one of the elements with the most allotropes next to sulfur and phosphorus, and many allotropes of selenium are of great fascination from the viewpoints of the structure and properties [7]. There are five crystalline configurations at ambient conditions: t-Se, α-, β- and γ-monoclinic Se, and rhombohedral Se, all of which are twofold coordinated that bonded into infinite chains, only with a varying nearest-neighbor-distance from 2.32 to 2.37 Å [8–10]. Due to the valence-electronic configuration s2p4, the structure of liquid Se in the thermodynamic equilibrium state at room temperature and atmospheric pressure is the semiconductor trigonal phase, in which the helical chains of atoms run along the axes of the hexagonal cell [7,11]. This structure may be interpreted by the Peierls-type distortion of a simple cubic lattice: two s electrons are at deep level and do not participate in the covalent bonding; the four p electrons occupy two-thirds of the p-band, so that the simple cubic lattice is unstable against the Peierls-type distortion and the six neatest-neighbor bonds are divided into two short intrachain bonds and four long interchain bonds [12,13]. Secondly, selenium is a good candidate to form glass and amorphous structure with an accessible glass transition temperature (Tg ~ 310 K). As shown in his work previously, Phillips has concluded that covalent materials with an average coordination number N ≤ 2.4
⁎ Corresponding author. Tel.: + 86 1 551 5591195. E-mail address: [email protected] (Y.B. Wang). 0022-3093/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2012.01.001
are potential glass formers [14]. Meanwhile, amorphous selenium exhibits wide applications in photocopying areas in virtue of its photoconducting and semiconducting properties [15,16]. The explicit elaboration of liquid structure is also of great application importance in the control the glass-forming process of Se. Upon heating the infinite chain structure with twofold coordination is retained in the liquid state near the melting point (494 K), which is consistent with experimental-observed high viscosity and low diffusion in liquid Se [17]. However, the magnitude of viscosity decreases more than one order while the diffusion coefficient shows a marked increase as temperature is increased to above 1000 K, suggesting that the chain structure may be strongly destroyed [18]. This point of view has also been advanced by several computer simulations [19–20]. Moreover, it was reported that the chain structure was disrupted in liquid Se under high pressure and the average chain length was reduced sharply to about only ten atoms near the critical point (TC = 1860 K, PC = 380 bar) [21]. These temperature/pressure induced structural variations were also witnessed by a semiconductor–metal (SC-M) transition in the measurements of conductivity [22]. Although a large number of experiments and simulations have been performed on liquid Se, the related structural-change mechanisms are still not fully understood. In present work, we have performed AIMD simulations based on density functional theory (DFT) and the pseudopotential method to study the atomic structural and electronic properties of liquid Se under different temperatures and pressures, in order to illustrate the deeply-inside structural-change behaviors. The paper is organized as follows: in Section 2, we describe the methods of present simulations; the results of our simulations and the corresponding discussions
874
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
are presented in Sections 3 and 4 respectively; and finally a short summary is given in Section 5. 2. Computational methods Our simulations were performed with the framework of the density functional theory (DFT) [23], with the generalized gradient approximation (GGA) formulated by Perdew and Wang (PW91) to the exchange-correlation energy [24], and the valence electron-ion interaction was modeled by the ultrasoft pseudopotential (USPP) of the Vanderbilt type [25,26], as implemented in the Vienna Ab initio Simulation Package (VASP) [27,28]. The 90-atom system was put in a simple cubic box with periodic boundary conditions. Only the Γ point was used to sample the Brillouin zone of supercell. The electronic wave functions were expanded in the plane-wave basis set, with an energy cutoff of 160 eV. Our canonical ensemble simulations were performed at 570, 870 and 1370 K, with a Nosé thermostat to control temperature [29]. The densities as measured in Refs. [30] and [31] were adopted in our simulations. The Verlet algorithm was used to integrate Newton's equations of motion and the time step of ion motion was 4 fs [32]. The Kohn–Sham energy function was minimized by the preconditioned conjugated-gradient method. The initial atomic configuration adopted was a random distribution of 90 atoms on the grid, which was constructed by dividing the supercell into 5 × 5 × 5 square segments. Then the system was heated up to 2000 K by rescaling the ionic velocities. After equilibration for 8 ps at this temperature, we gradually reduced the temperature to 1370 K. For other temperatures, we repeated this procedure and only changed the final temperature to 870 and 570 K, respectively. At each temperature, the physical properties of interest were obtained by averaging over 8 ps after an equilibration taking 4 ps. As is well-known, the increase of pressure corresponds to the increase of density in a normal liquid. We built a system with the density of 3.57 g/cm3 at 870 K from experiments data; after well equilibrating, we tentatively chose other densities 3.8, 4 and 4.5 g/cm3 at 870 K from low pressure to high pressure gradually. The physical parameters were also calculated by averaging over 8 ps after an equilibration of 4 ps.
Fig. 1. (a) The calculated pair correlation functions g(r) of liquid Se at three different temperatures; (b) The calculated pair correlation functions g(r) of liquid Se at 870 K with four different densities.
3. Results The pair correlation function plays a key role in the physics of liquids since it can be measurable directly, from which various properties can be estimated when coupled with appropriate theories. This quantity is used to describe three-dimensional spatial distribution of particles from one-dimensional point of view, and is proportional to the probability of finding a particle at a distance r from a reference particle. Using the atomic coordinates from simulations, the pair correlation functions under different temperatures and densities are displayed in Fig. 1. From which we can easily find two peaks, located at around 0.236 and 0.384 nm, in which the first peak presents a slight shift towards the right side with the rise of temperature and pressure, and the second peak moves to the larger r value at higher temperature but displays a different change tendency with increasing pressure. As expected, the heights of the two peaks decrease obviously under higher temperature and pressure, indicating that the correlations between the atoms become weaker. To gain a better understanding of the local structural properties around Se atoms, we have also calculated the contributions of the first, second…, sixth nearest-neighbor (NN) distances around a given atom as described in Ref. [33], and the results are shown in Fig. 2. It can be observed that the positions of the first and second NNs shift to the right, and the height of them decrease both with the increase of temperature and pressure. The positions and the heights of the third, fourth, fifth and sixth NNs present quite different temperature-dependence and pressure-dependence, and we will give a detailed discussion combined with g(r) in the following section.
In the physics of liquids, the structure factor S(Q) is also a fundamental quantity by which we can investigate the temperaturedependence and density-dependence of the atomic structure in the reciprocal space. S(Q) serves as a connection with experimental results and it can be obtained by two methods: one is calculated by Fourier transformation of pair correlation function g(r); the other is obtained directly from atomic coordinates considering Fouriertransformation calculation on a small simulation box may bring significant errors, SðQ Þ ¼
1 N
〈
→
→
∑ ∑ expð−i Q ⋅ r Þ i
j≠i
ij
〉
−NδQ ;0 ;
ð1Þ
where N is total number of atoms, rij is interatomic distance between atoms i and j. Calculated by these two methods, the total structure factors under different temperatures and densities are given in Fig. 3, in which the experimental results at three temperatures are also presented for a direct comparison. It can be found that S(Q) by two methods are in good agreement with each other on the whole; moreover, they correctly reproduce the peaks and temperaturedependence in experimental results. With rising temperature, one can observe that all peaks move towards the left side and grow broadened; particularly, the intensity of the first peak reduces significantly and it almost turns to be a shoulder at 1370 K. In Fig. 3, S(Q)s at different temperatures and pressures present relatively dramatically changes around the region from 17–27 nm − 1, corresponding to the intrachain bonds and the first two interchain bonds in liquid Se.
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
Fig. 2. Detailed contributions of the six neatest neighbors around a Se atom to g(r) under different temperatures and pressures.
875
876
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
coordination numbers present no significant variation and almost remain a constant of 2. We have also calculated the distributions of the number of neighbors in the intrachain bond when Rcutoff is chosen to be 0.286 nm, and the results are also presented in Fig. 4. We can find that the distribution of the number of neighbors is almost insensitive to the variation of temperature; however, two-neighbors-probability decreases apparently with the increase of pressure while threeneighbors-probability increases steadily. The bond-angle distribution function, g3(θ), is one of threebody distribution function and can be obtained from the atomic configuration in present simulations. The angle noted in is formed by a pair vectors drawn from a reference atom to any other atoms within a sphere of cutoff distance radius Rcutoff. The calculated bond-angle distributions of liquid Se under different conditions when Rcutoff = 0.286 nm are shown in Fig. 5. One can easily find a strong peak at around 103°, which consists with the previously-observed chain structure of crystalline modification [2]. With the increase of temperature and pressure, the height of 103°-peak decreases obviously while the position seems unmoved. As reported previously, the long-bond and short-bond induced by Peierls-type distortion in crystalline Se are 0.232 and 0.347 nm, respectively [13]. So we have calculated the bond-angle distributions with a longer cutoff distance Rcutoff = 0.350 nm, as displayed in Fig. 6. Compared to the bond-angle distributions with Rcutoff = 0.286 nm, the other two peaks can be found: the one peak, located at about 50°, is related to the close-packed configuration of atoms; the other peak, located at around 160°, is closed to the characteristics angles of Peierls-type distortion in trigonal phase Se, which
Fig. 3. (a) The calculated structure factors S(Q) of liquid Se at three different temperatures; (b) The calculated structure factors S(Q) of liquid Se at 870 K with four different densities, in which S(Q) calculated by Fourier transformation is plotted in solid line, calculated directly from coordinates is plotted in gray circle and line, compared with experimental data plotted in dash line.
Given the pair correlation function, the average coordination number can also be obtained by the integration of g(r) to its first minimum, R
2
N ¼ ∫o cutoff 4πr ρg ðrÞdr;
ð2Þ
where ρ0 is the number density, the cutoff distance Rcutoff is chosen to be the position of the first minimum of pair correlation functions, here is 0.286 nm. The temperature-dependence and density-dependence of N are shown in Table 1, from which one can find that the calculated Table 1 The calculated coordination numbers and viscosities under different temperatures and pressures. Temperature (K)
Density (g/cm3)
Coordination numbers
Calculated viscosities
570 870 870 870 870 1370
3.91 3.57 3.8 4.0 4.5 3.08
2.234 2.216 2.232 2.239 2.416 1.99
3.4754 2.7633 3.3274 3.9148 5.4 2.4208
Fig. 4. The distributions of the number of neighbors under different temperatures and pressures with Rcutoff = 0.286 nm.
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
877
Fig. 5. (a) The calculated bond angle distribution g3(θ) of liquid Se at three different temperatures with Rcutoff = 0.286 nm; (b) The calculated bond angle distribution g3(θ) of liquid Se at 870 K with four different densities with Rcutoff = 0.286 nm.
Fig. 6. (a) The calculated bond angle distribution g3(θ) of liquid Se at three different temperatures with Rcutoff = 0.350 nm; (b) The calculated bond angle distribution g3(θ) of liquid Se at 870 K with four different densities with Rcutoff = 0.350 nm.
may drop a hint of the Peierls-type distorted atomic configuration in the liquid state at high temperatures and pressure. As direct evidence of Peierls-type distorted local structures, the calculated angular limited bond–bond correlation functions are presented in Fig. 6. This function, P(r1, r2) is defined as the probability of finding an atom C at a distance r2 from an atom B, which is at a distance r1 far from the reference atom A. A constraint is placed → on the position of atom C. Namely, the angle θ formed by vectors BA and → BC is constraint in a small range. Here, θ is chosen between 150° and 170° in order to investigate the correlation of two bonds that form the angle around the third peak in g3(θ) with Rcutoff = 0.350 nm. For an undistorted structure, the maximum of the distribution should be on the diagonal (i.e. a maximum at r1 = r2) in the angular limited triplet correlation function; however, for a distorted structure, two maxima corresponding to the short-long and long-short correlations exist. In Fig. 6, the maxima are mainly located at two large regions centered at (0.240 nm, 0.335 nm) and (0.335 nm and 0.240 nm) at 570 K. As the temperature increases, the regions of short-long and long-short bonds correlations are reduced obviously. The locations of the maxima show no evident shift with the variations of temperature and pressure, which are very close to the characteristics bonds length of Peierls distortion in trigonal Se. Thermodynamic properties of liquid Se have also been studied in present work. As proposed by Iida et al., the viscosity of liquid Se can be estimated in terms of the pair distribution and average interatomic frequency,
the atom will stay in a oscillation state around a fixed coordination position, and a is the distance over which the transfer of momentum occurs [34]. It can be assumed that momentum interactions exist between the neighboring atoms. In Ref. [34], the nearest-neighbor distance is chosen as the minimum between the first and second peaks in g(r), which is used as the upper limit of integral in Eq. (3). According to the well-known Engkog theory of atomic transport, the collipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 sion frequency Г is given by Γ ¼ 4σ g ðσ Þρ πκ B T=m, in which m is the atomic mass and σ is the position of the first peak in g(r). Assuming that P(T) is inversely proportional to Г due to the fact that the atom will move to position once collision takes place, we can deduce Eq.(3) approximately as
2
a
4
η ¼ ð8π=9Þυ0 P ðT Þmρ ∫0 g ðr Þr dr;
ð3Þ
in which υ0 is a constant, corresponding to the frequency of oscillation when there is not net displacement, P(T) is the probability that
η∼
1 a 4 pffiffiffi ρ∫0 g ðrÞr dr: σ 2 g ðσ Þ T
ð4Þ
Calculated by Eq. (4), the relative viscosities under different temperatures and pressures are listed in Table 1. It can be observed that the calculated viscosities decrease significantly with rising temperature as reported in experiments, and they present an obvious increase under higher pressure as a normal liquid. The microscopic atomic structure is closely related to the electronic structure. Here, we have calculated the electronic density of state (DOS), i.e., the DOS for each atom species is decomposed into angular-momentum-resolved contributions. By projecting all wave functions in a sphere of radius R around atom i onto the spherical harmonic (l, m), we obtained the (l, m) angular momentum component of the atoms. The calculated DOS of liquid Se under different temperatures and densities are shown in Fig. 8. The most remarked characteristic of DOS is the appearance of a deep dip at the Fermi level EF, but the dip becomes shallower at high temperature and pressure.
878
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
Fig. 7. The angular limited bond-bond correlation functions under different temperatures and pressures.
Our simulations reproduce the semiconductor properties in the experiments, and the conductivity increases steadily with rising temperature. 4. Discussions Based on the experimental data, the atomic structural, dynamic and electronic properties have been investigated in present simulations. A good consistency can be found in these calculated results. The calculated structure factors in Fig. 3 correctly reproduce the positions and temperature-induced variations of peaks in experimental results, which indicate that our simulation could give an accurate picture of the local order and distance correlations in liquid Se under different circumstances. Likewise, although the calculated viscosities as displayed in Table 1 are not accurate results as obtained in experiments, they coincide well with the variation tendency as the temperature increases. Our simulations of DOS also reproduce the semiconductor properties in the experiments, and the conductivity increases steadily with rising temperature. The local atomic structure can be interpreted directly by the calculated pair correlation functions and the contributions of the first, second…, sixth nearest-neighbor (NN) distances, as displayed in Figs. 1 and 2, respectively. It can be found that the first and the second NNs contribute to the first peak of g(r), suggesting that they are the intrachain covalent bonds and their variations with rising temperature and pressure lead to the right shift of the first peak in g(r).
The third NN is the interchain bond, which locates at the position of the first trough of g(r) and could be considered as a result of the Peierls-type distortion, which will be discussed subsequently. The fourth, fifth and sixth NNs give rise to the second peak of g(r), which decrease obviously with rising temperature and shift towards the left side with increasing pressure. These variation tendencies fit well with those of the second peak of g(r), and suggest that the interchain bond is weakened significantly at high temperature and the bond length is reduced under high pressure. In other words, the interchain correlation is more sensitive to temperature and pressure compared with the intrachain correlation. As shown in Table 1, the calculated coordination numbers present no significant variation and almost remain a constant of 2, which clearly show that the chain structure with two neighbors reported previously still exists at high temperature and pressure. The distribution of the number of neighbors when Rcutoff = 0.286 nm has been presented in Fig. 4, from which we can observe that twofoldcoordinated atoms still prevail in liquid Se at high temperature but they could be disrupted under high density; meanwhile, threefoldcoordinated atoms increase gradually. In Ref. [13], the Peierls-type distorted local structure has been reported to exist in crystalline Se, with the short distance and the long distance being 0.232 and 0.347 nm, respectively, and this distortion could also be found in our calculated results. Compared to the bond-angle distributions with Rcutoff = 0.286 nm, the third peak of the bond-angle distributions with Rcutoff = 0.350 nm is located at
Y.B. Wang et al. / Journal of Non-Crystalline Solids 358 (2012) 873–879
879
dependence are in good agreement with available experimental data. Our results suggest that the chain structure in crystalline modification still exists in liquid state under high temperature and pressure. The fraction of twofold-coordinated atoms decreases obviously under high pressure while it remains nearly invariable with rising temperature. The Peierls-type distorted structures in trigonal Se still reserve in liquid state even under high temperature and pressure. The calculated DOS displays an obvious dip at EF, and the dip becomes shallower with rising temperature.
Acknowledgements This work is supported by the National Science Foundation of China (Grant no. 61072032 and Grant no. 60910005), the China Postdoctoral Science Foundation (Grant no. 2011M501072) and by the Center for Computational Science, Hefei Institutes of Physics Sciences.
References
Fig. 8. (a) The calculated electronic density of state DOS of liquid Se at three different temperatures; (b) The calculated electronic density of state DOS of liquid Se at 870 K with four different densities.
around 160°, which is closed to the characteristics angles of the Peierls-type distortion in trigonal phase Se. The calculated angular limited bond-bond correlation functions as shown in Fig. 7 can provide us direct signature of the Peierls-type distorted structures: two maxima corresponding to the short-long and long-short correlations are found under different temperatures and pressures. The lengths of the bond are also very close to the characteristic bond lengths of Peierls-type distortion in trigonal Se, which indicate that the distortions still retain in the liquid state even at high temperature and pressure. 5. Conclusions In summary, the local atomic structural, thermodynamic and electronic properties of liquid Se under different temperatures and pressures were studied using ab initio molecular dynamics simulations. The structure factor, pair-correlation function, average coordination number, bond-angle distribution function, angular limited bondbond distribution function, viscosity and DOS were calculated. The calculated structure factors and their temperature/pressure
[1] C. Bichara, A. Pelegatti, Phys. Rev. B 49 (1964) 6581. [2] D. Hohl, R.O. Jones, Phys. Rev. B 43 (1991) 3856. [3] F. Kirchhoff, M.J. Gillan, J.M. Holender, G. Kresse, J. Hafner, J. Phys. Condens. Matter 8 (1996) 9353. [4] N.G. Almarza, E. Enciso, F.J. Bermejo, J. Chem. Phys. 99 (1993) 6876. [5] F. Shimojo, K. Hoshino, M. Watabe, Y. Zempo, J. Phys. Condens. Matter 10 (1998) 1199. [6] Y. Katayama, T. Mizutani, W. Utsumi, O. Shimomura, K. Tsuji, Phys. Status Solidi (b) 223 (2001) 401. [7] J. Donohue, The Structure of the Elements, Wiley, New York, 1974. [8] Y. Miyamoto, Jpn J. Appl. Phys. 19 (1980) 1813. [9] K. Tamura, S. Hosokawa, H. Endo, S. Yamasaki, H. Oyanagi, J. Phys. Soc. Jpn 55 (1956) 528. [10] M. Inui, K. Tamura, M. Yao, H. Endo, S. Hosokawa, H. Hoshino, J. Non-Cryst. Solids 117&118 (1990) 112. [11] J.D. Joannopoulos, in: E. Gerlach, P. Grosse (Eds.), The Physics of Selenium and Tellurium, Springer, Berlin, 1979. [12] G. Kresse, J. Furthmüller, J. Hafner, Phys. Rev. B 50 (1994) 13181. [13] J.P. Gaspard, A. Pellagatti, F. Marinelli, C. Bichara, Philos. Mag. B 77 (1998) 727. [14] J.C. Phillips, J. Non-Cryst. Solids 34 (1979) 153. [15] R. Strudel, J. Non-Cryst. Solids 83 (1986) 63. [16] J. Robertson, Philos. Mag. B 51 (1985) 183. [17] R. Bellissent, G. Tourand, J. Non-Cryst. Solids 35&36 (1980) 1221. [18] A. Axmann, W. Gissler, A. Kollmar, T. Springer, Discuss. Faraday Soc. 50 (1970) 74. [19] C. Bichara, A. Pellegatti, J.P. Gaspard, Phys. Rev. B 49 (1994) 6581. [20] S. Balasubramanian, K.V. Damodaran, K.J. Rao, Chem. Phys. 166 (1992) 131. [21] W.W. Warren Jr., R. Dupree, Phys. Rev. B 22 (1980) 2257. [22] H. Hoshino, et al., Phil. Mag. B 33 (1976) 255. [23] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [24] Y. Wang, J.P. Perdew, Phys. Rev. B 44 (1991) 13298. [25] K. Laasonen, A. Pasquarello, A.R. Car, R. Lee, D. Vanderbilt, Phys. Rev. B 47 (1993) 10142. [26] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [27] G. Kresse, J. Furthmüller, Comput. Mater. Sci. 6 (1996) 15. [28] G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169. [29] S. Nosé, J. Chem. Phys. 81 (1984) 511. [30] J. Thurn, J. Ruska, J. Non-Cryst. Solids 22 (1976) 331. [31] S. Hosokawa, K. Tamura, J. Non-Cryst. Solids 117&118 (1990) 52. [32] L. Verlet, Phys. Rev. 159 (1967) 98. [33] G. Zhao, Yu'e Zhao, Wang Yubing, Ji Changjian, Phys. Scr. 82 (2010) 035603. [34] T. Iida, R.I.L. Guthrie, The Physical Properties of Liquid Metals, Clarendon, Oxford, 1988.