Molecular-dynamics Simulation of Structural and Thermodynamic Properties of Ternary Alloy ZnSxSe 1x

Molecular dynamics simulation is being widely applied in condensed matter physics, rather Monte Carlo simulation, whenever the dynamic properties are of much interest. Particulary of equilibrium and non-equilibrium phenomena, transport phenomena, diffusivity, specific heat and phase transitions etc. Therefore, we have used this method, in order to obtain more informations about the structural, elastic and thermodynamic properties of ZnS0.25Se0.75, ZnS0.5Se0.5 and ZnS0.75Se0.25. We have also predicted the disorder effect (bowing parameter). We have used molecular simulation coupled with the empirical interatomic potentials on the Tersoff which include the three-body interatomic interactions. We have proposed empirical potentials of Tersoff parameters of ZnSxSe1-x for x = 0.25, x = 0.50 and x = 0.75. We have tested the validity of our fitting of empirical interatomic potentials parameters by calculating the physical properties. Overall, the results are in good agreement with other calculations.


Introduction
In recent years, computer simulations have become an increasingly powerful tool for studying materials properties and their behaviour.While first principles quantum methods generally give the most accurate results, they cannot be applied to problems which require much larger system sizes or longer simulations.For these investigations, empirical many-body potentials can be used to provide useful first results The wide band gap II-VI semiconductors have been regarded as promising materials for the fabrication of visible-light-emitting devisces for decades.ZnS and ZnSe are the prototype II-VI semiconductors and their cubic phase, which occurs naturally as a mineral, has been called the zinc-blende structure (M. A. Haase, et al., 2001).On the theoretical side,many authors have studied these semiconductors to explore their electronic properties, using various computational methods (W.Xie, et al., 1993).ZnSxSe1-x is an interesting compound semiconductor for optoelectronic applications (M. A. Haase, et al., 2001, W. Xie, et al., 1993, M. J. Kastner, et al., 1996).The ternary compound has many advantages over binary compounds: ZnS and ZnSe, because the lattice constant, the bulk modulus, band gap and optical properties can be varied by changing the concentration (H.Jeon, et al., 1991).
However, no elastic and thermodynamic calculations have been performed on ZnS and ZnSe and their alloy, using a three-body potential coupled with Molecular Dynamics simulation, to our knowledge.
The mixed crystals were prepared and analysed as a complete solid solution and a core-shell type of crystal by the X-ray diffraction method.X-ray diffraction patterns of these alloy systems indicated a cubic zinc-blende structure and showed that the lattice parameter changed linearly with alloy concentration without changes in the crystal structure (Chin-Tsar Hsu, J., 1998).
In this paper, we present the results of molecular dynamics simulation of ZnSxSe1-x for x  0.25, x = 0.5 and x = 0.75; which are compared with experimental results (James E. Bernard and Alex Zunger, 1987).
Based on the Tersoff potential, molecular dynamics simulations are investigated in order to study the miscibility of ZnSxSe1-x and to calculate the structural properties of ZnS 0.25 Se 0.75 , ZnS 0.5 Se 0.5 and ZnS 0.75 Se 0.25 which are compared with experimental results (James E. Bernard and Alex Zunger, 1987).Various elastic and thermodynamic properties are also predicted, such as C 11 , C 12 , thermal expansion coefficient, heat capacity, Debye temperature and melting temperature.
The rest of the paper is organised as follows: in Section 2, we give a brief description of the model used, and calculation method.Then, the results of our calculation are presented and discussed in Section 3. Finally, a conclusion is given in Section 4.

