Evaluation of All-Atom Force Fields for Anthracene Crystal Growth

Article pubs.acs.org/crystal

Evaluation of All-Atom Force Fields for Anthracene Crystal Growth Peter Grančič,† Rita Bylsma,‡ Hugo Meekes,‡ and Herma M. Cuppen*,† †

Theoretical Chemistry and ‡Solid State Chemistry, Institute for Molecules and Materials, Radboud University Nijmegen, Heyendaalseweg 135, 6525 AJ Nijmegen, The Netherlands S Supporting Information *

ABSTRACT: Here we present a comprehensive evaluation of a set of wellknown all-atom force fields with the scope to model dynamic phenomena in molecular crystals composed of polyaromatic hydrocarbons. The capability of the force fields to reproduce experimentally and computationally available data is thoroughly scrutinized against anthracene molecular crystals that serve as a model system. First, the properties of the solid crystalline phase are investigated by employing geometry optimization using molecular mechanics. Because any inaccuracies can be easily overlooked in the constrained solid crystalline phase, the interaction energy of a variety of dimer conformations is obtained by employing an extensive local minima search algorithm. The larger configurational freedom in the dimer conformations better reflects the incorporation of molecules at the surface during crystal growth. The results are compared to known ab initio calculations as very little experimental data concerning the anthracene dimer are available. Finally, for three force fields with different performance in other tests, a polymorph prediction is carried out. Overall, we show that some of the selected all-atom force fields (BMM2, BMM3, W99, and isoPAHAP) perform remarkably well, whereas others (Amber, Bordat, Dreiding, DRESP, MM2, and MM3) fail to reproduce known computational data for a variety of reasons.

■

costs.4−6 As an example, dispersion corrections for the computationally feasible DFT method became generally accessible only recently, although its capabilities in this direction still need to be explored (see refs 11 and 12 and references therein). It is important to note that although many of the simplified isotropic all-atom force fields (where the constituent atoms are considered spherically symmetric) often perform remarkably well for aromatic systems without heteroatoms, there is an ongoing discussion whether such potentials should be replaced by models that take the anisotropic electron distribution into account.9,13,14 The quality of a given all-atom force field can be evaluated by its ability to reproduce experimentally and computationally available data with reasonable accuracy. In the present work, the ability of selected isotropic force fields (Amber,15 Bordat,16 Dreiding,17 MM2/3,18−23 isoPAHAP,24,25 and Williams26−28) and their combinations to predict properties of the anthracene molecular crystal is discussed in comparison with the anthracene gas-phase dimer structures and energies because the incorporation of anthracene molecules at the crystal interface will not necessarily be restricted to the configuration observed in the crystal structure. First, the crystal properties are determined after a simple geometry optimization via molecular mechanics. This includes the crystal unit cell parameters and the lattice energy. Second, the potential energy surface of a

INTRODUCTION Aromatic−aromatic interactions are at the core of many chemical and biological processes, including folding of DNA and proteins, packing of aromatic crystals, and conformational preferences of chain molecules. Among the aromatic molecules, polyaromatic hydrocarbons (PAHs) occupy a prominent position because of their specific physicochemical properties and potential applications. PAHs in the form of high-quality single crystals and thin films attract interest because of their promising semiconducting properties.1 Larger PAHs (number of C > 60) in stacked conformations are important to the astrochemical community as possible components of interstellar dust.2,3 Special attention was paid recently toward modeling the growth of single crystals of molecules belonging to the group of PAHs.4−6 To study molecular aggregation and crystallization, a realistic toolkit describing the intermolecular forces that determine the association of molecules into clusters is required. The complex nature of PAH molecules can be characterized by two types of bonds: weak intermolecular (controlling the crystal packing) and strong intramolecular (conserving the almost rigid shape of the molecules). Furthermore, the pronounced anisotropy of the surface energies is an important quantity controlling the orientation and shape of PAH crystals during growth.7−9 Although there is currently a trend to directly apply calculations from first-principles in prediction of intermolecular interactions in a crystal,10 an all-atom force-field approach remains widely used over decades because it provides sufficient level of accuracy at relatively low computational © XXXX American Chemical Society

Received: September 9, 2014 Revised: February 13, 2015

A

DOI: 10.1021/cg5013507 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Article

Crystal Growth & Design

