14290
J. Phys. Chem. B 2007, 111, 14290-14294
Vapor Pressure and Sublimation Rate of Molecular Crystals: Role of Internal Degrees of Freedom A. Maiti,* L. A. Zepeda-Ruiz, R. H. Gee, and A. K. Burnham Lawrence LiVermore National Laboratory, UniVersity of California, LiVermore California 94551 ReceiVed: July 30, 2007; In Final Form: October 8, 2007
It is a common practice to approximate the desorption rate of atoms from crystal surfaces with an expression of the form νeff exp(-∆E/kBT), where ∆E is an activation barrier to desorb and νeff is an effective vibrational frequency ∼1012 s-1. For molecular solids, however, such an approximation can lead to a many orders of magnitude underestimation of vapor pressure and sublimation rates due to neglected contributions from molecular internal degrees of freedom. Here, we develop a simple working formula that yields good estimates for a general molecular (or atomic) solid and illustrate the approach by computing equilibrium vapor pressure of three different molecular solids and an atomic solid, as well as the desorption rate of a foreign (inhibitor) molecule from the surface of a molecular solid.
1. Introduction Kinetics of atoms and molecules at surfaces involving adsorption, desorption, and diffusion govern many important processes like sublimation, evaporation, condensation, crystal and thin-film growth, and action of chemical/molecular sensors and is relevant to diverse disciplines of science and technology. Equilibrium between the condensed phase (solid, liquid) and the gas phase (vapor) is established when the rate of sublimation/ evaporation from the condensed phase is equal to the rate of condensation from the vapor phase. Computing either of these quantities accurately is challenging. While the rate of influx of vapor-phase atoms/molecules onto the condensed surface is given simply by the Knudsen’s formula1,2 p/x2πmkBT, where p is the vapor pressure and m is the mass of the atomic/molecular species, the probability of incorporation into the surface, known as the sticking coefficient, as well as the net rate of desorption, may significantly depend on atomic details of the surface, the presence of a boundary layer, the energy and impingement direction of the incoming particle, and so forth. Focusing on solid-vapor equilibrium of atomic crystals and assuming that all desorption processes effectively start from “kink” sites1 of exposed crystal faces, one could estimate a net site desorption rate using an Arrhenius type expression:3
φ˘ vsimple ) νeff exp(-∆E/kBT)
(1)
where νeff is an effective vibrational frequency and ∆E is the desorption energy from the kink site, which is given by the difference in molecular potential energy between the vapor phase (g) and the kink site (s); that is, ∆E ) g - s. The physical reasoning behind eq 1 is that each atom vibrates with an effective Einstein’s frequency νeff and during each oscillation has a probability exp(-∆E/kBT) to desorb from the surface.4 Equating the incoming flux (given by Knudsen’s formula1,2) with the outgoing flux estimated by multiplying eq 1 with θkink, the surface density of kink sites, one obtains the following formula for the equilibrium vapor pressure p: * Corresponding author. E-mail:
[email protected].
psimple )
θkinkνeff 2πmkBT exp(-∆E/kBT) Rκ x
(2)
where κ is an average sticking coefficient and R is the area per surface site. Unfortunately, for molecular crystals, eqs 1 and 2 are not applicable because they do not include contributions from the molecular internal degrees of freedom. Thus, for concreteness, let us consider two energetic materials, that is, pentaerythritol tetranitrate (PETN) and the β-polymorph of 1,3,5,7-tetranitro1,3,5,7-tetrazocane (β-HMX), which consist of 29 and 28 atoms with masses of 316.2 and 296.2 au, respectively (Figure 1). PETN in its common form (PETN-I) crystallizes in a bodycentered tetragonal structure5 with a heat of sublimation ∆H of ∼35.1 kcal/mol,6 while β-HMX crystallizes in a monoclinic structure7 with a heat of sublimation of ∼44.2 kcal/mol.8 By assuming a νeff in the range 1012-1013 s-1 and a maximum kink density of θkink ) 1,9 eq 1 yields sublimation rates lower by more than 10 orders of magnitude as compared with the experimentally measured value for either PETN6 or HMX, while as will be shown below, eq 2 (assuming κ and θkink of the order of unity10) underestimates the equilibrium pressure of PETN and HMX by 11 to 12 orders of magnitude as compared to experiments. 2. Equilibrium Vapor Pressure: An Approximate Working Formula There are standard techniques in statistical mechanics, for example, reaction rate theory,3,11 or its more modern variants12,13 which allow one, in principle, to compute the site desorption rate for molecular crystals. However, an accurate estimate for a real system would involve a very accurate description of the vibron and phonon modes and their coupling, which would mandate either a specially developed interatomic potential (or force field) that may not exist for a molecular system of interest (as is shown later in this work) or a high-level quantum computation, which may not be feasible given the size of these molecular crystals. With Monte Carlo simulations of surface kinetics in mind,14 we deemed it useful to develop a formula
10.1021/jp076038g CCC: $37.00 © 2007 American Chemical Society Published on Web 11/30/2007
Role of Internal Degrees of Freedom
J. Phys. Chem. B, Vol. 111, No. 51, 2007 14291
6 ln{Zvib(νeff)} ) [
∫ dν gph(ν) ln{Zvib(ν)} 3NM-6
∑ j)1
ln{Zvib(νj)}] (6)
In experimental measurements of vapor pressure, one typically fits the slope of ln(p) versus 1/kBT to the heat of sublimation ∆H (assuming that ∆H does not vary significantly within considered temperature ranges). Thus, we rewrite eq 5 as:
p ) p0 exp(-∆H/kBT) where p0 ) Zrot
kBT Λ {Zvib(νeff)}6 3
exp[{∆H - ∆E}/kBT] (7)
Figure 1. Molecular models of (a) PETN (C5H8N4O12, NM ) 29), (b) β-HMX (C4H8N8O8, NM ) 28), (c) water (H2O, NM ) 3), and (d) argon (Ar, NM ) 1). Color scheme: C (gray), H (white), N (blue), O (red).
that would integrate the vibrational and coupled degrees of freedom into a single effective frequency νeff. By studying how this frequency varies with molecule size for a given training set, such a formula could be useful in predicting equilibrium vapor pressure and sublimation rate for a general molecular crystal. Since the actual sublimation rate depends on the details of surface conditions, let us focus for the moment on the equilibrium vapor pressure p. To derive a general formula for p, we use equilibrium statistical mechanics; that is, we equate the chemical potential in the solid (µs) with that in the vapor phase (µg). Using the harmonic approximation, we obtain:15,16
µs ) s - kBT
∫ dν gph(ν) ln{Zvib(ν)}
(3)
and
The quantity in the exponent of eq 7 can be evaluated using the Clausius-Clapeyron equation ∂lnp/∂(1/kBT) ) -∆H, which yields:
∆H - ∆E ) 4kBT - {
∫
3NM-6
gph(ν) dν.fvib(ν) -
4kBT - 6fvib(νeff) 1 -
∑ j)1
ln{Zvib(νj)} (4)
In eqs 3 and 4, s and g are the reference potential energies of the two phases, NM is the number of atoms in a molecule, gph(ν) is the phonon density of states (DOS) for the solid with a normalization of ∫ dν gph(ν) ) 3NM, and Λ ) h/x2πmkBT the thermal (or de Broglie) wavelength. The vibrational and rotational partition functions are given by15 Zvib(ν) ) exp(hν/2kBT)/{1 - exp(-hν/kBT)} and Zrot ) π1/2(8π2kBT)3/2(det[I])1/2/σh3, where det[I] is the determinant of the moment of inertia tensor (taken about the molecule’s center of mass) and σ is a symmetry factor given by the number of proper rotational symmetry operations for an isolated molecule.17 At thermodynamic equilibrium, the chemical potentials µs and µg must be equal, leading to the following expression for p:
p ) Zrot
kBT Λ {Zvib(νeff)}6 3
)
∂ ln νeff ∂ ln T
(8)
p ) p0 exp(-∆H/kBT) kBT Λ {Zvib(νeff)}6 3
exp[{4kBT 6fvib(νeff)}/kBT] (9)
3NM-6
kB T
fvib(νj)} )
where, fvib(ν) ) hν/2 + hν/(exp(hν/kBT) - 1) is the vibrational free energy associated with a mode of frequency ν, including the zero-point contribution. From numerical simulations on realistic molecular systems using several different force fields (see below), we find that the variation of νeff with T, that is, ∂ ln νeff/∂ ln T, is small, with an absolute value
