J. Phys. Chem. B 2006, 110, 4111-4114
4111
Molecular Dynamics Study of the Thermal Entropy in Mixed Zinc Chalcogenides Karl Jug,* Nisanth N. Nair, and Igor P. Gloriozov† Theoretische Chemie, UniVersita¨t HannoVer, Am Kleinen Felde 30, 30167 HannoVer, Germany ReceiVed: December 8, 2005; In Final Form: January 9, 2006
Molecular dynamics (MD) calculations were performed to determine the vibrational contribution to the entropy of mixing and its importance for the mixing of ZnO/ZnS and ZnS/Zn3P2. These systems were modeled by cyclic clusters Zn48O48, Zn32S32, and Zn48P32. The mixed cyclic clusters considered were Zn48O47S, Zn32S31O, Zn33S30P2, and Zn47S2P30. For each of the clusters, the entropy was calculated in the range of the experimental temperature of the mixing process. The convergence of the entropy with respect to the number of MD steps was studied. Finally, the thermal part of the entropy of mixing was determined, and its dependence on the number of MD steps was investigated. It was found that the thermal entropy is important for the Gibbs free energy of mixing near the miscibility gaps.
1. Introduction
2. Property Calculations with Molecular Dynamics
Zinc chalcogenides are of technological importance because of their suitability for optoelectronic devices.1 Their band gaps classify them as semiconductors in the blue (ZnSe) and ultraviolet (ZnO, ZnS) light ranges. Modification of the band gap could provide taylor-made lasers in a desired range. Recently, anion doping has received much attention. Mixing of ZnS and ZnSe occurs without a miscibility gap, and the change of the band gap dependent on the mixing parameter has been extensively studied both experimentally2-4 and theoretically.5 It is more difficult to study the mixing of ZnO and ZnS or ZnS and Zn3P2, because a miscibility gap exists, which means that the thermodynamic stability of the solid solution exists only in a small range. With pulsed-laser deposition6 or radio frequency reactive sputtering,7 thin layers of type ZnO/S have been produced. But these systems appear to be only metastable. Locmelis and Binnewies8-10 generated solid solutions of mixed phases ZnS/P and Zn3P2/S by chemical transport at 900 °C. This technique leads to thermodynamically stable states. In our previous work, we studied the miscibility of ZnS/ZnSe and ZnO/ ZnS11 and found continuous mixing in the first case but no mixing in the second case at 1000 K. Mixing was predicted only at higher temperatures12 in contrast to the experiments. There could be two reasons for this discrepancy: (1) the small size of the cyclic cluster used for the simulation of the solids and (2) the neglect of the vibrational contribution to the entropy of mixing. The latter term lowers the Gibbs free energy of mixing. Initial molecular dynamics (MD) studies indicated a possible influence of this thermal contribution.13 A subsequent theoretical study on the mixing of ZnS and Zn3P2 14 revealed the same problem. We therefore decided to investigate this case by MD techniques as well and to give a detailed account of the influence of the vibrational contribution on the thermodynamic stability of these mixed compounds. In section 2, we present the formalism and in the following sections the results of the calculations.
Usually, the thermodynamic properties are calculated from the fluctuations of total energy.15 The molar heat capacity at constant volume CV is given by
* To whom correspondence should be addressed. E-mail: Jugthc@ mbox.theochem.uni-hannover.de. † Visiting scientist, permanent address: Department of Chemistry, Moscow State University, Moscow 119992, Russia.
CV )
(∂E∂T)
(1)
NV
CV can be estimated by numerical differentiation of eq 1 after performing a number of canonical ensemble simulations at different temperatures. The following analytical expression can also be used to obtain CV from a single NVT simulation at a certain temperature T.15
CV )
1 (〈E2〉 - 〈E〉2) 2 kBT
(2)
As CV is obtained from fluctuations rather than from a quantity directly obtained from the ensemble average, it is difficult to obtain accurate values for CV from an MD simulation. The entropy S is related to other thermodynamic properties as
dS )
µ dE P + dV - dN T T T
(3)
Here, P is the pressure, µ ) G/N is the chemical potential, and G is the Gibbs free energy. For a canonical ensemble, we have
∆S ) S(T2) - S(T1) )
∫TT 1
2
1 ∂E dT ) T ∂T NV
( )
∫TT 1
2
CV dT T
(4)
The integration in eq 4 is evaluated numerically by performing several canonical ensemble simulations at different temperatures between T1 and T2. A quantum-classical approach can be used for the calculation of thermodynamic properties of crystals. For crystals, the vibrational modes can be approximately considered as harmonic.16 The system is then viewed as a set of 3N harmonic oscillators.17,18 The total canonical partition function Q can be expressed in terms of the partition functions of the individual
10.1021/jp0571527 CCC: $33.50 © 2006 American Chemical Society Published on Web 02/08/2006
4112 J. Phys. Chem. B, Vol. 110, No. 9, 2006
Jug et al.
modes qj as
be used for crystals in a wide range of temperatures, while for molecular systems it is applicable only at very low temperatures. 3N
Q)
qj ∏ j)1
3. Mixing of Zinc Oxide and Zinc Sulfide Mixing of two different zinc chalcogenides ZnA and ZnB would result in a solid solution ZnA1-xBx.
Therefore
(1 - x)ZnA + xZnB f ZnA1-xBx
3N
ln Q )
ln qj ∑ j)1
(5)
If the normal frequencies are continuously distributed, then eq 5 can be written as
ln Q )
∫0∞ dν G(ν) ln q(ν)
(6)
Here, G(ν) is the distribution of the vibrational normal modes of the system. The effective vibrational intensity at frequency ν can be calculated as
G(ν) )
2
N
3
∑ ∑ MJ gKJ (ν)
(7)
kBT J)1 K)1
The spectral density gKJ (ν) of an atom J with mass MJ, in the Kth coordinate (K ) x, y, and z in the Cartesian coordinate system), is estimated from the square of the Fourier transform of the velocities as
gKJ (ν) ) lim
|
∫-τR˙ JK(t) exp(-i2πνt) dt| τ
(8)
Partition function q(ν) is given by
q(ν) )
exp(-βhν/2) 1 - exp(-βhν)
where β )
1 kBT
The entropy S can then be expressed as
S ) kB ln Q + β-1
(∂ ln∂TQ)
N,V
) kB
∫0∞ dν G(ν)Ws(ν)
(9)
where
Ws(ν) )
βhν - ln[1 - exp(-βhν)] exp(βhν) - 1
Similarly, the heat capacity at constant volume CV is given by
C V ) kB
∫0∞ dν G(ν)Wc(ν)
The stability of a solid solution is determined by its free energy of mixing ∆MG.
∆MG ) ∆MH - T∆MS
(12)
Here, ∆MH is the enthalpy of mixing, ∆MS is the entropy of mixing, and T is the temperature. The thermodynamic stability of the mixed compounds requires a negative ∆MG. The entropy of mixing ∆MS consists of configurational entropy ∆MSconf, which is related to the number of ways to arrange the distinguishable anions in the lattice, and thermal entropy ∆MSther, due to the vibrational motion of the atoms of the solid.
∆MS ) ∆MSconf + ∆MSther
(13)
By assuming that the mixing in reaction 11 leads to an ideal solution, ∆MSconf is estimated as
∆MSconf ) -R(1 - x)ln(1 - x) - Rx ln(x)
(14)
where R is the gas constant. ∆MSther is evaluated for reaction 11 as
2
2τ
τf∞
(11)
(10)
with
exp(βhν) Wc(ν) ) (βhν)2 (exp(βhν) - 1)2 In this approach, the entropy is obtained from a single NVT simulation, while the previous method which uses the thermodynamic integration eq 4 needs several NVT simulations. The quantum-classical method is applicable only in the temperature range in which the harmonic approximation is valid. It can thus
∆MSther ) S(ZnA1-xBx) - (1 - x)S(ZnA) - xS(ZnB)
(15)
The entropy S is temperature dependent and can be evaluated using the MD techniques decribed in the previous section. Solid solutions ZnS1-xSex produced by the mixing of ZnS and ZnSe were found to be stable for all possible values of x (0 < x < 1).19,20 The stability of ZnS1-xSex is due to the similar size of the anions. But due to the difference in the sizes of oxygen and sulfur, a miscibility gap is expected for ZnO and ZnS. Yoo et al.6 produced ZnO1-xSx with a very low concentration of S (x ≈ 0.03) at about 1000 K using laser ablation techniques. Locmelis et al.21 prepared ZnO/ZnS mixed crystals with low sulfur content (x ≈ 0.01) at about 1200 K. In previous work, the cyclic cluster model (CCM)22,23 implemented in MSINDO24 was used for the simulation of solids. Thus, the construction of extended periodic clusters (equivalent to supercells) from minimal clusters (equivalent to unit cells) is straightforward, facilitating the study of convergence with respect to system size. The use of a localized atomic basis for valence electrons together with pseudopotentials for the inner electrons circumvents the need to perform Brillouin zone (k-point) sampling, which would be required if a planewave basis set were employed. In addition, MSINDO uses only a few adjustable parameters based on the properties of the atoms and small molecules.11,14 With MSINDO-CCM, the miscibility of ZnS with ZnSe and ZnO was investigated.11 It was found that ZnS and ZnSe can mix for all possible values of x between 0 and 1. But a miscibility gap was observed for ZnS and ZnO at 1000 K. The latter result is in disagreement with the experimental observation that a small concentration of ZnS can mix with ZnO.21 In these calculations, vibrational contributions to the entropy were not included (∆MSther ) 0) and the enthalpy of mixing
Thermal Entropy in Mixed Zinc Chalcogenides
J. Phys. Chem. B, Vol. 110, No. 9, 2006 4113
Figure 1. Convergence of entropy S at 1000 K for ZnO with respect to the number of MD steps.
∆MH was approximated as the energy of mixing, ∆ME0.
∆MH ≈ ∆ME0
(16)
For reaction 11, ∆ME0 is given by
∆ME0 ) E0(ZnA1-xBx) - (1 - x)E0(ZnA) - xE0(ZnB) (17) Here, E0 is the total energy obtained from MSINDO. For the simulation of the ZnO and ZnS bulk three-dimensional (3D) cyclic clusters, Zn48O48 and Zn32S32 were chosen, respectively. In the previous study,11 these clusters were found to give converged bulk properties. MSINDO parameters were taken from there for the corresponding MD study in this section. The substituted clusters considered in the present work are Zn48O47S and Zn32S31O. They correspond to ZnO0.98S0.02 and ZnO0.03S0.97 systems, with x equal to 0.02 and 0.97, respectively. As the mixing of ZnO and ZnS was observed at about 1000 K in the experiments, calculations were also performed at 1000 K. To estimate the entropy, a canonical ensemble simulation was performed using the Nose´-Hoover chain thermostat.25 The time step of the simulation was 0.5 fs. A small time step increases the accuracy of the estimated spectral density, obtained by performing the numerical integration in eq 8. This was found necessary for an accurate determination of the entropy. The frequency of the thermostat was 750 cm-1, which is close to the maximum vibrational frequency of the bulk. The entropy of ZnO was estimated by using the 3D cyclic cluster Zn48O48. The system was equilibrated for 2000 steps at 1000 K. This was followed by an MD run for 12 000 steps (6 ps) from which the vibrational density of states of ZnO was obtained by the Fourier transform of the mass weighted atomic velocities using eqs 7 and 8. The entropy was evaluated using the vibrational density of states according to eq 9. Convergence of the entropy with respect to number of MD steps is shown in Figure 1. After 12 000 steps of the simulation, fluctuations in entropy were below