van der Waals interactions are replaced by the Buckingham pair potential functions taken from MM2 and MM3 force fields, respectively. Also, the electrostatic charges are taken from the MM2/3 force fields. The isoPAHAP force field24,25 derived from a transferable anisotropic potential for PAH molecules9 was shown to reproduce intermolecular interactions in a variety of dimer conformations of smaller PAH molecules accurately. Because of simplicity and low computational demands, the intramolecular terms are taken from the Bordat force field. The point charges for this model were taken from the Supporting Information of the original paper that described the anisotropic potential,9 which are based on ESP calculations. The Williams 99 (W99) force field, originally parametrized for hydrocarbon crystals, provides only the definition of the intermolecular van der Waals interactions. Therefore, the intramolecular part is taken from the Dreiding force field, and the intermolecular part was completed by the addition of ESP point charges mentioned above. Algorithms. All molecular mechanics calculations were carried out by using the software EON39,40 and a set of homemade scripts to generate the input files and visualize the output. EON allows a user-defined force field to be programmed and included in the package. In the actual implementation, all the intermolecular van der Waals and electrostatic interactions are switched off if the center-of-mass distance between two interacting molecules, rcm, belongs to the interval [14, 15] Å using the switching function f(x) = (2x−3)x2 + 1 with x = rcm − 14 and x ∈ [0, 1]. Because EON does not allow for a true NPT geometry optimization, we use a combination of two built-in EON methods: (i) steepest-descent method that includes box modification and (ii) conjugate-gradient method with a constant box. A 5 × 7 × 5 supercell containing 350 anthracene molecules was used. For each force field, the supercell geometry is optimized within a loop of both optimization methods starting from the experimental crystal structure.41 The loop terminates when the convergence criterion (maximum force of 10−3 kJ mol−1 Å−1) is fulfilled. For the box, periodic boundary conditions are applied. The lattice energy of the crystal is calculated from the energy of the optimized monomer, E1, and supercell, E350, as Elatt = (E350 − 350E1)/350. Furthermore, the internal order of the unit cell is expressed via the main molecular axis angle, φ, the molecular plane dihedral angle, ψ, and the intermolecular center-of-mass distance, δ (Figure 1). To obtain the complete set of dimer local minima, the following algorithm is employed. First, a single (reference) molecule geometry is optimized. Second, a copy of the molecule is made and positioned at a point at certain distance relative to the center of mass of the reference molecule using a set of three coordinates and rotated by using a set of three Euler angles. Third, the reference molecule is frozen and the geometry is optimized with the convergence criterion of maximum force applied to any atom below 10−3 kJ mol−1 Å−1 via conjugate-gradient method. Because of symmetry and defined intermolecular interaction cutoff, this procedure can be limited to 1/8 of a sphere with a maximum radius of 14 Å (on the basis of the center-of-mass distance between the molecules). The translational and rotational steps were 0.1 Å and 10°, respectively. Each unique dimer conformation is saved and later minimized again with decreasing convergence criteria in a stepwise manner from 10−1 to 10−3 kJ mol−1 Å−1 with both molecules flexible. During this procedure, the dimer con-

dimer is sampled, allowing for identification of possible stable dimer conformations and corresponding interaction energies. Finally, we run polymorph prediction calculations for representative force fields selected with respect to their performance in other tests to illustrate their ability to predict experimentally observed crystal packing from scratch. All results are compared to available experimental and computational data.

■

METHODOLOGY Force Fields. The force fields considered in the present work were selected on the basis of their popularity across different communities of computational chemists. For various reasons, a number of standard modifications have been made. To reproduce crystal packing realistically, in most cases electrostatic interactions were added.29 To decrease computational costs, point charges were chosen instead of dipoles and the intramolecular interactions were simplified when possible to harmonic functions by neglecting the contribution of higher order terms. Here we give a short description of each force field. For a complete list of parameter values and functional forms we refer to the Supporting Information. The widely used computationally cheap Amber15 force field has been used in combination with AM1BCC point charges, generated using Antechamber software.30,31 The van der Waals and the electrostatic 1−4 interactions are set to act with weights of 1/2 and 5/6, respectively. The Dreiding force field17 has been applied in two ways: (i) without electrostatic interactions (Dreiding), as used in the original paper, and (ii) a version accompanied by ESP point charges (DRESP), as recommended in the original paper. The ESP charges are obtained at MP2/aug-cc-pVTZ level of theory (using MOLPRO32 for geometry optimization and Molden33 for the ESP point charges). The MM218,19 and MM320−23 force fields are widely used in modeling of both aliphatic and aromatic hydrocarbons. Because the aim of the present article is to identify an appropriate force field to model processes involved in anthracene crystal growth, the MM2/3 force fields have been slightly modified. The intermolecular interactions are emphasized over the intramolecular ones. Therefore, the intramolecular terms representing stretch−bend in MM218 and the terms representing stretch−bend, torsion−stretch, torsion−bend, and bend− bend in MM321 force field have been omitted. Furthermore, as suggested by Pettersson and Liljefors,34 the original dipole moments of the C−H bonds were replaced by atom-centered point charges on the order of ±0.15 electron units. The bond lengths and torsional constants were obtained by the method of bond orders.35−38 Both the van der Waals and the electrostatic nonbonded 1−4 interactions are set to act with the weight of unity. For the desired application, the resulting force fields become computationally more feasible, and the accented intermolecular interactions that dominate molecular aggregation processes retain their realism. The force field by Bordat and Brown16 developed for modeling polyacene crystals (here referred to as the Bordat force field) was thoroughly scrutinized, giving rise to two modifications. Although the single molecule geometries (defined by the equilibrium bond lengths, bond angles, and dihedral angles) and the electrostatic point charges are very similar to MM2/3 force fields, a careful comparison of the van der Waals pair potential reveals significant differences, especially in the case of C−H interactions. Therefore, in these modified force fields denoted as BMM2 and BMM3, the B


Article


