CRYSTAL GROWTH & DESIGN
Molecular Crystal Global Phase Diagrams†
2004 VOL. 4, NO. 5 1009-1012
J. Brandon Keith, Jonathan A. Mettes, and Richard B. McClurg* Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota 55455 Received January 20, 2004;
Revised Manuscript Received June 9, 2004
ABSTRACT: Global phase diagrams are proposed as a tool for reverse engineering of molecular crystal structures. Reverse engineering is the delineation of the family of intermolecular potentials consistent with an observed crystal structure. The method is illustrated using a two-dimensional slice through the global phase diagram summarizing the crystalline phase behavior of tetrahedral molecules. The diagram is used to discuss the phase behavior of heavy methane. We envision using these diagrams as tools for data mining from crystallographic databases, for feedback in crystal design, and for predicting possible polymorphs. Introduction Molecular crystals are formed to purify a chemical species, to stabilize a product, or to produce a crystal with desired properties such as nonlinear optical response. Given the wide range of applications, there has been considerable effort devoted to crystal structure prediction from molecular structure, with only modest success. In two recent blind tests, computer codes succeeded in predicting the observed crystal structure in only about 17% of the cases, even when allowed to make three guesses each.1,2 The challenge in crystal prediction is to enumerate possible structures and then to estimate their free energies with sufficient accuracy to rank order them. The large number of nearly isoenergetic structures observed for many molecules makes this a very difficult problem. Gavezzotti has even suggested that the sensitivity of crystal structure to intermolecular potential is too strong to permit reliable predictions.3 Given this extreme parameter sensitivity, we have chosen to tackle the inverse problem. Rather than attempting to predict crystal structure from intermolecular potentials, we seek to delineate the family of intermolecular potentials consistent with an observed crystal structure. The same extreme parameter sensitivity that hinders structure prediction allows us to place tight constraints on the set of consistent potentials. The product of our research is a phase diagram that summarizes the crystal phase behavior of entire classes of molecules as a function of their intermolecular potential. A phase diagram in which potential parameters are treated as independent variables is called a global phase diagram. Such diagrams are well-established tools for classifying and rationalizing liquidvapor phase behavior.4 By representing the phase behavior of an entire class of molecules with the same point group symmetry, we highlight the trends and relationships among molecular crystals. The outline of the balance of our paper is as follows. In the methods section, we introduce the foundations †
Presented in poster form at the 6th International Conference on the Crystal Growth of Organic Molecules, held in Glasgow, Scotland, August 2003. * To whom correspondence should be addressed. Tel: (612) 6249056. Fax: (612) 626-7246. E-mail:
[email protected].
of our methods. A detailed description is to be published elsewhere.5 To illustrate our model, we consider singlecomponent molecular crystals constructed from molecules having tetrahedral point group symmetry (e.g., methane, carbon tetrachloride, adamantane, and white phosphorus). A two-dimensional slice through the global phase diagram is given in the results section. There is no fundamental barrier to application of the method to less symmetric point groups except that the global phase diagrams become high dimensional and correspondingly more difficult to visualize. In the discussion section, we elaborate on the strengths and weaknesses of our method and use the global phase diagram to discuss the phase behavior of methane. Finally, we discuss some of the uses of the global phase diagrams for data mining, for crystal design, and for rationalizing and/or predicting polymorphs. Methods The Helmholtz free energy for a crystal of N identical molecules can be written as
A ) -kT(ln zN int + ln Ztr + ln Zor)
(1)
where zint is the partition function for the internal vibrational modes of a typical molecule. Ztr and Zor are partition functions for the translational and orientational modes of the crystal, respectively. With the exception of certain torsional modes, the internal modes are largely decoupled from the translations and rotations that determine crystal structure.6 Therefore, we consider rigid molecules. For nonionic molecules, there is generally a weak translation-rotation coupling with the orientational modes driving phase transitions.7 Therefore, our model focuses on the orientational modes of the crystal. Relaxation of translational modes in response to a phase transition is addressed in the discussion section below. The orientational partition function is calculated as follows
Zor ≡
1 hNd
∫∫ exp[-
]
Hor(ω, pω) dω dpω kT
(2)
where d takes a value of two for linear molecules and three for nonlinear molecules, Hor is the orientational Hamiltonian, ω is the vector of all ωi ) {Ri,βi,γi}, which is the set of Euler angles for the orientation of molecule i, and pω is the vector of conjugate momenta. Since the Hamiltonian is a sum of kinetic and potential contributions,
10.1021/cg0499616 CCC: $27.50 © 2004 American Chemical Society Published on Web 08/17/2004
1010 Crystal Growth & Design, Vol. 4, No. 5, 2004 Hor ) K + V
Keith et al. (3)
we substitute for the kinetic energy for either linear or nonlinear molecules and carry out the integration over momenta.8 The result is
Zor )
(
)∫
zrot
2(2π)
N
d-1
exp(-V/kT) sin(β)dω
(4)
where zrot is the partition function for a free rotor of dimension d.8 The term in front of the integral in eq 4 contributes a temperature-dependent, but crystal structure independent, additive term to the free energy (A). Therefore, it is the integral in eq 4 that determines crystal structure. To evaluate the orientational integral in eq 4, we require a functional form for the potential energy (V). We choose to represent the potential using an expansion in a complete set of orthogonal functions,
V)
1 2
∑ν
nσ,nµ nσ,nµ li,l,lj (rij)Fli,l,lj (ωi,Ωij,ωj)
(5)
ij
In eq 5, i and j index each of the molecules in the crystal, ν(rij) is an expansion coefficient as a function of molecular center separation (rij), F(ωi,Ωij,ωj) is a function of the orientation of each molecule (ωi and ωj) and of the orientation of the vector from molecule i to j (Ωij), and summation over li, l, l j, nσ, and nµ is implied. The basis functions (F) in eq 5 can be constructed to be consistent with molecular point group symmetries, greatly reducing the number of expansion coefficients (ν) needed to describe highly symmetric molecules. Expressing the functions F in terms of trigonometric functions of its eight independent variables (three Euler angles for each molecule and two spherical angles) leads to an unreasonably large expression. Rather than present them here, we direct interested readers to the references below. The expansion coefficients serve as independent parameters in our work. Thus, we consider the crystal structures of all molecules of a given point group symmetry simultaneously. Similar potentials are described in greater detail elsewhere.5,9-12 The potential in eq 5 may be factored as follows,
V kT
)
1 2
∑U
li lilj lj mτnσ(ωi)KmτnσmFnµ(Ωij)UmFnµ(ωj)
(6)
ij
In eq 6, summation over li, lj, mτ, nσ, mF, and nµ is implied. The functions U are called rotator functions,9 the thermal averages of which serve as crystal order parameters.13 There is a strong similarity between eq 6 and the Ising model of phase behavior, with the rotator functions serving as a lattice of spins coupled by the function K. However, rotator functions are continuous functions, and there are multiple spins per lattice site, unlike the simplest Ising model. Nevertheless, some of the same techniques used for Ising models are useful here.14,15 We adopt a mean field approximation to decouple the integrations over each molecule in eq 4. A Fourier transform of eq 6 aids in introducing crystalline translational invariance and in identifying the resulting crystal structure. Detailed discussions of these transformations and the solution of the resulting equations will be published elsewhere.5
Results To illustrate our model, we consider tetrahedral molecular point group symmetry (Td) and retain two terms in the potential energy expansion, eq 5. The corresponding expansion coefficients (denoted ν363 and ν343) are important for determining the phase behavior of methane.16 In the high temperature limit all of the order parameters vanish except the first which is always unity. This
Figure 1. Two-dimensional slice of the global phase diagram for tetrahedral molecules. The origin is the plastic crystalline reference state. Rays emanating from the origin correspond to different molecules as a function of decreasing temperature. Areas correspond to phases bounded by phase transition loci. Triple points are denoted with dots.
corresponds to a fully disordered plastic crystalline state.6 The molecular centers sit on an FCC lattice, but there is no orientational order. This structure is observed for molecules such as methane17 and C60,18 whereas most molecular crystals melt, sublime, or decompose prior to reaching this state. Even if the plastic crystalline state is not observed for a particular molecular crystal, it serves as a useful reference state for our purposes. For a particular choice of expansion coefficients (ν), the order parameters are calculated as a function of decreasing temperature. For large but finite temperatures, a subset of the order parameters are nonzero. These order parameters are consistent with the symmetry of the FCC lattice and thus do not indicate a phase transition. Instead, they indicate that the molecules are not completely free to rotate, but preferentially sample a set of energetically equivalent orientations. Further reducing the temperature leads to a phase transition marked by a larger set of nonzero order parameters. The labels of the order parameters and their relative magnitudes are a signature for the lower symmetry phase.19 Further reduction of the temperature can lead to additional phase transitions. The loci of phase transition points as a function of the expansion coefficients (ν) is shown in Figure 1. The origin is the high temperature fully disordered plastic crystalline reference state. Rays emanating from the origin correspond to different tetrahedral molecules as a function of decreasing temperature. Areas correspond to phases bounded by phase transition loci and labeled with their space group designation. Triple points are denoted with dots. In addition to the seven labeled phases, there is a thin triangular region near {-1.5,9} in which Pmmn is stable. The phase transitions along the v363 axis are consistent with prior work using only one expansion coefficient.9,20 Since two of the expansion coefficients were retained in drawing Figure 1, it is a two-dimensional slice through the complete global phase diagram. Retaining three coefficients leads to a threedimensional section of the global phase diagram and additional phases.5 To illustrate the use of the global phase diagram, consider the phase behavior of deuterated methane.17
Molecular Crystal Global Phase Diagrams
From 90 K down to 27.0 K, the crystal is the fully disordered plastic crystalline structure, space group Fm3 h m. This is consistent with our choice for the hightemperature reference phase in the center of Figure 1. From 27.0 K down to 22.1 K, the crystal has eight molecules in the unit cell. Two of these are disordered and six are ordered with space group Fm3 h mc. From Figure 1 we conclude that the value of ν363 is negative and ν343 is either positive or very small and negative. If this were not the case, one of the other structures would have been observed. Below 27.1 K, the phase of deuterated methane is believed to be the orthorhombic structure Cmca with 16 molecules per unit cell. This phase is not observed in normal methane due to quantum rotational effects. Note that Cmca does not appear in Figure 1. This is because the particular expansion coefficients (ν363 and ν343) do not admit a space group representation consistent with this space group symmetry. At least one additional coefficient is required to permit Cmca. The translations of the molecules from their FCC lattice sites are 2% or less of the corresponding lattice parameter and are consistent with the symmetries of the Cmca space group. This is consistent with our assumption that translations relax following molecular ordering, but do not drive phase transitions. Note that the low temperature phases are not necessarily cubic, even though the molecular centers remain near a fixed FCC lattice, due to symmetry breaking by molecular orientations. Discussion One of the important features of our model is that it is based on the minimization of the free energy (A) rather than simply the potential (V). Many of the currently available software packages for predicting crystal structure are based on the minimization of the potential energy under the assumption that the molecules adopt single low energy orientations at low temperature.1,2 While this is appropriate for perfect crystals in the limit of absolute zero, it is inadequate at finite temperatures and for crystals with residual disorder. Our model includes contributions from the kinetic energy and from the entire distribution of orientations, both of which are important contributions to the free energy at finite temperature. The translational contributions to the free energy have been neglected however. This leaves open the possibility that differences in the acoustic density of states of the crystals could alter the phase behavior. Our assumption is that the free energy contributions due to the acoustic modes of the phases are nearly equal along the phase transition loci. Although there may be substantial differences in the details of the acoustic mode density of states, we expect that their integrated contributions are similar. Construction of low-dimensional global phase diagrams is dependent on the representation of the interaction potential using very few parameters. Previous authors have used potential energy expansions similar to eq 5 and have concluded that the convergence of the summations is too slow for practical applications.10 The same authors acknowledge that this is surprising given the early success in rationalizing phase behavior using only one parameter.9 We believe that the resolution of
Crystal Growth & Design, Vol. 4, No. 5, 2004 1011
this conflict lies in the method for determining the expansion parameters (ν). The standard technique is to multiply both sides of eq 5 by one of the basis functions and integrate over all orientations. Due to the orthogonality of the basis functions, this provides an estimate of one of the expansion coefficients. Repeating the process for each of the basis functions gives values for each of the coefficients. This method is known to provide the best estimates of the coefficients in the least squares sense.21 The drawback to this approach is that it equally weights all possible orientations, whereas crystal structures contain low-energy interactions among molecules. We are currently refining a process for evaluating the expansion coefficients by equating the multidimensional Taylor’s series expansions of the potential and its approximation in eq 5. The result is a potential that is very accurate in the vicinity of the minima of the potential with only a handful of parameters. The potential is less accurate for high-energy interactions, but these are irrelevant for molecular crystals. Details of the method will be published elsewhere.16 Molecular crystal global phase diagrams are tools for reverse engineering of crystals. Reverse engineering is the delineation of the family of intermolecular potentials that is consistent with an observed crystal structure. Given a crystal structure, we are able to locate the structure within the phase diagram and express strong constraints on the intermolecular potential that are consistent with the observed structure. This has multiple uses. First, we have a tool for data mining to gain insight into the intermolecular forces that form crystals. Second, we are able to provide guidance in crystal design to produce desired crystal structures. Current practice in the synthesis of nonlinear optical material and organic electronic materials is to produce a molecule that has desirable properties and to apply heuristic rules in the hopes of influencing the crystal phase behavior. Unfortunately, many of the resulting structures are incompatible with the desired physical properties, and little or no information is gained from the observed structure if it is not the desired one. Using our global phase diagrams, the potential energy coefficients can be estimated for both the observed and desired crystal structures. The differences between the sets of coefficients provides clues for changes to the molecular structure that may result in the desired crystal structure. Finally, we can use the global phase diagrams to rationalize and/or predict crystal polymorphs. While we have focused on the thermodynamically stable structures in constructing Figure 1, our model also produces lists of metastable structures and their free energies. Low free energy metastable structures are candidate polymorphs. Conclusions We have presented preliminary results of a method to produce global phase diagrams of the solid-phase behavior for molecules of a given point group symmetry. Global phase diagrams are, by construction, idealized representations of the phase behavior of entire classes of molecules. They are useful for understanding the relationship between the phase behavior of similar molecules. We envision using these global phase diagrams as tools for data mining from crystallographic
1012 Crystal Growth & Design, Vol. 4, No. 5, 2004
databases, for feedback in crystal design, and for predicting possible polymorphs. Acknowledgment. We thank 3M for partial funding of this work through a 3M nontenured faculty grant. The University of Minnesota provided additional funding. We also gratefully acknowledge the use of Minnesota Supercomputing Institute facilities. References (1) Lommerse, J. P. M.; Motherwell, W. D. S.; Ammon, H. L.; Dunitz, J. D.; Gavezzotti, A.; Hofmann, D. W. M.; Leusen, F. J. J.; Mooij, W. T. M.; Price, S. L.; Schweizer, B.; Schmidt, M. U.; van Eijck, B. P.; Verwer, P.; Williams, D. E. Acta Crystallogr. B 2000, 56, 697-714. (2) Motherwell, W. D. S. et al. Acta Crystallogr. B 2002, 58, 647-661. (3) Gavezzotti, A. Acc. Chem. Res. 1994, 27, 309-314. (4) Van Konynenburg, P. H.; Scott, R. L. Philos. Trans. Royal Soc. (London) 1980, A298, 495-540. (5) Mettes, J. A.; Keith, J. B.; McClurg, R. B. Acta Crystallogr. A. 2004, in press. (6) Sherwood, J. N. The Plastically Crystalline State; WileyInterscience: New York, 1979. (7) Lynden-Bell, R. M.; Michel, K. H. Rev. Mod. Phys. 1994, 66, 721-762. (8) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000.
Keith et al. (9) James, H. M.; Keenan, T. A. J. Chem. Phys. 1959, 31, 1241. (10) Briels, W. J. J. Chem. Phys. 1980, 73, 1850-1861. (11) Stone, A. J.; Tough, R. J. A. Chem. Phys. Lett. 1984, 110, 123-129. (12) van der Avoird, A.; Wormer, P. E. S.; Moszynski, R. Chem. Rev. 1994, 94, 1931-1974. (13) Landau, L. D.; Lifshitz, E. M. Statistical Physics; Pergamon: New York, 1969. (14) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (15) Chaikin, P. M.; Lubensky, T. C. Principles of Condensed Matter Physics; Cambridge University Press: New York, 1995. (16) Missaghi, M. N.; Mettes, J. A.; McClurg, R. B. 2004, submitted. (17) Neumann, M. A.; Press: W.; Noldeke, C.; Asmussen, B.; Prager, M.; Ibberson, R. M. J. Chem. Phys. 2003, 119, 15861589. (18) Michel, K. H.; Copley, J. R. D. Z. Phys. B 1997, 103, 369376. (19) Stokes, H. T.; Hatch, D. M. Isotropy Subgroups of the 230 Crystallographic Space Groups; World Scientific: Teaneck, NJ, 1988. (20) Nagamiya, T. Prog. Theor. Phys. (Kyoto) 1951, 6, 702-713. (21) Davis, H. Fourier Series and Orthogonal Functions; Dover: New York, 1989.
CG0499616