Details of calculations
A molecular dynamics calculation is a dynamic study of a many-body system in which the equations of motion are solved explicitly.In the MD calculation one solves the equations motion numerically and from the phase space trajectories one can follow the motions of the atoms in an isolated system (constant N, V and E) or (N, V and T).The potentials of interactions of the atoms, and hence the force, are necessary input parameters; This trajectory contains all of the microscopic information relevant to the system.We present an atomistic simulation based on a three-body potential, and we mainly use the Tersoff model to calculate dynamic quantities and, to predict the different thermodynamical properties.
Among the many empirical model potentials that have been developed for tetrahedral semiconductors, that of Tersoff has been the most successful in that it reproduces many of the elemental semiconductors properties particularly for silicon ( J. Tersoff, 1988, Rev. B. 38), and carbon ( J. Tersoff, 1988, Phys. Rev. Lett. 61).
Another form is developed for multicomponent systems (J.Tersoff.1989, Phys. Rev. B. 39) to treat mixtures of these elements.This potential developed by Tersoff, are based on the concept of bond order.The strength of a bond between two atoms is not constant, but depends on the local environnement (A.Laref, et al., 1999).
Because of the crucial role of bond order and its dependence upon local geometry, it seems attractive to include an environment-dependent bond order explicity into the potential in the following way.The interatomic potential is taken to have the form (J. Tersoff.1989, Phys. Rev. B. 39): Here E is the total energy of the system, which is decomposed for convenience into a site energy Ei and a bond energy Vij.The indices i and j run over the atoms of the system, and r ij is the distance from atom i to atom j.
The function f R represents a repulsive pair potential and f A represents a attractive pair potential associated with bonding.f C is a smooth cutoff function, to limit the range of the potential, since for many applications short-ranged functions permit a tremendous reduction in computational effort.
where b ij is the many-body order parameter describing how the bond-formation energy is affected by local atomic arrangement due to the presence of other neighbouring atoms ( the k atoms).It is many-body function of the positions of atoms i, j and k.
where ζ ij is called the effective coordination number and g( ijk ) is a function of the angle between r ij and r ik that has been fitted to stabilize the tetrahedral structure.
We note that  3 and  are set equal to zero (J.Tersoff., 1998), and  is a parameter which strengthens or weakens the heteropolar bonds, relative to the value obtained by simple interpolation.
On the basis of the transferability of the Tersoff potential (W.Sekkal, et al., 1998, W. Sekkal, et al., 19981), all the parameters for Zn S 0.25 Se 0.75 , and ZnS 0.75 Se 0.25 have been kept constant except the cutt-off R and , which are adjusted to ensure the stability of each structure.
All parameters are listed in Table 1 Using Molecular dynamic calculations, the interactions are described for different concentrations by the the Tersoff potential.In our simulation, we assume that the system is treated as two components (ZnS and ZnSe).Indeed, ZnS and ZnSe are considered as an equivalent one-component system.Initially, ZnS and ZnSe are always located at the nodes of a zincblende structure with interaction via the Tersoff potential.The integral energy of the system is a sum of interactions between pairs of ZnS-ZnS, ZnSe-ZnSe, ZnS-ZnSe.
In order to obtain more information about the structural, elastic and thermodynamic properties of ZnS x Se 1-x we run molecular dynamics simulation using the Tersoff potential with 216 molecules.This simulation is carried out within the canonical NVT-ensemble where the control is investigated using Andersen's method (H.C. Andersen, 1980).
The MD cell is formed in a cube of side L with 3 x3 x 3 FCC unit cells, where 216 molecules with 108 ZnSe and 108 ZnS (molecules) are included.The periodic boundary conditions are applied.The molecular dynamics routine is based on a fifth-order gear-predictor-corrector algorithm (G.W. Gear, 1971) of the Newtonian equations of motion using a neighbor list technique with a time step h  5.02 fs.After 100 ps, different properties are computed along the trajectory of the system in phase space.