identical structures), (iii) energy minimization of the structures using the Mechanics minimizer of Cerius2 and (iv) cluster analysis of the optimized structures.43 All long-range interactions were treated using Ewald summation with a constant value for epsilon. The algorithm was repeated 12 times to obtain reasonable statistics. The resulting structures (number of low energy polymorphs stored was 50) were compared to the experimental structures using the software MERCURY.44

■

RESULTS AND DISCUSSION Lattice Energy. The bulk phase of the anthracene crystal at ambient temperature and pressure can be characterized by a herringbone packing (where two molecules in the ab plane form layers perpendicular to the crystallographic c direction) within the monoclinic space group P21/a. The molecular mechanics geometry optimization of the 5 × 7 × 5 supercell reveals that most of the force fields perform rather well in terms of the lattice energy prediction (Table 1). Although direct comparison of the lattice energy Elatt to experimental ° should be done with care because sublimation enthalpy ΔHsub vibrations play an important role when converting calculated lattice energies into experimental sublimation enthalpies,7 one would expect to obtain similar absolute values. The lattice energy obtained from the original Bordat force field significantly underestimates this value. Similarly, the MM2/3 force fields give values that are a little bit lower than what is expected but for a different reason, as will be explained toward the end of this section. Also, the W99 force field underestimates the lattice energy. On the other hand, the Dreiding force field accompanied by ESP point charges overestimates the lattice energy. The unit cell parameters a, b, c, α, β, and γ give the first indication to the source of the problems. Although the unit cell angles (α, β, and γ) are in all present force fields relatively wellpreserved, the unit cell lengths (represented by the parameters a, b, and c) differ more significantly compared to the experimental values. To illustrate the problem at this stage of the discussion, let us compare the performance of two force fields, for example BMM2 and Bordat. The values of parameters a, b, and c are well preserved in the case of the former one, whereas the latter exhibits significant deviations from the experimentally observed values. The intermolecular interactions in the case of the Bordat force field are too attractive in the b direction and too repulsive in the other two crystallographic directions (a and c). To accommodate this, the geometry of the unit cell changes in the optimization procedure.

Figure 1. Definition of internal order parameters: φ, main molecular axis angle; ψ, intermolecular plane dihedral angle; and δ, the mutual center-of-mass distance. The dashed line indicates the unit cell.

