J. Phys. Chem. C 2008, 112, 9501–9509
9501
Estimate of Benzene-Triphenylene and Triphenylene-Triphenylene Interactions: A Topic Relevant to Columnar Discotic Liquid Crystals Giorgio Cinacchi*,† and Giacomo Prampolini‡ Dipartimento di Chimica, UniVersita di Pisa, Via Risorgimento 35, I-56126 Pisa, Italy ReceiVed: September 25, 2007; ReVised Manuscript ReceiVed: March 26, 2008
The interaction potential energy of several benzene-triphenylene and triphenylene-triphenylene dimer arrangements has been calculated by a Møller-Plesset second order perturbation theory with a suitably modified 6-31G* basis set. This approach, already tested and successfully employed on the benzene dimer, is presently shown to reproduce, with a level of accuracy comparable to other computationally affordable methods, the interaction potential energy of a few representative configurations of naphthalene dimer, for which high-level ab initio reference data exist. Thus, the method has been confidently applied to the interaction energy computation of the above-mentioned more complex dimers, for which high-level ab initio calculations are yet unfeasible. Comparisons have been made with the results obtained with an empirical force field and, as regards the benzene-triphenylene dimer, also with a dispersion corrected DFT method. The computed set of data may contribute to comprehending how aromatic interactions evolve with the ring size, a subject of general interest and particularly relevant to the field of columnar discotic liquid crystals, where the central cores of the mesogenic molecules often consist of triphenylene or triphenylene-derived units. In this case, attention has been paid to the determination of the favorable relative disposition of the monomers when the aromatic planes are parallel to each other, which is requisite information for their potential use as organic semiconductors. 1. Introduction Polycyclic aromatic hydrocarbons constitute a class of molecules whose members may be seen as formed by appropriately condensing benzenoid units. The first member of this particular series is naturally benzene, while graphene may be seen as the ultimate member. In between them, there exist polycyclic aromatic hydrocarbons having a disk-like shape, such as triphenylene, coronene, and hexabenzacoronene (Figure 1). These polycyclic aromatic disks are particularly interesting because they form the central core of columnar discogenic molecules. Typically, the structure of the latter consists of attaching to the border of one of these cores a number of pendant alkyl chains. These disk-like molecules frequently organize spontaneously to give rise to columnar liquid-crystalline phases. In such mesophases, the disk-like molecules are piled on top of each other, and the resulting columns are arranged in a regular lattice.1,2 Columnar discotic liquid crystals have attracted a great deal of recent attention because they are potentially very useful in a number of applications. All applications proposed up to now essentially rely on the possibility that these columnar liquid crystals behave as one-dimensional conductors.1,2 While lateral chains are certainly of importance in establishing the liquidcrystalline behavior and properties of these substances, charge transport properties are primarily determined by the aromatic character of the cores and the interactions among them. Therefore, it would be of importance to have a good estimate * Corresponding author. E-mail:
[email protected]. † Present affiliation: ARCES, Universita di Bologna, Via Vincenzo Toffano 2, I-40125 Bologna, Italy. Present contact address: Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy. ‡ E-mail: giacomo @ dcci.unipi.it.
of these interaction potential energy curves, especially in those pair configurations where the planes of the disks lie parallel to each other.3 Charge transport along the columns is very sensitive to the reciprocal position and orientation of two nearest-neighbor disks.4 It is conceivable that the face-to-face configurations, in which the two disks are exactly superimposed, are optimal for charge transport, but they could not correspond to the most favorable intermolecular interaction potential energies.4 Rather than the calculation of absolute interaction potential energy, which is a quantity very difficult to evaluate accurately, what is appearing more helpful for the design of columnar discotic liquid-crystalline semiconductors is a proper estimate of the relative interaction potential energies among the various dimer configurations and, hence, of the minimum energy arrangements. To this aim, different approaches may be employed: ab initio, semiempirical, and empirical methods, to mention them in order of computational feasibility. All of them must take into account the known fact that intermolecular interactions arise from a detailed balance of exchange-repulsion, electrostatic, induction, and dispersion forces.5 For aromatic hydrocarbon dimers, the formation of weakly bonded complexes is essentially due to dispersion interactions,6–10 a purely quantal electron correlation effect. For this reason, within the ab initio methods, one has to employ post-Hartree-Fock (HF) techniques. Among these, computational convenience would suggest to adopt Møller-Plesset second order perturbation (MP2) theory. Unfortunately, recent calculations on benzene dimer have shown that MP2 method with large basis sets overestimates the interaction energy by a nonnegligible factor.6,7,11–14 In the past 5 years it has indeed been shown 6–13 that highly correlated and computationally expensive methods, such as fourth order Møller-Plesset perturbation (MP4) theory or coupled cluster (CC) approaches, such as CCSD(T), must be used for accurately estimating
10.1021/jp0776917 CCC: $40.75 2008 American Chemical Society Published on Web 05/31/2008
9502 J. Phys. Chem. C, Vol. 112, No. 25, 2008
Cinacchi and Prampolini
Figure 1. Representative disk-like polyaromatic hydrocarbons: triphenylene (left), coronene (center), and hexabenzacoronene (right).
aromatic interaction energies. Furthermore, the adoption of large, diffuse, and polarizable basis sets is recommended.8–10 Consequently, such calculations are unfeasible, at least for standard computational resources, for polycyclic aromatic dimers much larger than benzene dimer, and whenever the aim is to construct entire interaction potential energy curves rather than an accurate estimate of the energy of just a few representative configurations. The semiempirical methods based on HF or density functional theory (DFT) are generally less expensive. However, in both cases dispersion corrections are indispensable. Electron correlation is not accounted for by the HF method, while none of the standard density functionals, based on the local density approximation, generalized gradient, and hybrid approximations, has been able to reproduce even benzene dimer interaction curves at a quantitative level.15,16 Previous works have shown how an empirical correction to the HF17 or DFT18,19 energy is capable of accounting for most of the dispersion interaction energy of aromatic dimers. The best results reported in ref 18, based on the BLYP density functional,20,21 underestimate the reference interaction energy values by =30%. These results are significantly ameliorated in ref 19, where a new density functional of the B97 type22 and new empirical corrections for a range of standard density functionals have been appropriately devised. In the present context, particularly relevant are the estimates of the binding energies of benzene and coronene to a graphene layer, and of the interaction between two graphene sheets very recently performed in ref 23 by this dispersion corrected DFT (DC-DFT) method. The computed values of interaction energy have been found in good agreement with the corresponding experimental estimates. One further interesting and potentially useful DC-DFT approach may be that based on optimizing atom-centered effective potentials.24 In particular, this method has been used in a very recent paper25 to calculate interaction potential energy data for disk-like polyaromatic hydrocarbons up to hexabenzacoronene. This approach is particularly attractive because it would permit, in the end, simulation of columnar liquid-crystalline materials by DFTbased molecular dynamics (MD) techniques,26 thus having access to not only structural and dynamic but also electronic properties. The method of refs 24 and 25 necessitates calibration of the correction on potential energy data already available and believed to be reasonably accurate. Unfortunately, this has not been yet the case. The first calibration performed on the benzene dimer24 overestimates the reference values7–10 by, on average, more than 50%,24,25 while the second calibration has been performed on the sole face-to-face benzene dimer configuration, and then is far from being satisfactorily completed.27 In addition, it would appear to be essential to include in the calibration process the correct long-range behavior of the dispersion forces, which scale with a distance R between the centers of the two interacting monomers as R-6. Among DFT-based techniques, it is finally worth mentioning that combination with symmetry adapted perturbation theory (SAPT) has been shown28–30 to yield results in good agreement with the most accurate ab initio values
for the benzene dimer. However, although its scaling properties with both the molecular dimensions and basis sets are better than CC methods,28 the DFT-SAPT approach seems still computationally too expensive for its straightforward application to the calculation of large polycyclic aromatic dimers. Turning to the less expensive empirical approaches, one route could be to extend literature force fields developed for benzene to polycyclic aromatic hydrocarbons. By employing force-field parameters, one can model atomistically large molecules and calculate the interaction between a pair of them instantaneously. The case of disk-like polycyclic aromatic hydrocarbons seems particularly fortunate: it would appear that the parameters of the optimized potentials for liquid simulations (OPLS), empirically determined for benzene,33 are transferable to the higher members of the series. The comparison of the potential energies calculated through transferred OPLS parameters with the available experimental data, namely adsorption energies of benzene and coronene molecules onto a graphene sheet and the interaction energies between two graphene sheets, is favorable.34 Moreover, it is here noticed that the OPLS potential energy curve for a pair of aromatic carbon sites is not so different from the same functions obtained from experimental data on fullerene-graphene and graphene-graphene interactions.35,36 Furthermore, the benzene dimer interaction energies computed via the OPLS parameters have been recently compared, for different configurations, with accurate quantum chemical results.37 It has turned out that the OPLS force field fails in predicting the overall shape of the parallel-displaced curves, but is capable of reproducing very well the benzene dimer interaction energy curve for the face-to-face geometries. This whole amount of facts suggests that the contact distance and well-depth OPLS parameters for a pair of aromatic carbon sites, determined originally for benzene, might be of a larger applicability. The possibility that they “interpolate” in a certain effective way between benzene and graphene, passing through polycyclic aromatic disks of intermediate size, is worth being explored. The quantum chemical intermolecular potential curves for the benzene dimer reported in ref 37 were obtained with a reduced computational effort by an appropriate choice of the method and the basis set. In fact, it has been reported 6,37–39 that an MP2 theory together with a modified 6-31G* basis set is able to provide results for the interaction potential energy of selected representative configurations of the benzene dimer in good agreement with the most accurate theoretical data available nowadays.7–10 This accord was achieved by appropriately calibrating the exponent of the Gaussian d functions centered on the carbon atoms.6,37–39,44 Notwithstanding it has been very recently criticized and perhaps too hastily dismissed as inadequate,29 the MP2 theory in conjuction with the modified 6-31G* basis set was previously used to reasonably calculate geometries and the interaction potential energies for selected configurations of not only the benzene dimer but also larger aromatic dimers up to the pyrene dimer.40–42 Noteworthy are also the reliable
Benzene- and Triphenylene-Triphenylene Dimers
J. Phys. Chem. C, Vol. 112, No. 25, 2008 9503
Figure 2. Selected arrangements for the benzene and naphthalene dimers.
estimates given by this method on the interaction potential energy of stacked nucleic acid base dimers.31,32 Finally, going to condensed matter properties, the force field generated from the above-mentioned MP2 set of dimer energies37 was proven to well reproduce a large collection of thermodynamic, structural, and dynamic data on benzene crystal and liquid in a wide range of temperatures at ambient pressure. In view of the seeming aptitude of this level of theory to treat aromatic hydrocarbon interactions, it has appeared sensible to extend the application of this MP2/6-31G* method to compute interaction potential energies for two of the next possible dimers in the polycyclic aromatic series, namely benzene-triphenylene and triphenylene-triphenylene. The computational burden has compelled a selection among the possible pair arrangements to investigate. Because of the columnar character of the discotic mesogenic materials, the planes of the two interacting molecules have been kept parallel to each other, that is, the face-to-face and parallel-displaced arrangements have been considered. A fortiori, the same computational resources presently available do not permit considering more complex dimers. The details of the calculations for benzene-triphenylene and triphenylene-triphenylene dimers are given in the next section. Section 3 contains the presentation and discussion of the results obtained, as well as the comparison with the predictions of the OPLS force field. Section 4 contains a conclusive summary. 2. Computational Details The optimization of all monomer geometries has been performed by the well-tested density functional B3LYP method43 with the correlation consistent basis set cc-pVDZ. The resulting optimized structures have been then employed for all quantum chemical and OPLS binding energy calculations. The quantum chemical energies have been computed in the supermolecule approach with the MP2 method. This method yields accurate interaction energies for the benzene dimer if carefully calibrated basis sets are employed.6,37–39 Here, the basis set is the standard 6-31G* modified for the low exponent of the Gaussian d functions centered on the carbon atoms: Rd ) 0.25 versus the usual values of 0.7-0.9, as already used in prior works.6,37–42 All energies have been computed using the counterpoise correction scheme, to take care of the basis set superposition energy error (BSSE).45 Up to now, although often
questioned, this is the accepted standard correction9 for the incompleteness of the basis set. Naphthalene-naphthalene and benzene-triphenylene dimer interaction energies have been also computed by the DC-DFT method of ref 19. In this case BLYP density functional20,21 with the TZVP basis set52 has been used, correcting for dispersion the always repulsive interaction energies with the empirical expression and parameters provided in ref 19. All MP2 and DFT calculations have been performed with the Gaussian 03 package.46 The OPLS interaction energies have been computed with the standard equation NA NB
EOPLS )
∑∑ i)1 j)1
{
[( ) ( ) ]}
qiqj σij + 4ij rij rij
12
-
σij rij
6
(1)
where NA and NB are the sum of all carbon and hydrogen atoms in monomers A and B, respectively. Point charges, qi, and Lennard-Jones parameters, i and σi, have been all taken from the all-atom OPLS force-field parameters optimized for benzene.33 Interaction parameters ij and σij were computed through standard OPLS mixing rules:33
ij ) √ij ;
σij ) √σiσj
(2)
3. Results and Discussion It is useful to commence this section by recalling the main results known for the benzene dimer. Representative arrangements of a pair of benzene molecules are the face-to-face (FF), parallel-displaced (PD), and T-shaped (TS) configurations (Figure 2). Table 1 shows the values of the optimized intermolecular distances and minimum interaction energies obtained for these arrangements employing several methods. The highest level calculations of refs 7–10 suggest that the PD and TS configurations are almost isoenergetic. These arrangements are supposed to be the most stable, with a binding energy of about 2.5-2.8 kcal/mol, whereas the FF arrangement has a binding energy ca. 1 kcal/mol lower. This trend, as well as that in the dimer geometries, is reproduced by the MP2/6-31G*(0.25) method with an average error of ca. 5%. The DC-DFT method of ref 19 is capable of reproducing benzene dimer optimal geometries and binding energies with
9504 J. Phys. Chem. C, Vol. 112, No. 25, 2008
Cinacchi and Prampolini
TABLE 1: Intermolecular Distances (R0 and R1, in Å) and Minimum Interaction Energies (E, in kcal/mol) for the Selected Arrangements of Figure 2 As Computed with the Indicated Methods FF
PD
TS
reference
method
R0
E
R0
R1
E
R0
E
7 8–10 37, 38 19 28 33
MP2/CCSD(T) MP2/CCSD(T) MP2 DC-DFT DFT-SAPT OPLS
3.8 4.0 3.9 3.9 3.9 3.8
-1.48 -1.81 -1.73 -1.77 -1.81 -1.69
3.5 3.6 3.5 3.5 3.5 3.5
1.8 1.6 1.8
-2.46 -2.78 -2.56 -2.75 -3.03 -2.07
5.0 5.0 5.0 4.9 5.0 5.1
-2.48 -2.74 -2.28 -2.99 -2.84 -2.15
1.8 2.7
TABLE 2: Comparison of Reference Data for the Naphthalene Dimer Interaction Potential Energy (in kcal/mol) As Obtained in Ref 48 with the Results of the MP2/6-31G*(0.25) and OPLS Methods for the Same Dimer Geometriesa
a
geometry
ref 48 (6-31G)
ref 48 (6-31G*)
MP2/6-31G*(0.25)
DC-DFT
OPLS
A B E F
-4.16 -5.73 -6.16 -4.46
-3.78 -5.28 -5.73 -4.34
-5.10 -6.85 -7.30 -4.78
-4.47 -6.07 -6.64 -6.37
-3.99 -4.51 -4.32 -2.61
The basis sets 6-31G and 6-31G* are those of medium size used in ref 48 for the CCSD(T) estimates at the basis set limit.
an average error of ca. 10%. The DFT-SAPT28 results also well agree with the ab initio reference values, except for the overestimation of the PD binding energy. The average error of the DFT-SAPT method is ca. 10%. Finally, the OPLS model, despite its empirical origin, is already known to reproduce the overall energetics of the benzene dimer interactions. Its average error with respect to the reference values is ca. 15%. This brief perusal of the results for the benzene dimer demonstrates that the MP2/6-31G*(0.25), as well as the empirical OPLS force field, is at least as promising as other computationally affordable methods proposed so far to deal with aromatic interactions. It is now of interest to compare the predictions of the quantum chemical and empirical methods also on the interaction potential energy of a few configurations of the naphthalene dimer, the other aromatic dimer for which high-level ab initio computations exist.48 These specific configurations are of the face-to-face, cross, parallel-displaced, and T-shaped type. They were labeled as A, B, E, and F, respectively, in ref 48 and are shown in Figure 2. The data of ref 48 were computed with the same method used before for the benzene dimer7 and provide an estimate of the CCSD(T) interaction potential energies at the basis set limit, with an accuracy of 0.5-0.7 kcal/mol. These estimates using the method of ref 48 depend on MP2 calculations performed with a medium-size basis set, such as 6-31G and 6-31G*. Despite being the second homologue of the polycyclic aromatic series, calculations on naphthalene dimer are already a highly demanding task. Indeed, quite a few studies can be found in the literature reporting geometry and energy data of this dimer to be compared with the most accurate results of ref 48. With the exception of ref 41, where the MP2/6-31G*(0.25) and MP2/ 6-31+G* methods were used, other MP2 predictions49–51 on the naphthalene dimer interaction energy, performed without adequately adapted basis sets, compare poorly (roughly 40% error, at least) with the reference values. Besides giving the reference data, Table 2 provides the values of the interaction potential energy for the four selected naphthalene dimer configurations computed with the MP2/6-31G*(0.25) and OPLS methods, as well as the DC-DFT method of ref 19. Notwithstanding the fact that the MP2/6-31G*(0.25) overestimates of ca. 20%, on average, the reference data, especially those evaluated making use of 6-31G* as a medium-size basis set, it is still able to satisfactorily reproduce the trend of the energy data, with the stabler dimers being ranked in the following order: E > B > A > F. That is, only the positions of
configurations A and F have been switched with respect to the order provided by the ab initio reference data. The DC-DFT method, although able to provide energies for the three configurations A, B, and E, having the planes of the two monomers parallel, within 10% of the reference data, substantially overestimates the binding energy of configuration F, corresponding to a T-shape arrangement. The DC-DFT method has proven not able to properly rank the four configurations investigated, with the F configuration being found almost the stablest. Finally, and, to a certain degree, interestingly, the OPLS force field underestimates by almost the same averaged percentage amount the CCSD(T) reference data. However, it cannot reproduce entirely the trend of energy data: the B configuration is found to be slightly stabler than the E configuration and the A configuration stabler than the F configuration. In light of these results on naphthalene dimer, the reliability of the MP2/631G*(0.25) method, based on the encouraging old data for the benzene dimer, is overall confirmed, if one adopts, as a criterion, the capability of a method to reproduce relative, rather than absolute, interaction energies. This section then continues with the presentation and discussion of the new results for the benzene-triphenylene and triphenylene-triphenylene interactions obtained with the MP2/ 6-31G*(0.25) method. They will be compared to those obtained applying the OPLS force field and, as regards the former dimer, also compared with DC-DFT estimates. The investigated arrangements of the dimer formed by the benzene and triphenylene molecules are depicted in Figure 3. The one labeled FF consists of placing the aromatic planes of the benzene and triphenylene molecules parallel to and the benzene molecule cofacial with the inner hexagon of the triphenylene molecule. The FF interaction potential energy curve is then generated by moving the benzene molecule along the axis perpendicular to the aromatic planes. This curve is shown in the left panel of Figure 4. The value of the minimum of the interaction potential energy is -7 kcal/mol and occurs at a distance R ) 3.4 Å. It may be interesting to note that at this distance the FF benzene dimer interaction energy is positive.37,38 In the benzene-triphenylene dimer, the decrease of the optimal plane separation derives from the contribution of the outer triphenylene rings, which are in a PD arrangement with respect to the benzene ring. The three PD curves considered for the benzene-triphenylene dimer are obtained by fixing the benzene coordinate on the above-mentioned axis at 3.4 Å, and displacing
Benzene- and Triphenylene-Triphenylene Dimers
Figure 3. Directions of the displacement applied to benzene (central bold hexagon) in the interaction with triphenylene (inner and outer hexagons). Direction FF is perpendicular to the benzene and triphenylene planes. When moving benzene along the PD1, PD2, and PD3 directions, the distance between the aromatic planes is fixed.
Figure 4. Interaction potential energy curves for the benzene-triphenylene dimer calculated with the MP2/6-31G*(0.25) (left panel), with the DCDFT (central panel), and with the OPLS force field (right panel). The directions of the R coordinate (FF, PD1, PD2, and PD3) are defined in Figure 3.
the phenyl ring along the three directions indicated in Figure 3. The corresponding interaction potential energy curves are also shown in the left panel of Figure 4. It is noticed that the three PD curves are superimposed up to a distance of 1 Å, where they all reach a minimum value of -7.5 kcal/mol, i.e., a bit lower than that of the FF configuration. This result confirms that displacing the benzene ring from being cofacial with a triphenylene hexagon is energetically advantageous. Nevertheless, the percentage difference between the FF and the three PD minima seems to diminish with the increase of the aromatic ring size. For distances beyond 1 Å, the PD1 curves rises most steeply, while the other two curves remain almost superimposed up to 6 Å. On the basis of the results known for the benzene dimer, the steeper rise of the PD1 curve is comprehensible since, along this direction, the benzene molecule sensibly departs from two of the outer hexagons of triphenylene, thus losing any favorable, PD-like, interaction with them, and moves toward the remaining outer hexagon. When the distance R reaches ca. 2.5 Å, i.e., twice the apothem of a benzene hexagon, the latter results cofacial with an outer triphenylene hexagon, and the interaction potential energy curve shows a local maximum, due to the less favorable FF-like arrangement. The successive local minimum is due again to a
J. Phys. Chem. C, Vol. 112, No. 25, 2008 9505 more favorable PD arrangement between the benzene molecule and the above-mentioned outer ring. The mean value of the minimum energy for the FF and PD arrangements is ca. 3.5 times larger than the corresponding quantity for the benzene dimer, although the number of atoms constituting the benzenetriphenylene supermolecule is 2.5 times larger than the number of atoms constituting the benzene dimer. By looking at the central panel of Figure 4, it appears that the DC-DFT curves behave with distance R in the same way as MP2/6-31G*(0.25) estimates. In fact, both methods agree in locating the FF minimum energy at 3.4 Å, in predicting the gain in binding energy by displacing the benzene monomer, and in the oscillatory behavior of the PD1 curve. As expected, the two methods do not agree in the absolute energy values, with the MP2 data being ∼2 kcal/mol lower than the DFT ones. Unfortunately, the lack of high-level, accurate reference values does not allow us to make absolute comparisons. The right panel of Figure 4 shows the analogous calculations performed employing the OPLS force field.33 This means that the triphenylene molecule has been modeled by assigning to the carbon and hydrogen atoms the same values of the well depth and contact distance given to these atoms in the OPLS model for benzene. For the inner carbon atoms a null charge has been assumed, while to the outer carbon atoms, as well as to the hydrogen atoms, the same value of the charge determined in the benzene case has been assigned.47Binding energies are sensibly smaller than those calculated with the MP2/631G*(0.25) method but comparable with the DC-DFT. Conversely, differently by quantum chemical findings, the most stable arrangement is the FF one, with an interaction energy of -5.0 kcal/mol at a distance of 3.6 Å. Note that in the evaluation of the PD interaction potential energy curves the two aromatic planes have been kept separated by this distance. The PD curves exhibit an overall smaller variance than those calculated with the MP2/6-31G*(0.25) and DC-DFT methods; in fact, the PD1 curve is not so different from the other two. They remain essentially constant up to a distance of 1 Å; that is, the OPLS force field predicts that the configurations obtained by these small displacements are almost isoenergetic with the FF arrangement. The displacements along the PD1 direction are anyway predicted to be the less favored, in agreement with the quantum chemical results. The other two PD curves are better resolved, with the PD3 displacement being the more favored. The value of -5 kcal/mol is ca. 2.5 times the mean value of the minimum energy for the FF and PD arrangements of the benzene dimer. Differently from the quantum chemical results, OPLS energies scale with the number of atoms of the dimer: this is expected since the OPLS-derived interaction energies are obtained as a sum of site-site contributions. Since the geometric information and the relative energies were almost identical between the two quantum mechanical methods, for computational saving the calculations for the largest dimer have been performed employing only one of the two methods, namely MP2/6-31G*(0.25). The computed interaction potential energy curves for the triphenylene-triphenylene dimer are shown in the left panel of Figure 5. These curves have been generated by taking two triphenylene molecules superimposed, rotating one with respect to the other of an angle γ around the common C3 axis, and separating them a distance RFF along this axis. The abscissa and, especially, the ordinate of the minimum energy decrease considerably, rotating one of the interacting molecules from 0° up to 30°, while it remains essentially constant after an additional rotation of 15°. The value of the minimum energy then gently increases after
9506 J. Phys. Chem. C, Vol. 112, No. 25, 2008
Cinacchi and Prampolini TABLE 3: Optimized Intermolecular Distances (R0, Å) and Minimum Energies (E0, kcal/mol) As Computed with the MP2/6-31G*(0.25) and OPLS Methods for the Face-to-Face Arrangement of the Three Aromatic Hydrocarbon Dimers MP2/6-31G*(0.25)
Figure 5. Interaction energy curves of the triphenylene-triphenylene dimer computed with the quantum chemical method and with the OPLS force field reported in the left and right panels, respectively. The angle γ defines the rotation of one triphenylene molecule around its Cˆ3 axis, while the RFF coordinate is the distance between the two planes. γ is null if the two molecules are superimposed at RFF ) 0.
an ulterior rotation of 15°, which makes the molecule rotated 60° in a star-like configuration. This behavior can again be qualitatively rationalized by decomposing the total interaction energy in a sum of hexagon-hexagon contributions. With γ ) 0°, the less favorable interaction energy arises from the four FF-like arrangements realized by the hexagons of the two triphenylene molecules. Conversely, departing from γ ) 0°, the contribution to the interaction energy of FF-like arrangements lessens in favor of more attractive PD-like arrangements. Figure 5 also shows how the differences among the curves tend to disappear rather quickly on increasing RFF: no difference is observed for distances larger than 4.5 Å. The mean minimum energy, ca. -16 kcal/mol, is 9.2 times more negative than the minimum energy for the FF arrangement of the benzene dimer, and ca. 2.3 more negative than the minimum energy for the FF arrangement of the benzene-triphenylene dimer. If these energy values were scaled with the number of atoms of the dimer, one would expect that the ratios of the mean minimum energy of the triphenylene dimer to the energy of the FF arrangement of the benzene-triphenylene and benzene dimers would be 2.5 and (2.5)2, respectively. While it appears that the first of these expectations is indeed realized, the second is not. This result, together with the already discussed findings about the scaling of the energy in the benzene-triphenylene dimer, suggest that the binding energies will scale with the number of atoms of the dimer as its size becomes larger. The analogous calculations performed employing the OPLS parameters are displayed in the right panel of Figure 5. In analogy with the benzene-triphenylene case, these empirical potential energy curves show a smaller value of the binding energy and a smaller variance with respect to the corresponding quantity evaluated with the quantum chemical method. There is little difference about the abscissa and the ordinate of the minimum energy among the curves corresponding to the different rotation angles. The values of these coordinates gently decrease by rotating one of the molecules from 0° to 60°. This could be expected, since the OPLS force field underestimates the energy gain in passing from the FF to a PD benzene dimer configuration. Moreover, the 45°-60° rotated star-like configurations are predicted to be the most stable, in contrast to the
OPLS
dimer
R0
E0
R0
E0
benzene-benzene benzene-triphenylene triphenylene-triphenylene
3.9 3.5 3.7
-1.73 -7.05 -12.78
3.8 3.6 3.7
-1.69 -4.99 -10.22
quantum chemical results which predict the 30°-45° rotated configurations as the most stable. The difference in the binding energy of these most favored configurations is significant: 11.8 kcal/mol for the OPLS force field versus 18 kcal/mol for the MP2/6-31G*(0.25). Expectedly, the binding energies of the triphenylene dimer are predicted by the OPLS force field to be ca. 2.5 larger than the corresponding quantity of the benzenetriphenylene dimer. The OPLS force field leads to interaction energies which scale with the number of atoms of the dimer. To reach a more comprehensive understanding of the evolution of the interaction potential energy with the size of the two monomers, as well as of the differences among the values of the quantum chemical and OPLS derived energies, it is useful to examine the FF curves for the benzene-benzene, benzenetriphenylene, and triphenylene-triphenylene dimers. Table 3 reports the value of the optimized distance and minimum energy as computed with the two methods. Figure 6 shows the three FF curves resolved in the two components relevant for the quantum chemical and OPLS methods. In the former case, the MP2 energy of interaction, EMP2 is decomposed in the HF contribution, EHF, and the dispersion contribution, defined as
Edisp ) EHF - EMP2
(3)
In the latter case, from the total interaction potential energy, EOPLS, the contribution due to sum of the attractive r-6 tails, E6, is extracted, while the rest of the energy has been labeled Eres. Qualitatively, the two sets of curves show a certain similarity. For instance, both methods agree in stating that the smallest value of the abscissa, for which the interaction potential energy is zero, is observed for the benzene-triphenylene dimer, while the largest value is found for the benzene-benzene dimer. In addition, the two methods exhibit a similar trend in the evolution of the respective dispersion contributions with dimer dimensions. For each dimer, the agreement between Edisp and E6 can be made almost quantitative by a slight horizontal shift of a few tenths of an angstrom. The most salient difference is found among the repulsive curves. The quantum chemical method predicts that the repulsive curve of the benzene-triphenylene dimer exhibits the less steep rise, whereas this happens for the benzene-benzene dimer among the OPLS curves. The relatively low values of the repulsive component of the quantum chemical curve for benzene-triphenylene are responsible for the emergence of an attractive energy well deeper than would be expected if the binding energies scaled with the number of atoms of the dimer. This may again be connected with the OPLS lack in reproducing the smaller repulsive contribution between two hexagons in the PD arrangement with respect to the FF arrangement. One further similarity between quantum chemical and OPLS derived sets of curves for the FF arrangement is that, once distance and energy values have been scaled with the respective abscissa and ordinate of the energy minimum, all curves lie
Benzene- and Triphenylene-Triphenylene Dimers
Figure 6. Ab initio and OPLS computed interaction potential energy curves of the benzene-benzene (B-B), benzene-triphenylene (B-T), and triphenylene-triphenylene (T-T) dimers in the face-to-face arrangement. In the left panel, the quantum chemical dispersion contributions (dotted lines) have been computed as the difference between total interaction energy (solid lines) and HF energy (dashed lines). In the right panel, dotted lines represent the dispersion contributions while the dashed lines are the remaining parts of the interaction energy computed by subtracting the dispersion contributions to the total interaction energies (solid lines).
J. Phys. Chem. C, Vol. 112, No. 25, 2008 9507 This simulation was based on the OPLS force field, which predicts almost equal minimum energies for two rotations of 45° and 60° (Figure 5), thus suggesting that condensed matter effects and the pendant aliphatic chains tend to favor smaller rotation angles. In an NMR experiment of benzene dissolved in triphenylenebased discotic liquid crystals, the nonmonotonic temperature dependence of the experimental deuterium splittings was interpreted in terms of a two-site model, where the solute was supposed to populate two different regions of the discotic solvent. In ref 56, the authors report that the intracolumn region between two triphenylene cores is favored at lower temperatures. Conversely, at higher temperatures, the benzene solute molecules are more likely to be found in the hydrocarbon matrix, formed by the pendant aliphatic tails. This picture has been partly contradicted by an atomistic simulation of benzene dissolved in a columnar phase of a triphenylene-based discogen.57 This numerical study, again based on the OPLS force field, has found that, at least at the single temperature considered, benzene is never able to populate the intracolumn region; rather it is always located within the hydrocarbon matrix. The nonmonotonic trend of the deuterium splittings with temperature can be explained by considering the nonisotropic character of the hydrocarbon matrix, where three different sites could be identified in which benzene adopts different orientations. The quantum chemical computed interaction energies of Figures 4 and 5 provide support to the simulation-derived scenario in that they predict that the intercalation of a benzene molecule within two triphenylene cores, presumably accompanied by a loss of entropy, is also an endothermic process. 4. Conclusions
Figure 7. Interaction energy master curve for the considered polyaromatic dimers in the face-to-face arrangement. The R0 and E0 values are reported in Table 3.
approximately on the same master curve (Figure 7). The existence of a master curve was previously reported in ref 25. For both methods, the curve that slightly departs from the other two is that corresponding to the benzene dimer. Having presented the results for the benzene-triphenylene and triphenylene-triphenylene interaction potential energies, comparisons with certain experimental and simulational observations are in order. The rotation angle corresponding to the minimum interaction energy of the triphenylene-triphenylene dimer predicted by the MP2/6-31G*(0.25) method, 30°-45°, is in agreement with the results of X-ray studies on triphenylenebased columnar discotic liquid crystals,53,54 where a value of 33°-45° is reported for two nearest-neighbor molecules in the same column. The rather broad distribution function of this rotation angle, calculated in an atomistic simulation of columnar hexakis(pentyloxy)triphenylene discogen, also peaks at 36°.55
The interaction potential energy for benzene-triphenylene and triphenylene-triphenylene dimers have been evaluated with quantum chemical and empirical methods. In the dimer configurations considered, the aromatic planes have been kept parallel. For the benzene molecule, two motions are considered with respect to the triphenylene molecule: the phenyl ring was either kept cofacial with the inner hexagon of the triphenylene and moved along the common symmetry axis or displaced in a plane separated from that of triphenylene by a fixed distance. The triphenylene dimer interaction energy was sampled by varying the distance between the aromatic planes at several rotation angles around the common C3 axis. By employing a MP2/6-31G*(0.25) approach, previously validated on the benzene dimer, the interaction potential energy values for selected configurations of the naphthalene dimer as well as the above-mentioned configurations have been computed and the results compared to those obtained with a DC-DFT method and OPLS empirical force field. The overall good agreement between all methods in interaction potential energies observed for benzene dimer deteriorates for the naphthalene dimer and no longer exists for the larger dimers investigated in the present work. As regards the OPLS force field, the origin of these differences probably resides in the underestimation, with respect to the quantum chemical methods, of the contributions to the interaction potential energy coming from two phenyl hexagons in a PD-like configuration. Hence, the methods do not reinforce reciprocally, nor can the quantum mechanical methods be used to support the hypothesis that the OPLS optimized for benzene is safely transferable to polycyclic aromatic dimers of intermediate size. While the quantum chemical and OPLS methods reproduce the reference data for benzene dimer energetics, the former
9508 J. Phys. Chem. C, Vol. 112, No. 25, 2008 overestimates and the latter underestimates the binding energies of the naphthalene dimer. However, the MP2/6-31G*(0.25) method satisfactorily reproduces the order of the stabler configurations, while DC-DFT provides too attractive energies for the T-shaped dimer and the OPLS method is at least able to predict the parallel-displaced, E, configuration significantly stabler than the T-shape, F, configuration. The available experimental data do not contradict either method. However, their paucity and very indirect nature do not permit truly singling out which one performs better. On the basis of the results for benzene and naphthalene dimers and all available experimental data, the quantum chemical methods should be given more credit as far as the prediction of dimer optimal geometries and relative energetics is concerned, while the OPLS force field could be still considered a valid source of qualitative information for aromatic hydrocarbon interaction energy, especially when large aromatic dimers are considered. If the absolute energetics of the interaction is concerned, the fact that, on one hand, the quantum chemical and OPLS methods give comparable predictions in the case of benzene dimer, while, on the other hand, the MP2/6-31G*(0.25) overestimates the interaction potential energies of naphthalene dimer by almost the same amount as the latter underestimates them may suggest that the two methods may be taken, respectively, as upper and lower bounds for binding energies of aromatic dimers. In this sense the two methods might be used synergistically to generate a set of interaction potential energy data which serves different purposes. This set can be used to parametrize coarse-grained58–61 or atomistic55,57,62 pair potential for use in computer simulations of columnar discotic liquid crystals to cite just one application in the field which motivated this work. More generally, the above-mentioned set of data may contribute to the comprehension of the aromatic-aromatic interactions, a long-debated subject of great importance for chemistry, biochemistry,63–65 and nanotechnologies1,66 with repercussions on other diverse fields such as astrophysics,67 astrobiology,68 and medicine.69 References and Notes (1) Bushby, R.; Lozman, O. Curr. Opin. Colloid Interface Sci. 2002, 7, 343. (2) Laschat, S.; Baro, A.; Steinke, N.; Giesselmann, F.; Hagele, C.; Scalia, G.; Judele, R.; Kapatsina, E.; Sauer, S.; Schreivogel, A.; Tosoni, M. Angew. Chem., Int. Ed. 2007, 46, 4832. (3) Mulder, F. M.; Stride, J.; Picken, S. J.; Kouwer, P. H. J.; de Haas, M. P.; Siebbeles, L. D. A.; Kearley, G. J. J. Am. Chem. Soc. 2003, 125, 3860. (4) Lemaur, V.; Da Silva Filho, D. A.; Coropceanu, V.; Lehmann, M.; Geerts, Y.; Piris, J.; Debije, M. G.; Van de Craats, A. M.; Senthilkumar, K.; Siebbeles, L. D. A.; Warman, J. M.; Bredas, J. L.; Cornil, J. J. Am. Chem. Soc. 2004, 126, 3271. (5) Stone, A. J. The Theory of Intermolecular Forces; Oxford University Press: Oxford, 1996. (6) Hobza, P.; Selzle, H. L.; Schlag, E. W. J. Phys. Chem. 1996, 100, 18790. (7) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. J. Am. Chem. Soc. 2002, 124, 104. (8) Sinnokrot, M.; Valeev, E.; Sherrill, C. J. Am. Chem. Soc. 2002, 124, 10887. (9) Sinnokrot, M.; Sherrill, C. J. Phys. Chem. A 2004, 108, 10200. (10) Sinnokrot, M.; Sherrill, C. J. Phys. Chem. A 2006, 110, 10656. (11) Jaffe, R. L.; Smith, G. D. J. Chem. Phys. 1996, 105, 2780. (12) Tsuzuki, S.; Uchimaru, T.; Matsamura, K.; Mikami, M.; Tanabe, K. Chem. Phys. Lett. 2000, 319, 547. (13) Tsuzuki, S.; Uchimaru, T.; Sugawarea, K.; Mikami, M. J. Chem. Phys. 2002, 117, 11216. (14) Donchev, A.; Galkin, N.; Pereyaslavets, L.; Tarasov, V. J. Chem. Phys. 2006, 125, 244107. (15) Meijer, E.; Sprik, M. J. Chem. Phys. 1996, 105, 8684. (16) Tsuzuki, S.; Lu¨thi, H. J. Chem. Phys. 2001, 114, 3949.
Cinacchi and Prampolini (17) Gonzalez, C; Lim, E. C. J. Phys. Chem. A 2003, 107, 10105. (18) Grimme, S. J. Comput. Chem. 2004, 25, 1463. (19) Grimme, S. J. Comput. Chem. 2006, 27, 1787. (20) Becke, A. Phys. ReV. A 1988, 38, 3098. (21) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (22) Becke, A. J. Chem. Phys. 1997, 107, 8554. (23) Grimme, S; Mu¨ck-Lichtenfeld, C.; Antony, J. J. Phys. Chem. C 2007, 111, 11199. (24) von Lilienfeld, O.; Tavernelli, I.; Rothlisberger, U.; Sebastiani, D. Phys. ReV. Lett. 2004, 93, 153004. (25) von Lilienfeld, O.; Andrienko, D. J. Chem. Phys. 2006, 124, 054307. (26) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (27) Lin, I. C.; Coutinho-Neto, M. D.; Felsenheimer, C.; von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Phys. ReV. B 2007, 75, 205131. (28) Hasselmann, A.; Jansen, G. J. Chem. Phys. 2005, 122, 014103. (29) Podeszwa, R.; Bukowski, R.; Szalewicz, K. J. Phys. Chem. A 2006, 110, 10345. (30) Tekin, A.; Jansen, G. Phys. Chem. Chem. Phys. 2007, 9, 1680. (31) (a) Hobza, P; Sˇponer, J. J. Am. Chem. Soc. 2002, 124, 11802. (b) Jurecka, P.; Hobza, P. J. Am. Chem. Soc. 2003, 125, 15608. (32) Hobza, P.; Zahradnik, R; Mu¨ller-Dethlefs, K. Collect. Czech. Chem. Commun. 2006, 71, 443. (33) Jorgensen, W.; Severance, D. J. Am. Chem. Soc. 1990, 112, 4768. (34) Cinacchi, G. J. Chem. Phys. 2006, 125, 057101. (35) Girifalco, R. L. L. A.; Hodak, M. Phys. ReV. B 2000, 62, 13104. (36) Ulbricht, T. H. H.; Moos, G. Phys. ReV. Lett. 2003, 90, 095501. (37) Cacelli, I.; Cinacchi, G.; Prampolini, G.; Tani, A. J. Am. Chem. Soc. 2004, 126, 14278. (38) Cacelli, I.; Cinacchi, G.; Prampolini, G.; Tani, A. J. Chem. Phys. 2004, 120, 3648. (39) Amovilli, C.; Cacelli, I.; Campanile, S.; Prampolini, G. J. Chem. Phys. 2002, 117, 3003. (40) Lee, N. K.; Park, S.; Kim, S. K. J. Chem. Phys. 2002, 116, 7902. (41) Lee, N. K.; Park, S.; Kim, S. K. J. Chem. Phys. 2002, 116, 7910. (42) Lee, N. K.; Kim, S. K. J. Chem. Phys. 2005, 122, 031102. (43) Becke, A. J. Chem. Phys. 1993, 98, 5648. (44) Note that a typographical error occurred in both refs 37 and 38, whereRd was erroneously put equal to 0.28 instead of 0.25, the correct value. Nevertheless, no appreciable differences in the dimer geometry and binding energy arise upon the use of these so slightly different exponent values. (45) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (46) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision 0.1; Gaussian, Inc.: Pittsburgh, PA, 2003. (47) Analogous calculations have been performed using, for triphenylene, the partial charges obtained from the EPS fit of the monomer charge distribution as results from the B3LYP/cc-pVDZ computation. No appreciable differences have been observed in the calculated binding energies. (48) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M. J. Chem. Phys. 2004, 120, 647. (49) Gonzalez, C.; Lim, E. C. J. Phys. Chem. A 2000, 104, 2953. (50) Walsh, T. R. Chem. Phys. Lett. 2002, 363, 45. (51) Fomine, S.; Tlenkopatchev, M.; Martinez, S.; Fomina, L. J. Phys. Chem. A 2002, 106, 3941. (52) Schaefer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829. (53) Levelut, A. M. J. Phys. Lett. 1979, 40, 81. (54) Fontes, E.; Heiney, P. A.; de Jeu, W. H. Phys. ReV. Lett. 1988, 61, 1202. (55) Cinacchi, G.; Colle, R.; Tani, A. J. Phys. Chem. B 2004, 108, 7969. (56) Goldfarb, D.; Luz, Z.; Zimmermann, H. J. Phys. (Paris) 1982, 43, 1255. (57) Cinacchi, G. J. Phys. Chem. B 2005, 109, 8125. (58) Bates, M. A.; Luckhurst, G. R. J. Chem. Phys. 1996, 104, 6696. (59) Zewdie, H. B. Phys. ReV. E 1998, 57, 1793. (60) Cinacchi, G.; Tani, A. J. Chem. Phys. 2002, 117, 11388. (61) Caprion, D.; Bellier-Castella, L.; Ryckaert, J. P. Phys. ReV. E 2003, 67, 041703.
Benzene- and Triphenylene-Triphenylene Dimers (62) Andrienko, D.; Marcon, V.; Kremer, K. J. Chem. Phys. 2006, 125, 124902. (63) McGaughey, G.; Gagne´, M.; Rappe´, A. J. Biol. Chem. 1998, 273, 15458. (64) Hobza, P.; Sponer, J. J. Am. Chem. Soc. 2002, 124, 11802. (65) Meyer, E.; Castellano, R.; Diederich, F. Angew. Chem., Int. Ed. 2003, 42, 1210. (66) van de Craats, A.; Warman, J.; Mu¨llen, K.; Geerts, Y.; Brand, J. AdV. Mater. 1998, 10, 36.
J. Phys. Chem. C, Vol. 112, No. 25, 2008 9509 (67) Henning, T.; Salama, F. Science 1998, 282, 2204. (68) Ehrenfreund, P.; Rasmussen, S.; Cleavesand, J.; Chen, L. Astrobiology 2006, 6, 490. (69) Perera, F. P. Science 1997, 278, 1068.
JP0776917