Results and discussion
We display, in Figure .1 the pair distribution function g(r) for ZnS 0.25 Se 0.75 , ZnS 0.5 Se 0.5 and ZnS 0.75 Se 0.25 alloy.The results for the fourth peaks listed in Table 2 are in good agreement with the values in Ref. (W.Sekkal, et al., 2000, K. Endo, et al., 1993).Thus, we confirm the stability of the structure during the simulation.In order to calculate the structural and elastic properties of ZnS x Se 1-x , we plot in Figure .2 the variation of the cohesive energy with the volume for ZnS 0.25 Se 0.75 , ZnS 0.5 Se 0.5 and ZnS 0.75 Se 0.25 in the zincblende structure.The curves are fitted to the Murnaghan equation of state (F.D. Murnaghan, 1944) from which the equilibrium cohesive energy, the lattice parameter, the bulk modulus, and its derivative are obtained.
From the results listed in Table 3, we notice the agreement of the lattice parameter with experiments (James E. Bernard and Alex Zunger, 1987), the accuracy is between 0.1 and 0.8.Our results of the bulk modulus are in good agreement with other calculations based on the Cohen's relation (M.L. Cohen, 1985): where d is the nearest-neighbor distance in A°, and B is the bulk modulus given in Mbar.Using this equation, we found good values of the bulk modulus for ZnS 0.25 Se 0.75 (0.78 Mbar), ZnS 0.5 Se 0.5 (0.82 Mbar), ZnS 0.75 Se 0.25 (0.85 Mbar) with an accuracy about 17 , 17  and 15  respectively As we see in Figure .3,the variation of B with concentration is quadratic and is given by the following equation: The quadratic term stands for the bowing parameter.This bowing is known to be an effect caused by disorder.Indeed, the disorder effect contributes in two different ways: i) The chemical disorder (contribution due to the compositional disorder in the anionic sublattice). ii) The volume deformation effect (contribution due to the compression and the dilatation of the two binary alloys into the ternary alloy volumes).
According to van Vechten and Bergstresser (J. A. Van Vechten and J. K. Bergstresser, 1970), the disorder contribution (which is related to the difference of electronegativities of the alloyed atoms) plays a dominant role in determining the bowing parameter.
In ZnS x Se 1-x, the difference of electronegativities of S and Se atoms is equal to 0.1.Moreover, our results show that the bowing parameter is equal to 0.068.This value is close to the value found for Cu Ag x I 1-x (0.08) (W.Sekkal, et al., 2000).In Figure .4,we plot the lattice parameter for ZnS x Se 1-x solid solution.We notice that for all concentrations (x = 0.25, x = 0.5, x = 0.75), the calculated lattice parameters are in good agreement with experimental results (James E. Bernard and Alex Zunger, 1987) for which the variation is linear and follow Vegard's law (Vegard L. Z, 1921) very closely.In conclusion, we can say that there is an ideal mixing between ZnSe and ZnS.
To compute elastic constants, we will follow the treatment established in Ref.( M.J. Mehl, 1993) which we shall outline here briefly.To compute C 11 -C 12 , we use a volume-conserving strain matrix: which yields an energy change E =  2 (C 11 -C 12 ) V/3 + O( 3 ) where V is the volume.We compute the total energy at several values of  and drop the terms of order  3 and higher.Plotting E versus  2 yields a straight line with slope = (C 11 -C 12 ) V/3.
Based on the fact that our results for C 11 and C 12 agree well with experimental results for ZnS and ZnSe (see Table 3), we have extended our calculations to predict these quantities for ZnS 0.25 Se 0.75 , ZnS 0.5 Se 0.5 and ZnS 0.75 Se 0.25 .
We are interested now in calculating the thermodynamic properties of ZnS x Se 1-x for x = 0.25, x = 0.5, and x = 0.75.
In MD simulation, the linear thermal expansion coefficient  can be computed directly from the definition: where a is the lattice parameter.Therefore, we consider the temperature variation of the lattice constant at zero pressure (see Figure .5)A molecular dynamics simulation is performed in an NVT ensemble at each temperature to equilibrate the system and then to determine the corresponding zero pressure lattice constant.
From the slope of the total energy versus temperature curve (see Figure .6),we can estimate the specific heat of the system (A.Garcia and C. Penland, 1991) according to the following equation: Our calculations can be used to determine an approximate melting.Fine et al (M. E. Fine, et al., 1984) noticed empirically that the melting temperature and elastic constants of cubic materials are related by the approximate expression: Debye temperature is calculated from the bulk modulus using the empirical relation taken from Ref. (W.Sekkal, et al., 1998): All the predicted thermodynamic properties are listed in Table 3.
Based on the fact that the values of  D and T f agree well with experimental calculations for ZnSe and ZnS, we have predicted these quantities for ZnS x Se 1-x at x= 0.25, x = 0.5, and x = 0.75 using Equations ( 18) and ( 19).According to our results, we see that the Debye temperature increases with concentration (x) according to the following equation (see Figure.The specific heat for ZnS 0.25 Se 0.75 , ZnS 0.5 Se 0.5 and ZnS 0.75 Se 0.25 alloy is around 24.6 J K -1 .mol - .This constant value corresponds to the classical Dulong-Petit result ( 24 J K -1 .mol - ) which is obtained at high temperatures for all solids (W .A. Harisson, 1980).

Conclusion
We have presented the results of molecular dynamics calculations, based on the well-tested Tersoff potential.Structural, elastic and thermodynamic properties are performed at 0.25, 0.5 and 0.75.We found also, that the lattice parameter and the bulk modulus varies linearly with alloy concentration and the results are in good agreement with experiments.However, the apparition of a bowing in the variation of the bulk modulus with concentration can be due to the volume deformation effect and the structural relaxation has a minor effect.We have predicted elastic and thermodynamic properties such C 11 , C 12 , heat capacity, melting point, Debye temperature and linear thermal expansion coefficient for ZnS 0.25 Se 0.75 , ZnS 0.5 Se 0.5 and ZnS 0.75.