J. Phys. Chem. B 2008, 112, 13231–13237
13231
Investigating the Polymorphism in PR179: A Combined Crystal Structure Prediction and Metadynamics Study T. Zykova-Timan,* P. Raiteri,† and M. Parrinello‡ Computational Science, Department of Chemistry and Applied Biosciences, ETH Zu¨rich, USI Campus, Via Giuseppe Buffi 13, CH-6900 Lugano, Switzerland ReceiVed: April 6, 2008; ReVised Manuscript ReceiVed: July 23, 2008
The polymorphism of an industrial important pigment (PR179) was studied with a combination of standard crystal structure prediction and metadynamics. The former provided a starting set of candidate polymorphs whose structural and thermal stability were then probed by metadynamics. Moreover, metadynamics allowed for exploring the free energy surface to look for other possible polymorphs that were not included in the original set. Our calculations indicate that two structures have a high structural stability and are therefore good candidates to be found in experiments. The lower energy phase of the two indeed corresponds to the known polymorph, and we suggest that the other might be a metastable polymorph not yet experimentally discovered. I. Introduction Perylenes are well-known pigments and are widely utilized in the art industry and for the production of electrophotographic photoreceptors,1 photoconductors,2 and optical disks.3 These compounds are important due to their crystallochromy, i.e. the dependence of their color on crystal packing.4,5 However, the industrial synthesis of organic materials is often complicated by the existence of several metastable structures (polymorphs) of a single compound. In the photographic industry, a careful control on the crystallization procedure is therefore crucial for the production of sensitizing dyes,6 where a minor modification in the structure might drastically affect the resulting color. Crystallochromy and polymorphism are therefore closely related issues and a systematic search for polymorphs in organic pigments has an industrial relevance as well as fundamental importance. In the last decades, the increasing interest in understanding polymorphism stimulated computational crystallographers to develop advanced methods to tackle the crystal structure prediction (CSP) problem. Strategies for CSP generally involve several stages that include packing together of molecular units to form possible 3D tessellations, followed by exploration of the potential energy surface, i.e. at zero temperature, based on techniques such as genetic algorithms7,8 or Monte Carlo simulation. Such methods often yield many possible candidates within a narrow energy range, and even with the inclusion of more sophisticated criteria for candidate selection, it is not always possible to identify unambiguously the true structure.9 The state-of-the-art in the field is tested every few years through a blind structure prediction exercise,10-12 and the difficulty of achieving success is illustrated by the results of this process. Because of the computational demands of CSP, it is necessary to employ approximate force fields during the polymorph search. Hence, it would be easy to blame this factor for any failures. * Electronic address:
[email protected]. † Present address: Nanochemistry Research Institute, Dept. of Applied Chemistry, Curtin University of Technology, GPO Box U1987, 6845, Perth, Western Australia. ‡ Present address: Institut fu ¨ r Physik, Johannes-Gutenberg-Universita¨t Mainz, Staudinger Weg 7, 55099 Mainz, Germany.
However, another important contribution to the crystal stability is given by entropy, which is usually neglected in standard CSP or included in a later stage via the harmonic or quasi-harmonic approximations.13-15 These approximations have shown to be valid when the working conditions are far below the melting temperature,13-15 which is usually not the case for organic materials that have a low thermal resistance. Therefore, a new methodology is needed to change our standpoint from a quest on the potential energy surface to a search on the free energy surface (FES). A possible solution to this can be given by the metadynamics algorithm,16,17 which naturally includes entropic effects. The core of metadynamics is the recognition of the important collective variables (CVs) and the addition of a time dependent potential acting on these variables which discourages the system from remaining in the same state, allowing a fast exploration of the FES. This algorithm has shown a great efficacy in simulating solid-solid phase transitions of inorganic17-20 and simple organic systems.21,22 In particular, metadynamics successfully reproduced the benzene polymorphism by starting from a simple parallel stacking of the molecules in a cubic box and eventually finding all the experimental polymorphs.18,21 In this article, we apply a further extension of this approach which combines the classical lattice energy minimization (LEM) algorithms used in CSP with metadynamics22 to the study of the polymorphism of the Pigment Red 179 (PR179). PR179 (N,N′-dimethylperylene-3,4:9,10-bis(dicarboximide), molecular formula C26H14N2O4) is a molecule characterized by a flat skeleton with symmetric CH3 tails (see Figure 1). The experimental crystalline structure of PR179 has a monoclinic unit cell (V ) 873.19 Å3) with experimental lattice parameters: a ) 3.874 Å, b ) 15.58 Å, c ) 14.597 Å, β ) 97.65°, space group P21/ C, and Z ) 2.23 Figure 2 shows the projection of the pigment structure on the (a, c) and (b, c) planes. The molecules are packed in stacks along the a-axis with a spacing of ∼3.5 Å and an angle of ∼13° between the stacks. II. Method A complete description of how CSP and metadynamics can be successfully combined is reported elsewhere,22 and here, we
10.1021/jp802977t CCC: $40.75 2008 American Chemical Society Published on Web 09/30/2008
13232 J. Phys. Chem. B, Vol. 112, No. 42, 2008
Zykova-Timan et al.
Figure 1. Ball and stick representation of Pigment Red 179. A label is given to all nonsymmetry equivalent atoms to define the atomic charges in Table 2.
give only a brief description of this approach to guide the reader through the simulation procedure. The basic idea is quite simple, rather than applying metadynamics on a random or ad hoc chosen arrangement of molecules, as it was done for benzene,21 a more reasonable starting point might be given by LEM structures. Therefore, an initial set of putative structures was constructed with a CSP run and then metadynamics simulations were performed to probe their stability at ambient conditions and to explore the FES for new stable polymorphs that were not included in the original set. A. Force Field. Before entering into a discussion about the generation of the candidate crystal structures, we briefly describe the force field which was used in all our calculations. Due to the planarity and symmetry of the molecule, we decided to treat PR179 as a rigid body in the simulations. The intermolecular interactions were modeled by the empirical Williams99 (W99) force field:24
E(rij) )
qiqj A - 6 + B exp(-Crij) rij r
(1)
ij
where rij is the distance between atoms i and j and q is the atomic charge; see Table 1 for a list of the parameters used. The molecule’s geometry was optimized at the self-consistent field level with Gaussian using HF/6-31G** basis set25 and the atomic charges were successively fitted on the ab initio molecular electrostatic potential with the foreshortening of the C-H bonds, consistent with the force field (Table 2). The foreshortening of the C-H bonds was also used to determine the position of the H masses. This causes a slight change in the momentum of inertia of the molecule; however, we believe this to be negligible due to the size of the molecule. B. Crystal Structure Prediction. An initial set of putative structures was generated by Crystal Predictor.26 We limited our search to the most frequently encountered symmetry groups in organic compounds27 (P1, P1j, P21, P21/c, P21212, Pna21, Pca, Pbca, Pbcm, C2/c, Cc, C2, Pc, Cm, P21/m, C2/m, P2/c, C2221), and the candidate crystal structures were systematically generated via the application of proper symmetry operations on a set of suitable molecular clusters. The variation of cell angles and lengths were restricted in the range of 30°-130° and 2-50 Å,
Figure 2. Experimental structure projected along b (a) and a (b). Hydrogen atoms are not shown.
TABLE 1: Intermolecular W99 Force Field for H, C, N, and O Atoms (kJ/(mol Å))a atom type
A (kJ/(mol Å6))
B (kJ/mol)
C (Å-1)
C (sp3) C (sp4) O N H
1701.73 978.36 1260.73 1398.15 278.37
270363 131571 284623 102369 12680
3.60 3.60 3.96 3.48 3.56
a The standard geometric and arithmetic combination rules are used for interactions between atoms of different kind. See ref 24 for a detailed description.
TABLE 2: Atomic Charges Fitted from Ab Initio Computations atom type
atomic charge
O N C1 C2 C3 C4 C5 C6 C7 C8 H1 H2 H3 H4
-0.563 -0.236 -0.304 0.697 0.218 -0.328 0.035 -0.268 0.068 -0.007 0.153 0.132 0.142 0.173
respectively; the lower boundary for cell density was chosen as 300 kg/m3. We used a real space cutoff of 20 Å and the Ewald convergence parameter of 10-16. We performed a LEM (at 0 K) of 50 000 structures with Crystal Predictor by preserving the crystal group symmetry, but only the structures within ∼10 kJ/mol from the global minimum were considered for the metadynamics study. After a cluster refinement, 18 different putative structures remained (see Table 3). The global minimum is separated by a 4 kJ/mol energy gap from all the other structures, Figure 3a, and therefore, it is the best candidate to match the experimental structure. In order to
Investigating the Polymorphism in PR179
J. Phys. Chem. B, Vol. 112, No. 42, 2008 13233
TABLE 3: List of the Initial Putative Structures within 10 kJ/mol from the Global Minimum As Generated by Crystal Predictor26 with Their Lattice Energies (E) and Structural Informationa num
E (kJ/mol)
space group
Z
a (Å)
b (Å)
c (Å)
R (deg)
β (deg)
γ (deg)
V (Å3)
I IIb III IV V VI VIIb VIII IX X XI XIIb XIIIb XIVb XVb XVIc XVII XVIII
0.00 4.22 4.38 5.51 6.81 7.45 7.85 7.85 8.50 8.62 9.06 9.22 9.23 9.38 9.61 9.61 9.84 10.02
P21/c P21 P21 P21 Pna21 P21/c C2/c P1j P21/c P21/c P21 P21 P21/c Pna21 Pna21 P1j P21/c P1j
2 2 2 2 4 4 8 2 4 4 2 2 4 4 4 2 4 2
4.02 13.74 22.49 16.19 29.50 15.67 16.88 4.28 6.60 16.95 16.52 3.92 6.94 31.67 29.42 8.62 12.58 16.83
15.76 16.28 15.77 14.66 3.81 16.39 16.23 24.97 16.25 15.57 13.87 16.16 16.79 4.16 3.92 11.95 17.59 8.14
14.11 4.04 3.96 3.92 16.03 7.36 16.90 9.55 16.90 15.16 3.95 14.48 19.45 13.72 15.83 11.81 8.32 19.54
90.0 90.0 90.0 90.0 90.0 90.0 90.0 72.5 90.0 90.0 90.0 90.0 90.0 90.0 90.0 95.7 90.0 64.2
97.2 98.2 140.4 105.7 90.0 73.9 51.6 98.3 90.1 153.3 87.7 86.9 53.4 90.0 90.0 53.6 97.7 81.6
90.0 90.0 90.0 90.0 90.0 90.0 90.0 110.6 90.0 90.0 90.0 90.0 90.0 90.0 90.0 77.0 90.0 139.5
443.2 446.9 446.7 447.4 450.5 453.8 453.2 454.5 453.2 450.4 452.3 457.6 454.2 451.5 456.4 455.5 456.2 455.8
a The energies and volumes (V) are reported per molecule. b Structures that were unstable at ambient conditions after the free energy minimization are shown in italics (see text for details). c Structures that were unstable at ambient conditions during the metadynamics are shown in italics (see text for details).
global minimum, it was worth performing a further search on the FES with metadynamics to understand whether other stable polymorphs might be observed in nature. C. Metadynamics Framework. Metadynamics is a method based on the choice of a suitable set of CVs and on the application of a time dependent bias which boosts the exploration of the FES and eventually may lead to a reconstruction the free energy profile projected on the chosen CVs.16,28 Metadynamics has been successfully applied to many different fields ranging from chemical reactions29 to biophysics30 and solid-solid phase transitions.17-22 Following the procedure of Martonˇa`k et al.,17 we decided to use the lattice vectors as CVs and arrange them in an upper triangular matrix h in order to remove the unphysical rotation of the cell. In this formulation of metadynamics, the CVs are evolved in a steplike fashion and the driving force acting on the CVs is given by the derivative of the Gibbs free energy G(h) plus a time dependent potential which is constructed as a sum of Gaussians placed along the trajectory in the CVs space. The thermodynamic contribution to the driving force can be written as17
∂G/ ∂ hij ) -V[h-1(p - PI)]ji
Figure 3. (a) Scatter plot of the lattice energies at 0 K and (b) enthalpies at ambient conditions with respect to the volume for the 18 putative structures. The asterisk symbol is used for the structures that were unstable after the free energy minimization (see text for details).
verify this hypothesis and the accuracy of the force field, we performed a LEM run starting from the experimental structure where the experimental molecules were replaced by our rigid units. Finally, a powder diffraction pattern comparison of the LEM experimental structure against the one of our global minimum showed a perfect correspondence. Even though the CSP correctly identified the experimental structure with the
(2)
where V ) det h is the volume of the system, p is the internal pressure tensor, and P is the external hydrostatic pressure at which the simulations are performed. The tensor p can be readily evaluated during constant temperature and volume (NVT) molecular dynamics (MD) simulations, and it is evident that, given its statistical nature, long averaging runs are necessary to obtain a good accuracy on ∂G/∂h. However, due to the presence of hard and soft modes in organic crystals, the free energy is highly anisotropic and the use of spherical Gaussians would cause an inhomogeneous filling of the free energy basin and a reduction of the efficiency of the method. We therefore perform a change of coordinates20 to rescale the FES according to the stiffness of the elastic modes of the crystals. To use a more compact notation, we arrange the CVs in a six dimensional vector h˜ ) (h11, h22, h33, h12, h13, h23)T and define a new set of scaled coordinates as
13234 J. Phys. Chem. B, Vol. 112, No. 42, 2008
si ) √λi
∑ Oji(h˜j - h˜j0)
Zykova-Timan et al.
(3)
j
where λi are the eigenvalues of the Hessian matrix, Aij ) ∂2G(h˜)/ ∂h˜i∂h˜j|h˜0, which is diagonalized by the orthogonal matrix O. In this way, we obtained a free energy well that is spherical close to the minimum, h˜0. The thermodynamic force is then given by
∂G ) ∂si
∑ ∂G ∂h˜ j
j
Oji
(4)
√λi
The equations of motion ruling the CV’s evolution during a metadynamics run can therefore be written as
st+1 ) st + δs G (s) ) G(s) + t
∑
∂G t/ ∂ s |∂G t/ ∂ s|
(
ht+1 ) ht + δh
(5)
|s - st′|2 W exp -
t′
