close

Вход

Забыли?

вход по аккаунту

?

Exploring Polymorphism The Case of Benzene.

код для вставкиСкачать
Angewandte
Chemie
been investigated for almost 40 years, but its structures have
been only partially elucidated.
A small and simple molecule, benzene, has provided a
challenging testing ground for many theorists in spite of its
limited significance in chemical practice. Only three of the
experimentally observed phases have been correctly reproduced by computer simulations.[8–13] Here, by using a novel
theoretical algorithm[14] we describe that at finite temperature
there are at least seven possible stable phases. Our results
help to answer a number of pending questions and identify six
of the experimental structures. For the seventh, with limited
experimental data, we are able to propose a complete
structure.
A tentative phase diagram showing all the benzene phases
identified in the literature is shown in Figure 1. At room
Theoretical Chemistry
Exploring Polymorphism: The Case of Benzene**
Paolo Raiteri,* Roman Martoňk, and
Michele Parrinello
The prediction of crystal structures from the molecular
structure is of great importance for the chemical and
pharmaceutical industries as well as for fundamental science,
and it represents a great challenge for theoretical chemists
and physicists.[1–7] An even harder task is the computational
exploration of all possible polymorphs of a given compound.
Crystalline benzene is a case in point: its phase diagram has
[*] Dr. P. Raiteri, Dr. R. Martoňk,+ Prof. M. Parrinello
Computational Science
Departement fr Chemie und Angewandte Biowissenschaften
ETHZ, USI Campus
Via Buffi 13, 6900 Lugano (Switzerland)
Fax: (+ 41) 919-138-817
E-mail: paolo.raiteri@phys.chem.ethz.ch
[+] Permanent address:
Department of Physics
Faculty of Electrical Engineering and Information Technology
Slovak University of Technology
Ilkovičova 3, 812 19 Bratislava (Slovakia)
[**] We thank A. Laio for many interesting discussions, J. D. Dunitz for
his critical remarks, and A. Gavezzotti for many helpful suggestions
for improving the manuscript.
Angew. Chem. Int. Ed. 2005, 44, 3769 –3773
Figure 1. Proposed phase diagram of solid benzene. For phases I, II,
III, III’, and IV we follow the nomenclature of ref. [15]. As a consequence phase IV of ref. [16] becomes here phase V. The liquid–solid
boundary and all the points reported in the picture are taken from
refs. [15, 16]. Hypothetical phase boundaries are represented with
dashed lines.
temperature X-ray and Raman scattering[15, 16] have determined the existence of five phases, I, II, III, III’, and IV, but a
good agreement between theory and experiment has been
obtained only for the first three,[8–13] while phases III’ and IV
are more controversial. It has in fact been argued that
benzene III’ may be only a distortion of benzene III,[15] while
benzene IV has been hypothesized as being a polymer-like
compound, where new chemical bonds have been formed due
to the very high pressure.[16] However, more recent farinfrared spectroscopic data support the existence of both
phases as true benzene crystals.[17] Less detailed evidence is
available for the remaining phases I’ and V. The former has
been identified only on the basis of a discontinuity in the
lattice constants as a function of temperature, and it has been
suggested that the transformation benzene I!benzene I’ is a
second order transition.[15] Finally, in the high-temperature
region only Raman data can be called on to support of the
existence of another phase, here called benzene V.[16]
The practical relevance of crystal structure prediction has
stimulated the development of many computational protocols,[1–7] among which the most common approach is to
generate large numbers of structures and select the one with
the lowest enthalpy by means of different minimization
DOI: 10.1002/anie.200462760
2005 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
3769
Communications
techniques. Entropy is usually neglected or evaluated in the
harmonic lattice approximation. Here we take a different
approach and apply a new general method, called metadynamics,[14] where the search for stable polymorphs is guided by
the Gibbs free energy, taking fully into account the role of
temperature and pressure. In its application to solid–solid
phase transitions, metadynamics[18] shares the basic idea of the
Parrinello–Rahman technique,[19, 20] where the simulation box
edges act as an order parameter to favor a structural
transition. The technique allows a fast and accurate exploration of the Gibbs free energy surface as a function of this
order parameter. Early applications of this method to
inorganic crystals have been very encouraging.[18, 21]
The problem of determining all the possible polymorphs
of a given compound, consisting of a search for stable or
metastable crystal structures at a given temperature T and
pressure P, can be formulated as a search for minima of the
Gibbs free energy G(h) in the space of a suitable order
parameter, h. In the simulation of crystals, periodic boundary
conditions are commonly used, and if no defects are present,
the simulation box defined by its edges a, b, and c has to be
commensurate with the crystalline unit cell. Therefore, the
simulation box vectors, arranged in a 3 3 matrix, h = (a,b,c),
can be used as an order parameter that discriminates between
different crystal structures in the same way as in the
Parrinello–Rahman technique.[19, 20] It is convenient to take
the box matrix h as upper triangular, which effectively freezes
the rotations of the box and reduces the dimensionality of the
order parameter to 6. This can be easily achieved by a rotation
of the simulation cell so that the a vector is aligned with the x
Cartesian axis and the b vector lies in the xy Cartesian plane.
Metadynamics explores the Gibbs free energy G(h) in the
subspace of the order parameter h in a stepwise fashion by
means of a steepest-descent-like dynamics [Eq. (1)], where dh
htþ1 ¼ ht þ dh
t
jt j
ð1Þ
is the maximum allowed variation of the box dimensions. The
driving force, ft = @Gt/@h, is a gradient of the Gibbs free
energy plus a time-dependent potential that discourages the
system from going into previously visited configurations. This
potential is obtained as a sum of Gaussians placed on all
previously visited points in the space of the order parameter h
[Eq. (2)].
Gt ðh,tÞ ¼ GðhÞ þ
X
t0 <t
0
ðhht Þ2
W exp
2
2 dh
ð2Þ
The parameters W and dh determine the resolution in
energy and order parameter and can be chosen following the
empirical guidelines provided in ref. [18]. The derivative of
the Gibbs free energy can be calculated using Equation (3),
@G
¼ V½h1 ðpPji
@hij
ð3Þ
where P is the external pressure and p is the averaged
pressure tensor. The latter can easily be evaluated from
3770
2005 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
molecular dynamics runs at constant h and temperature.
Entropy is thus automatically accounted for through its
volume–pressure dependence. As the dynamics proceeds, the
added Gaussians modify the free energy landscape by filling
the potential energy well relative to the initial crystal
structure so that eventually the system must escape to a
new basin of attraction, corresponding to a new crystal
structure.
Each metadynamics step consists of: a) the equilibration
of the system at the desired h and temperature, b) the
calculation of the average pressure tensor p, and c) the
evolution of h according to Equation (1) followed by the
rescaling of the atomic coordinates. Each of steps a and b
consists of a 1-ps molecular dynamics run at constant
temperature (300 K) and h. The external pressure, P, is
fixed at 2 GPa.
Many force fields have been developed for studying
organic crystals (see, for example, refs. [22–24]). Because it is
known that force fields optimized on static crystal structures
may perform poorly in dynamic calculations, here we adopted
the GROMOS96 force field.[25] The molecules are fully
flexible and hydrogen atoms are explicitly considered. The
partial atomic charges were obtained by the standard RESP[26]
method and are 0.11 e for the carbon atoms and 0.11 e for the
hydrogen atoms. The long-range electrostatic interactions are
computed with the Particle Mesh Ewald (PME) summation.
The simulation supercell is replicated in space by means of the
periodic boundary conditions. All the molecules in the
simulation box are free to move and to find a suitable
periodic arrangement without any restriction on the number
of molecules in the unit cell. For the time-dependent potential
we chose the Gaussian width, dh, between 2 and 3 and the
Gaussian height, W, between 500 and 2500 kJ mol1, while the
W/dh2 ratio fulfills the guidelines given in ref. [18].
This method offers a number of competitive advantages
over standard approaches to crystal structure prediction. At
finite temperature the Gibbs free energy is much smoother
than the potential energy surface at T = 0. This eliminates
most of the shallow minima, leaving only the relevant ones.
Another consequence of explicitly considering a finite
temperature is that we can discover phases that are stabilized
by entropy. In our application to benzene, only a limited
number of thermodynamically stable phases are found, some
of which (benzene I and V) are unstable at T = 0 and P = 0. A
further advantage of this method is that it can be implemented on modest computer equipment like standard PCs. In
our case, a whole metadynamics run took a couple of hours on
a 2-GHz Pentium IV workstation. Clearly, however, the cost
of the method depends mainly on the intermolecular potential used for the MD simulations.
We tackled the problem of determining all the polymorphs of benzene assuming only the knowledge of its
molecular structure. As an initial step in our procedure, we
generated a periodic lattice with two molecules in the unit
cell, where lattice vectors and molecular orientations were
taken at random. Then we performed a pure steepest descent
relaxation of the box and of the atomic coordinates using a
simulation cell that contained no more than eight molecules.
Once the system relaxed to an energy basin, the temperature
www.angewandte.org
Angew. Chem. Int. Ed. 2005, 44, 3769 –3773
Angewandte
Chemie
was raised to room temperature and a short conventional
constant-pressure molecular dynamics simulation was performed. If the system so generated is unstable, the procedure
was started again.
We used tens of randomly generated crystal structures as
starting points for the metadynamics. Most of them have no
symmetry and are highly metastable, and a few metadynamics
steps were often sufficient to escape from the shallow energy
basin and obtain a more stable structure, which most of the
time is similar to benzene I or III. This structure was then
used as a building block for larger simulation cells used in
subsequent metadynamics runs. In a series of trial runs it was
observed that if we chose a MD simulation cell that was too
small, no structural transition was observed. On the other
hand, in large MD cells stacking faults were created at low
energy cost, thus preempting the crystal-to-crystal phase
transition. With a cell of 192 independent molecules it was
possible to observe smooth crystal-to-crystal transformations.
The structure factor Fhkl provides a means to monitor the
occurrence of a phase transition [Eq. (4)].[18] When a proF hkl ¼
X
! !
ei G hkl r j
ð4Þ
j
nounced change in a selected reflection of Fhkl was observed,
the dynamics was stopped and the new structure was
optimized with respect to the box edges and the internal
coordinates, both at T = 0 and at the experimental condi-
Table 1: Comparison of experimental[15] and theoretical data of the seven benzene polymorphs.[a]
Theory
Lattice paramenters
Expt
Space group
Theory
Expt
I
a
b
c
b
Vm
7.42
9.74
7.27
90.0
131.4
7.49
9.71
7.07
90.0
128.5
Pbca
Z=4
Pbca
Z=4
I’
a
b
c
b
Vm
7.48
9.47
7.14
90.0
126.4
–
–
–
–
–
Cmca
Z=4
–
–
–
–
–
II
a
b
c
b
Vm
5.72
5.72
14.6
90.0
119.4
5.54
5.54
15.3
90.0
117.4
P43212
Z=4
P43212
Z=4
III
a
b
c
b
Vm
7.51
5.50
5.50
110.0
106.7
7.44
5.20
5.31
109.4
96.88
P21/c
Z=2
P21/c
Z=2
III’
a
b
c
b
Vm
9.01
5.90
7.90
109.0
99.3
5.15
4.96
7.23
110.9
86.3
C2/c
Z=4
P21/c
Z=2
IV
a
b
c
b
Vm
9.27
5.87
6.40
90.0
87.1
9.13
4.96
6.46
101.8
71.6
Pbam
Z=4
Pbam
Z=4
V
a
b
c
b
Vm
5.60
4.04
9.52
95.0
107.3
–
–
–
–
–
P21
Z=2
–
–
–
–
–
[a] The cell lengths (a, b and c) are in , the angle b is in degrees, and the molecular volume (Vm) is in 3. All the crystalline phases have a = g = 90. The
theoretical data are computed at room temperature and experimental pressure,[15] except for benzene I’ and V, which are computed at T = 0 K and P =
0 GPa and T = 600 K and P = 5 GPa, respectively. The experimental data relative to benzene II are from ref. [11], and the experimental crystallographic
assignments of benzene III’ and IV are rather tentative.[15]
Angew. Chem. Int. Ed. 2005, 44, 3769 –3773
www.angewandte.org
2005 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
3771
Communications
tions.[15] This was the key step in the identification of the new
stable structures. The space group and the number of
molecules in the unit cell for the equilibrium structure were
determined by the program PLATON.[27] Then we simulated
the powder diffraction spectra S(2q) for all the equilibrated
structures. Comparison of the structural information and of
S(2q) with the available experimental data allowed the
matching of theoretical and experimental phases. In particular, the difference in the interplanar distances is quite low
( 6 %), and there is a good correspondence in the relative
peak heights of the S(2q). We label the structures found by the
metadynamics with the names of the experimental phases to
which they correspond. A most noticeable result is that
Table 1 exhausts the list of all the minima that are found by
metadynamics; it is not a selection of those structures that
better fit the experimental data.
Benzene I is found to be stable at room temperature, but
it converts into a different structure, benzene I’, upon
relaxation at 0 K. This supports the hypothesis of the
existence of a such a phase suggested in ref. [15], and
benzene I is an example of a phase stabilized at room
temperature by an entropic contribution. The two other
well-characterized phases, benzene II and III,[10, 11] are also
found correctly by the metadynamics. It is interesting to note
that benzene II, which is stable between 2 and 5 GPa,[15] has
the same packing as its inorganic analogue, borazine, at
ambient pressure.[28]
Moving to the high-pressure side of the phase diagram,
the metadynamics is able to answer some open questions.
Benzene III’ indeed exists, and it is not a modification of
benzene III as proposed in ref. [15] but a new monoclinic
phase in the space group C2/c, which may better match the
experimental data when one cell parameter and the number
of molecules in the unit cell are halved. At even higher
pressure, benzene IV comes into play, and we are able to
determine that it belongs to the Pbam space group. Finally, we
assign the seventh structure discovered, whose space group is
P21, to benzene V. Owing to the lack of structural information
the last assignment is tentative and should be confirmed by
experiment; moreover, this structure appears to be unstable
at pressures lower than 5 GPa. The quantitative agreement
between theory and experiment decreases with increasing
pressure. This is not surprising since the potential has been
modeled on low-pressure data.
In view of the very delicate balance between the different
energy contributions we find our results remarkably good.
The accuracy of our predictions, the absence of spurious
minima, and a comparison between the lattice energies
computed with the GROMOS96 and with the W99 force
fields (Table 2) give an a posteriori validation of the capability
of the force field used. It is worth noting that the
GROMOS96[25] force field was not optimized for studying
crystal structures. We were able to fill the experimental gaps
in the understanding of the phase diagram of benzene by
determining all the stable crystal structures. This success
demonstrates that our methodology is a new and powerful
tool for the fast and accurate determination of the relevant
polymorphs of an organic molecular crystal, and that it may
help solve the problem of crystal structure prediction. The
3772
2005 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Table 2: Lattice energies in kJ mol1 for all the benzene phases that are
stable at T = 0 and P = 0 computed with both the GROMOS96[25] and the
W99[22] force fields.[a]
II
III
III’
IV
GROMOS96
W99
0.00
0.25
1.25
2.07
0.00
0.43
0.84
0.37
[a] For both force fields we considered the molecules as rigid units and
we minimized the energy with respect to the lattice vectors and
molecular orientations. We observed that upon cooling benzene I
converts into I’ for the GROMOS96 force field, while the opposite
occurs for the W99 force field. The reference structure is benzene II. The
atomic charges are those used for the molecular dynamics simulations.
most promising aspect of our work is that it provides a
workable new paradigm in which the Gibbs free energy and
not the potential energy surface plays the determining role.
Since in most practical applications the relevant structures are
those which are stable at ambient conditions, this appears to
us a major step forward. This is especially true for substances
like most organic molecules for which T = 300 K is very close
to the melting temperature. Once the thermodynamically
relevant phases are found, the calculations can be further
refined with the use of more accurate potentials and thus the
relative enthalpy of the different phases can be determined.[11, 21] By using methods of statistical mechanics it is
also possible to draw a full phase diagram.[29] These extensions
are, however, outside the scope of the present work.
Received: November 30, 2004
Published online: May 11, 2005
.
Keywords: crystal structure prediction · molecular dynamics ·
phase transitions · polymorphism
[1] J. P. M. Lommerse, W. D. S. Motherwell, H. L. Ammon, J. D.
Dunitz, A. Gavezzotti, D. W. M. Hofmann, F. J. J. Leusen,
W. T. M. Mooij, S. L. Price, B. Schweizer, M. U. Schmidt, B. P.
van Eijck, P. Verwer, D. E. Williams, Acta Crystallogr. Sect. B
2000, 56, 697 – 713.
[2] W. D. S. Motherwell, H. L. Ammon, J. D. Dunitz, A. Dzyabchenko, P. Erk, A. Gavezzotti, D. W. M. Hofmann, F. J. J.
Leusen, W. T. M. Mooij, S. L. Price, H. Scheraga, B. Schweizer,
M. U. Schmidt, B. P. van Eijck, P. Verwer, D. E. Williams, Acta
Crystallogr. Sect. B 2002, 58 647 – 661.
[3] S. L. Price, Adv. Drug Delivery Rev. 2004, 56, 301 – 319.
[4] S. Datta, D. J. W. Grant, Nat. Rev. Drug Discovery 2004, 3, 42 –
57.
[5] R. J. Wawak, J. Pillardy, A. Liwo, K. D. Gibson, H. A. Scheraga,
J. Phys. Chem. A 1998, 102, 2904 – 2918.
[6] Y. Yonetani, K. Yokoi, Mol. Phys. 2001, 99, 1743 – 1750.
[7] J. Pillardy, Y. A. Aranautova, C. Czaplewski, K. D. Gibson,
H. A. Scheraga, Proc. Natl. Acad. Sci. USA 2001, 98, 12 351 –
12 356.
[8] T. Shoda, K. Yamahara, K. Okazaki, D. E. Williams, THEOCHEM 1994, 313, 321 – 334.
[9] T. Shoda, K. Yamahara, K. Okazaki, D. E. Williams, THEOCHEM 1995, 333, 267 – 274.
[10] M. M. Thiry, C. Rrat, J. Chem. Phys. 1996, 104, 9079 – 9089.
www.angewandte.org
Angew. Chem. Int. Ed. 2005, 44, 3769 –3773
Angewandte
Chemie
[11] B. P. van Eijck, A. L. Spek, W. T. M. Mooij, J. Kroon, Acta
Crystallogr. Sect. B 1998, 54, 291 – 299.
[12] V. E. Bazterra, M. B. Ferraro, J. C. Facelli, J. Chem. Phys. 2002,
116, 5984 – 5991 and V. E. Bazterra, M. B. Ferraro, J. C. Facelli, J.
Chem. Phys. 2002, 116, 5992 – 5995.
[13] Y. Yonetani, K. Yokoi, Mol. Phys. 2002, 100, 3915 – 3919.
[14] A. Laio, M. Parrinello, Proc. Natl. Acad. Sci. USA 2002, 99,
12 562 – 12 567.
[15] M. M. Thiry, J. M. Lger, J. Chem. Phys. 1988, 89, 4255 – 4271.
[16] F. Cansell, D. Fabre, J.-P. Petitet, J. Chem. Phys. 1993, 99, 7300 –
7304.
[17] L. Ciabini, M. Santoro, R. Bini, V. Schettino, J. Chem. Phys.
2001, 115, 3742 – 3749.
[18] R. Martoňk, A. Laio, M. Parrinello, Phys. Rev. Lett. 2002, 90,
75 503.
[19] M. Parrinello, A. Rahman, Phys. Rev. Lett. 1980, 45, 1196 – 1199.
[20] M. Parrinello, A. Rahman, J. Appl. Phys. 1981, 52, 7182 – 7190.
[21] C. Ceriani, A. Laio, E. Fois, A. Gamba, R. Martoňk, M.
Parrinello, Phys. Rev. B 2004, 70, 113 403.
[22] D. E. Williams, J. Comput. Chem. 2000, 21, 1154 – 1166.
[23] Y. A. Aranautova, A. Jagielska, J. Pillardy, H. A. Scheraga, J.
Phys. Chem. B 2003, 107, 7143 – 7154.
[24] J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J.
Comput. Chem. 2004, 25, 1157 – 1174.
[25] D. van der Spoel, A. R. Buuren, D. P. Tieleman, H. J. C. Berendsen, J. Biomol. NMR 1996, 8, 229 – 238.
[26] C. Bayly, P. Cieplak, W. Cornell, P. A. Kollman, J. Phys. Chem.
1993, 97, 10 269 – 10 280.
[27] A. L. Spek, J. Appl. Crystallogr. 2003, 36, 7 – 13.
[28] G. Raabe, Z. Naturforsch. 2004, 59, 609 – 614.
[29] E. Sanz, C. Vega, J. L. F. Abascal, L. G. MacDowell, Phys. Rev.
Lett. 2004, 92, 255 701.
Angew. Chem. Int. Ed. 2005, 44, 3769 –3773
www.angewandte.org
2005 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
3773
Документ
Категория
Без категории
Просмотров
0
Размер файла
368 Кб
Теги
exploring, polymorphism, case, benzenes
1/--страниц
Пожаловаться на содержимое документа