formations are classified by means of the value of the convergence criterion down to which they remain stable. The dimer interaction energy is calculated from the energy of the optimized monomer, E1, and dimer, E2, as Eint = E2 − 2E1. For the comparison with conformations previously considered by other authors, the same script is used. First, the geometry is generated with a given distance and rotation. Second, a single point energy calculation is carried out with varying distance in a given direction. Predictions of possible crystal structures from scratch, without any experimental input, were carried out by employing the Polymorph Predictor module within the software package Cerius2.42 Because it is technically not possible to include all previously mentioned force fields into the package since there is a limited set of functional forms that can be used, this study is restricted to three force fields (Bordat, DRESP, and W99) that have shown to exhibit interesting behavior in the dimer tests. User-defined force-field input was used to ensure compatibility with EON. Because this work focuses on the intermolecular interactions, the geometry-optimized molecules were kept rigid during the full prediction procedure. The polymorph prediction procedure consists of the following steps: (i) search for reasonable crystal structures using a Monte Carlo simulated annealing method (given the molecular symmetry, the space group was specified as P21 instead of P21/a with the glide plane dropped, allowing for one independent molecule in the unit cell), (ii) cluster analysis of the resulting structures (removal of

Table 1. Unit Cell Parameters and Lattice Energies Obtained from the Supercell Geometry Optimization force field

a (Å)

b (Å)

c (Å)

α (deg)

β (deg)

γ (deg)

Elatt (kJ mol−1)

φ (deg)

ψ (deg)

δ (Å)

Amber BMM2 BMM3 Bordat Dreiding DRESP MM2 MM3 isoPAHAP W99 Exp.41,45

8.26 8.36 8.57 8.71 8.00 8.22 8.04 8.20 8.34 8.52 8.56

6.05 5.99 5.97 5.57 6.08 5.83 5.98 6.01 5.98 5.92 6.04

11.07 11.21 11.34 11.57 10.91 10.99 11.25 11.32 11.17 11.06 11.18

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

126.7 126.7 127.2 127.0 125.9 126.8 124.9 125.4 126.7 123.5 124.7

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

−105.1 −102.9 −98.1 −77.2 −104.5 −117.0 −92.0 −88.2 −102.7 −90.8 101.9 ± 1.3 (ΔHsub ° )

14.0 12.6 10.4 6.1 16.2 13.5 12.9 11.5 14.1 28.7 14.4

135.8 132.3 130.2 114.4 140.2 133.7 133.1 131.9 134.9 130.0 128.9

5.1 5.1 5.2 5.2 5.0 5.0 5.0 5.1 5.1 5.2 5.2

C


Article


Figure 2. Stable dimer conformations. The displayed dimer geometries correspond to all observed stationary points on the different potential energy surfaces tested as summarized in Table 2.

Figure 3. Comparison of the dimer interaction energies between selected force fields and SAPT(DFT) calculation by Podeszwa and Szalewicz,51 with dip denoting the interplanar distance and dcm denoting center-of-mass distance, respectively. The corresponding dimer geometries are displayed in Figure 2a−c and e.

force field is questionable because the molecules are forced by their neighbors to keep the correct positions. Furthermore, the initial configuration is taken from the crystallographic data, and such initial geometry offers sufficient restriction of the possible rearrangement of the molecules even if the geometry optimization procedure ignores the corresponding space group. There is very little data (experimental or computational) available on the stability of anthracene dimers in the gas phase. Careful study of the history of the anthracene-dimer problem reveals experimental evidence of two dimer conformations simultaneously present in the gas phase.46,47 Two possible conformations, crossed (symmetry D 2d) and T-shaped (symmetry C2v), were initially suggested (Figure 2a,b).46 It has been shown that the perpendicular T-shaped conformation must be considerably less stable than any stacked conformation. Therefore, the list of possible dimer geometries was enriched by adding the slipped-parallel conformation (symmetry C2h, Figure 2c).48,49 Although the interaction energy calculations were carried out at the MP2 level of theory, the initial geometries were identified by performing molecular dynamics simulations

One can investigate the geometry of the optimized supercell in more detail by looking at the internal order, represented by the parameters φ, ψ, and δ. Although the intermolecular distance δ within a unit cell is preserved in each case, the mutual main axes angles φ and molecular plane dihedrals ψ indicate significant geometry changes in the case of the Bordat force field. Furthermore, in the case of the W99 force field, the surprisingly large value of the order parameter φ is associated with reduced interlayer distance in the crystallographic c direction. Finally, it is worth noting that although the geometry optimization procedure does not take the crystallographic space group into account, the monoclinic character of the anthracene crystal is preserved for each of the selected force fields. Dimer Conformations and Interaction Energies. The capability of a given force field to reproduce the dimer conformation correctly is considered to be one of the most important characteristics when evaluating its quality. The reason for this lies in the fact that the crystal structure can be preserved upon geometry optimization even if the quality of the D


E

−17.09b −28.69c −32.03d −23.13d −22.38b −31.57b −31.03d −24.36d −16.91b −15.59d −11.86c −30.92d 3.37 5.17 7.37 3.97 3.38 3.60 3.63 4.48 4.76 5.51 7.37 7.13 7.30

−17.70b −28.10c −31.32d −21.26d −20.41c −32.49b,e −31.97d −25.01c −17.30b −15.84d −11.96d −17.09d −31.74d 3.27 5.12 7.33 3.88 3.28 3.51 3.54 4.54 4.73 5.47 7.33 7.13 7.30 5.54 4.57 15 5/3/7 4.61 14 6/2/6

−32.26b,e −27.04b

BMM3

−31.63b −27.96b

BMM2

13 5/5/3

3.79 4.05 4.11 5.38 5.47 5.14 7.01 6.83 6.94 5.36

7.03

3.29 4.77

−3.50c 1.00d 1.34c −20.13c −20.76c −25.28d −20.52b −19.63b −14.88b −15.39c

−20.14d

−3.59b −28.45b,e

Bordat

dcm (Å)

−26.35c −17.03b −16.22d −13.22c

−22.30d −14.70b −13.74d −12.81d

10 6/1/3

13 4/6/3

5.44 7.26 7.07 7.13

5.46 7.27 7.07 7.13

3.39 3.50 3.51

7.27

3.37 5.04 3.70 7.25 3.83 3.79 3.39 3.56 3.57

3.38 5.07

−14.83b

−48.35c −51.89b,e −51.84b

DRESP −46.83b,e −29.19b −41.10c −17.71c −42.34b −42.27c −46.55d −40.66d −40.08c

−48.36b −24.48b

Dreiding

4.57 15 6/1/8

7.42 4.72 3.87 3.28 3.51 3.54 4.50 4.73 5.43 7.43 7.14 7.29

3.27 5.06

−18.73d

−12.66b −14.30d −5.67d −4.12d −5.15d −5.53b −15.86b −17.29d −17.43c −10.71b −8.17d −6.45d

−1.72b −23.53b,e

MM2

7.44 4.76 3.95 3.37 3.61 3.63 4.49 4.76 5.43 7.45 7.15 7.29 5.54 4.61 16 5/3/8

3.37 5.12

−12.34b −16.27d −10.57c −7.76d −11.28d −10.40c −17.65b −19.36d −18.12c −11.08b −8.51d −6.29d −9.07d −20.40d

−7.69b −23.69b,e

MM3

9 5/2/2

7.36

7.28

9 5/2/2

5.64 7.45

3.57

5.58 7.38

3.62

3.38 5.26 3.76 7.45 3.91

−10.01d

−11.48d

3.38 5.17 3.82 7.38 3.98

−20.41c −13.66b

−31.57d

−38.78b,e −22.53b −33.41b −14.23c −34.64b

W99

−22.75c −15.09b

−30.49d

−36.99be −25.44b −32.56b −15.87c −33.56b

isoPAHAP

Dimer conformations are classified by means of the value of the convergence criterion down to which they remain stable, i.e., maximum force below b10−3, c10−2, d10−1 kJ mol−1 Å −1. eConformation with the lowest energy for a given force field.

a

−1

Amber

Conformation Eint (kJ mol ) a −32.50b b −29.18b c d −18.52b e −27.25d f −30.23d g −32.27d h −25.46d i −24.73b j −32.27c k −31.92d l −25.69b m −17.90b n −16.22d o −12.69d p −17.57d r −32.86b,e Conformation a 3.39 b 5.10 c d 7.32 e 3.83 f 3.91 g 3.39 h 3.55 i 3.58 j 4.44 k 4.76 l 5.53 m 7.34 n 7.14 o 7.21 p 5.53 r 4.57 no. of DSPs 16 with a/b/c 7/1/8

force field

Table 2. List of Dimer Stationary Points (DSPs)a

Crystal Growth & Design Article


Article


the MM2 and MM3 force fields, the true minima for a given configuration are completely different (Table 2). The reason behind this discrepancy is the additional effect of the intramolecular terms. This can be revealed by comparing two force fields that share identical definition of the intermolecular interactions and similar equilibrium geometries of isolated molecules, e.g., BMM2 and MM2 force fields. The geometry of single molecules in the BMM2 force field is enforced much more strongly than in the MM2 force field. Because of the interactions with the other molecule in the dimer, the weaker MM2 force field allows displacements of atoms from their equilibrium positions even at T = 0 K, increasing the intramolecular potential energy. It should be emphasized that the significant intramolecular contributions to the dimer energy for MM2/3 force fields do not represent a true physical effect but rather an artifact arising from errors in the force-field parametrization. Although the predictive ability of all-atom force fields to identify possible potential energy minima for a given system is questionable, we carried out an advanced sampling of the potential energy surface with the goal to identify all stationary points that our set of force fields can predict and to compare the corresponding geometries and energies to the known computational data (because no experimental values are available). All observed relatively stable dimer conformations are listed in Figure 2 and Table 2 (together with the corresponding convergence criteria of maximum force applied to any atom). For the Amber force field, the most stable conformation is the tilted dimer with the geometry close to that observed in the crystal structure. Very similar character and interaction energy values are observed in the cases of the BMM2 and BMM3 force fields. As already pointed out, the original Bordat force field performance is conceptually somewhat upside down. By looking at the Bordat force field dimer interaction energies, it becomes obvious that the T-shaped conformations are an order of magnitude more stable than the stacked ones. Furthermore, two parallel conformations have repulsive interaction energies. The Dreiding-based force fields show the same trend as the Amber and BMM2/3 force fields with much stronger binding. The DRESP, isoPAHAP, and W99 force fields are the only three that allow the existence of a stable slipped-parallel conformation. The slipped-parallel conformation is not stable in any other force field. The interaction energies obtained in the case of isoPAHAP force field are very similar to the reference calculations by Podeszwa and Szalewicz even if a complete all-atom optimization procedure is employed. By looking at the interaction energy values (Table 2), the stable dimer conformations can be divided into two classes: the stacked geometries and the perpendicular ones. Except for the original Bordat, MM2, and MM3 force fields, the general trend is to favor stacked conformations. This is in good agreement with the SAPT(DFT) calculations by Podeszwa and Szalewicz.51 The stability of the MM2 and MM3 dimer conformations shows a similar trend to the Bordat force field but for a different reason. Because of the internal energy contributions, the dimer interaction energies of the MM2 and MM3 force fields are significantly lower, compared to those of the BMM2 and BMM3 force fields. The fact that the graphite conformation cannot be stable in the case of BMM2 and BMM3 force fields even though it represents the most stable conformation in the case of MM2 and MM3 force fields points toward the

and geometry optimization using the MM3 force field within the Tinker package. To complete the list of previously discovered dimer conformations at this stage, one should not omit another T-shaped conformation with inplane long axes of the two monomers perpendicular to each other and a partially displaced center-of-mass (Figure 2d).49 The same authors later carried out another study using the Hartree−Fock dispersion model, considering yet two other conformations, the graphite conformation (symmetry Ci, Figure 2e) and another stacked conformation (symmetry approximately C2, similar to our Vshaped dimer, Figure 2f).50 Some of these conformations were later adopted by Podeszwa and Szalewicz to perform SAPT(DFT) calculations without a priori knowledge of whether they really represent local/global minima or saddle points on the potential energy surface.51 To the best of our knowledge, no comprehensive optimization was ever carried out that would confirm the existence of these minima. Figure 3 shows the comparison of the interaction energies between the selected all-atom force fields and the SAPT(DFT) calculation by Podeswa and Szalewicz51 (see Figure 2a−c and e for the corresponding geometries). A figure containing a comparison against all force fields studied is provided in the Supporting Information. It becomes obvious that the Bordat force field is absolutely unable to describe the interaction energy of any of the four dimer conformations. Furthermore, Dreiding and DRESP force fields both overestimate the well depth for the crossed, graphite, and slipped-parallel conformations and underestimate the well depth of the T-shaped conformation. On the other hand, force fields that include the MM2/3 intermolecular terms (BMM2, BMM3, MM2, and MM3) together with Amber and W99 force fields perform rather well. The isoPAHAP force field best resembles the SAPT(DFT) calculations, which is not too surprising because this force field is based on SAPT(DFT) calculation on a large set of different PAH structures. For all four conformations, the positions of all four minima in comparison to the SAPT(DFT) calculations are well represented in the case of isoPAHAP, BMM2 and BMM3 force fields. This is followed by Amber, Dreiding and DRESP force fields being shifted slightly toward lower intermolecular distances and W99 being slightly off in the case of T-shaped dimer conformation. The minima shifts are more pronounced in the case of MM2 and MM3 force fields and completely off in the case of the Bordat force field. The isoPAHAP force field, reflecting the anisotropic character of the PAH molecules, is the only force field capable of correctly reproducing the interaction energies of both stacked and perpendicular conformations simultaneously. All other force fields are capable of fitting only one group of the dimer conformations. As shown by Totton and co-workers9 in the case of isoPAHAP and W99 force fields, this represents a usual problem when dealing with anisotropy of PAH interactions. We note that except for the energy difference between the selected conformations both force fields predict the same set of dimer minima with the same convergence criteria. Totton et al.25 attribute the better performance of isoPAHAP to a combination of the presence of the Tang− Toennies damping term and to fact that it has been fitted specifically to PAH molecules. W99 was simply fitted to a much broader group of molecules. It must be pointed out that the curves in Figure 3 represent single-point energy calculations; as it turns out, in the case of F


Article

Crystal Growth & Design importance of the intramolecular contributions. A very tiny difference in the single-molecule geometry can result in a completely different potential energy. To quantify this, one can look at the energy contributions arising from the internal terms in the case of the crossed dimer by MM2 and MM3 force fields. The potential energy difference between the monomer when optimized isolated or in the presence of another monomer molecule is 15.79 and 13.02 kJ mol−1 for MM2 and MM3, respectively. Although almost invisible when inspected by the naked eye, the slightly misshapen molecular geometries contribute to the overall energy significantly. When discussing the intramolecular contributions, it should be pointed out that the Amber force field also exhibits a tendency toward misshapen geometries (not shown here). In contrast to MM2 and MM3 force fields, this only affects the flexibility of the C−C−H angle, probably as a consequence of the general purpose of the Amber force field. Pointing our attention toward the mutual center-of-mass distances between the dimer molecules reveals general agreement, except for the original Bordat force field. Polymorph Prediction. To check whether our dimer analysis can be used as a quality assessment of a force field to use in the solid state, we carried out some polymorph predictions. For this purpose, we built the crystal from scratch using the Polymorph Prediction module in the package Cerius2.42 As already mentioned, it is not possible to include all the considered force fields into the package because of the missing Tang−Toennies description. Therefore, three representative force fields (DRESP, W99, and Bordat) that exhibited different behavior during the dimer test have been selected. Figure 4 plots the lattice energy per fragment versus the obtained density for the predicted structures of five independent prediction runs for each of the applied force fields; including more runs does not lead to visually different graphs. The black circles correspond to structures with a similarity larger than 0.95 as compared to the experimentally determined structure. This similarity is determined by the program MERCURY on the basis of the simulated powder X‑ray diagram. From Figure 4, it is clear that a large number of polymorphic forms are predicted by the Bordat force field, many of which have a very low density. The experimental density at room temperature is 1.28 g cm−3. Furthermore, the predicted structures are compared against the optimized experimental structure using the same Bordat force field because structure optimization resulted in differences from the experimental structure too great for MERCURY to recognize them to be the same. Using the optimized experimental structure, MERCURY selects many structures to be similar to the optimized structure, covering a relatively larger area of the energy versus density plot. In agreement with the experimental supercell-geometryoptimization procedure, the corresponding lattice energies of these matching structures significantly underestimate the experimental sublimation enthalpy, Elatt ≈ −56 kJ mol−1. Figure 4 further shows that the Polymorph Predictor finds another set of structures with higher density and more stable lattice energies using the Bordat force field. Because the Tshaped dimer conformation represents the global minimum on the potential energy surface for the Bordat force field, it follows that when a Bordat crystal is constructed from scratch, it resembles this perpendicular pattern (with the internal order parameters φ = 0.4°, ψ = 89.5°, and δ = 4.4 Å). Consequently,

Figure 4. Lattice energy per fragment versus the obtained density for each of the predicted structures using the Bordat, DRESP, and W99 force fields. The black circles correspond to structures with a similarity larger than 0.95 as compared to the experimentally determined structure or, in the case of the Bordat force field, the optimized experimental structure.

the unit cell becomes significantly misshapen (a = 6.28 Å, b = 6.29 Å, c = 12.03 Å, α = γ = 90°, and β = 100.5°). The DRESP force field found a number one ranking in all prediction runs with a structure very similar to the experimentally observed ones, and the obtained similar structures all correspond to a well-defined lattice energy and density combination that is well separated from the other structures in the graph. The tilted pattern is obtained (with the internal order parameters φ = 16.4°, ψ = 134.2°, and δ = 5.2 Å). However, the resulting unit cell is slightly contracted in the crystallographic c direction and misshapen (a = 8.42 Å, b = 6.07 Å, c = 9.39 Å, α = γ = 90°, and β = 103.9°). As shown during the experimental supercell-geometry-optimization procedure, the lattice energy is overestimated, Elatt = −115.1 kJ mol−1. G



■

Finally, polymorph prediction runs using the W99 force fields proved to be difficult without any preprocessing because of Buckingham potential used for the van der Waals interaction without any damping term. Several very-low-energy and veryhigh-density structures were obtained, where the atoms of different molecules collided because of the very attractive interaction at close interatomic distances. The H−H interaction is particularly problematic because its maximum already lies at 1.1 Å and it has a relatively low repulsive barrier of roughly 200 kJ mol−1. However, if we exclude these unphysical structures, then the performance of the W99 force field in the prediction of the anthracene crystal structure is very similar to that of the DRESP force field as can be seen from Figure 4, which does not include the unphysical structures. The structure with the mostnegative lattice energy corresponds to the experimental structure, is well separated from the other predicted structures, and corresponds to a well-defined lattice energy and density combination.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This project is supported by the European Research Council (ERC-2010-StG, Grant Agreement No. 259510-KISMOL) and The Netherlands Organization for Scientific Research, NWO (VIDI Research Programme 700.10.427). We are grateful to Ad van der Avoird and Leendertjan Karssemeijer for their contribution to this work.

■

REFERENCES

(1) Kim, K. S.; Tarakeshwar, P.; Lee, J. Y. Chem. Rev. 2000, 100, 4145−4185. (2) Rapacioli, M.; Calvo, F.; Spiegelman, F.; Joblin, C.; Wales, D. J. J. Phys. Chem. A 2005, 109, 2487−2497. (3) Rapacioli, M.; Calvo, F.; Joblin, C.; Parneix, P.; Toublanc, D.; Spiegelman, F. Astron. Astrophys. 2005, 460, 519−531. (4) Cuppen, H. M.; Graswinckel, W. S.; Meekes, H. Cryst. Growth Des. 2004, 4, 1351−1357. (5) Massaro, F. R.; Moret, M.; Bruno, M.; Rubbo, M.; Aquilano, D. Cryst. Growth Des. 2011, 11, 4639−4646. (6) Massaro, F. R.; Moret, M.; Bruno, M.; Aquilano, D. Cryst. Growth Des. 2012, 12, 982−989. (7) Reilly, A. M.; Tkatchenko, A. J. Chem. Phys. 2013, 139, 024705. (8) Nabok, D.; Puschnig, P.; Ambrosch-Draxl, C. Phys. Rev. B 2008, 77, 245316. (9) Totton, T. S.; Misquitta, A. J.; Kraft, M. J. Chem. Theory Comput. 2010, 6, 683−695. (10) Sancho-García, J. C.; Aragó, J.; Ortí, E.; Olivier, Y. J. Chem. Phys. 2013, 138, 204304. (11) Klimeš, J.; Bowler, D. R.; Michaelides, A. J. Phys.: Condens. Matter 2010, 22, 022201. (12) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92, 246401. (13) Price, S. L.; Stone, A. J. J. Chem. Phys. 1987, 86, 2859. (14) Ambrosch-Draxl, C.; Nabok, D.; Puschnig, P.; Meisenbichler, C. New J. Phys. 2009, 11, 125010. (15) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179−5197. (16) Bordat, P.; Brown, R. J. Chem. Phys. 2009, 130, 124501. (17) Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III, J. Phys. Chem. 1990, 94, 8897−8909. (18) Allinger, N. L. J. Am. Chem. Soc. 1977, 99, 8127−8134. (19) Sprague, J. T.; Tai, J. C.; Yuh, Y.; Allinger, N. L. J. Comput. Chem. 1987, 8, 581−603. (20) Allinger, N. L.; Lii, J.-H. J. Comput. Chem. 1987, 8, 1146−1153. (21) Allinger, N. L.; Yuh, Y. H.; Lii, J.-H. J. Am. Chem. Soc. 1989, 111, 8551−8566. (22) Lii, J.-H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 111, 8566− 8575. (23) Lii, J.-H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 111, 8576− 8582. (24) Totton, T. S.; Misquitta, A. J.; Kraft, M. Chem. Phys. Lett. 2011, 510, 154−160. (25) Totton, T. S.; Misquitta, A. J.; Kraft, M. Phys. Chem. Chem. Phys. 2012, 14, 4081−4094. (26) Williams, D. E. J. Mol. Struct. 1999, 485, 321−347. (27) Williams, D. E. J. Comput. Chem. 2001, 22, 1−20. (28) Williams, D. E. J. Comput. Chem. 2001, 22, 1154−1166. (29) Williams, D. E.; Starr, T. L. Comput. Chem. (Oxford, U.K.) 1977, 1, 173−177.

■

CONCLUSIONS A relatively simple set of tests allows us to scrutinize the quality of a given force field for potential use in modeling anthracene crystal growth. Problems that do not become obvious by looking at the crystal properties can be revealed by sampling the potential energy surface of small clusters such as dimers and/or by predicting the possible crystal polymorphs from scratch. To generalize our findings, the overall best performance is exhibited by the isoPAHAP force field, which is closely followed by the performance of the W99, BMM2, BMM3 and Dreiding-based force fields (Dreiding and DRESP). The Amber force field, originally developed for a large variety of organic molecules, is too general to be used specifically for PAH systems. As such, it suffers from an orientational freedom of the C−H bond that is too large. Although completely perfect when predicting crystal values, it fails to reproduce the dimer conformations in a very specific manner. One can only speculate precisely why only the C−H intermolecular pair potential constants in the Bordat force field are so significantly problematic. However, the intramolecular description of the molecules in the Bordat force field makes it among the most feasible intramolecular potential because of the use of harmonic functions with sufficient stiffness that ensure the preservation of the rigid shape of the molecules. Although originally considered a generalized force field for any PAH molecule, its use in this direction is not straightforward. The MM2/3 force fields provide a rather good description of the intermolecular potential and forces. According to our experience with the anthracene molecule, when accompanied by stiffer definitions of intramolecular forces its only drawback lies in its isotropic character. The same holds for the W99 force field, which by definition requires the user to select appropriate intramolecular and electrostatic potentials. As expected, the performance of the isoPAHAP force field is the best. When accompanied by an appropriate intramolecular potential, it is expected to serve well in a large variety of problems involving walks on the potential energy surface.

■

Article

ASSOCIATED CONTENT

S Supporting Information *

All force field parameters used in the calculations and full version of Figure 3. This material is available free of charge via the Internet at http://pubs.acs.org. H


Article

Crystal Growth & Design (30) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. J. Comput. Chem. 2000, 21, 132−146. (31) Jakalian, A.; Jack, D. B.; Bayly, C. I. J. Comput. Chem. 2002, 23, 1623−1641. (32) (a) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. WIREs Comput Mol Sci 2012, 2, 242−253. (b) Werner, H.J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M. MOLPRO, version 2012.1; University College Cardiff Consultants Limited: Cardiff, U.K., 2012; a package of ab initio programs. http://www.molpro.net. (33) Schaftenaar, G.; Noordik, J.H. Molden: a pre- and postprocessing program for molecular and electronic structures. J. Comput.Aided Mol. Design 2000, 14, 123−134 http://www.cmbi.ru.nl/molden/ molden.html. (34) Pettersson, I.; Liljefors, T. J. Comput. Chem. 1987, 8, 1139− 1145. (35) Allinger, N. L.; Sprague, J. T. J. Am. Chem. Soc. 1973, 95, 3893− 3907. (36) Kao, J.; Allinger, N. L. J. Am. Chem. Soc. 1977, 99, 975−986. (37) Liljefors, T.; Tai, J. C.; Li, S.; Allinger, N. L. J. Comput. Chem. 1987, 8, 1051−1056. (38) Allinger, N. L.; Li, F.; Yan, L.; Tai, J. C. J. Comput. Chem. 1990, 11, 868−895. (39) Pedersen, A.; Jónsson, H. Math. Comput. Simul. 2010, 80, 1487−1498. (40) Chill, S. T.; Welborn, M.; Terrell, R.; Zhang, L.; Berthet, J.-C.; Pedersen, A.; Jonsson, H.; Henkelman, G. EON: Software for long time simulations of atomic scale systems. Model. Simul. Mater. Sci. Eng. 2014, 22, 055002. http://www.theochem.org/EON. (41) Mason, R. Acta Crystallogr. 1964, 17, 547−555. (42) Cerius2, v. 4.7; Accelrys: San Diego, CA, 2002 (43) Karfunkel, H. R.; Gdanits, R. J. J. J. Comput. Chem. 1992, 13, 1171−1183. (44) Macrae, C. F.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Shields, G. P.; Taylor, R.; Towler, M.; van de Streek, J. J. Appl. Crystallogr. 2006, 39, 453−457. (45) Roux, M. V.; Temprado, M.; Chickos, J. S.; Nagano, Y. J. Phys. Chem. Ref. Data 2008, 37, 1855−1996. (46) Chakraborty, T.; Lim, E. C. J. Phys. Chem. 1993, 97, 11151− 11153. (47) Matsuoka, T.; Kosugi, K.; Hino, K.; Nishiguchi, M.; Ohashi, K.; Nishi, N.; Sekiya, H. J. Phys. Chem. A 1998, 102, 7598−7602. (48) Gonzalez, C.; Lim, E. C. Chem. Phys. Lett. 2000, 322, 382−388. (49) Gonzalez, C.; Lim, E. C. J. Phys. Chem. A 2000, 104, 2953− 2957. (50) Gonzalez, C.; Lim, E. C. J. Phys. Chem. A 2003, 107, 10105− 10110. (51) Podeszwa, R.; Szalewicz, K. Phys. Chem. Chem. Phys. 2008, 10, 2735−2746.

I


Evaluation of All-Atom Force Fields for Anthracene Crystal Growth

Recommend Documents