Binding Interactions in Dimers of Phenalenyl and ... - ACS Publications

Apr 1, 2013 - Phenalenyl, an open-shell neutral radical that can form both π-stacked dimers and conducting molecular crystals, has gained attention f...
2 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCA

Binding Interactions in Dimers of Phenalenyl and Closed-Shell Analogues Brian Kolb,† Miklos Kertesz,‡ and T. Thonhauser*,† †

Department of Physics, Wake Forest University, Winston-Salem, North Carolina 27109, United States Department of Chemistry, Georgetown University, Washington, D.C. 20057, United States

‡

S Supporting Information *

ABSTRACT: Phenalenyl, an open-shell neutral radical that can form both πstacked dimers and conducting molecular crystals, has gained attention for its interesting and potentially useful electrical and magnetic properties. The properties of this complex physical system are fairly well understood, making it an ideal testing ground for the newly developed van der Waals density functional (vdW-DF). We invoke a simple approximation, allowing the vdWDF to be used within spin-polarized density functional theory and test this approximation on the π-stacked phenalenyl dimer. The results indicate that the vdW-DF is capable of qualitatively describing the interaction between two neutral radicals in the πstacked configuration, producing, in line with experiment, binding distances that are significantly below the sum of the van der Waals radii. This is a nontypical distance range where most other theories fail. We then investigate two hypothetical closed-shell analogues of this dimer, one formed by replacing the central carbon of phenalenyl with a nitrogen atom and the other formed by replacing the central carbon with a boron atom. In these cases, relatively strong interaction energies are obtained at more typical equilibrium distances for van der Waals dimers. The nitrogen-substituted dimer shows an unexpected rotational barrier that is dictated by the electronic kinetic energy within the system. The torsional curve of the boron-substituted dimer also exhibits a rotational barrier, but this is found to disappear when exact exchange is used in place of a local or semilocal exchange functional.

■

INTRODUCTION Interest in phenalenyl (C-Phy) and various derivatives has grown recently both because of their unusual binding characteristics1 and because they can form molecular crystals with tunable electric, magnetic, and optical properties.2 Even the isolated phenalenyl molecule has been shown to exhibit charge-transfer capabilities that may allow it to be used effectively within molecular circuits.3,4 When two C-Phy molecules bind in a π-stacked fashion, the contact distances can be surprisingly short, significantly under the sum of the van der Waals radii.1 This distance can be made even shorter by judicious placement of electron withdrawing substituents. This tunable distance is important because a recent report shows that the dimer can be made to conduct effectively if the contact distance can be made short enough.5,6 Overall, the short binding distance of this systemtogether with its orientational specificitymake it very intriguing, as it is quite different from the typical π−π stacking case. It has long been known that, when C-Phy molecules are covalently linked, they can be crystallized to form solids with varying conductivities, depending on the precise form of the linker and the side chains present.7−15 Even relatively minor changes in the attached ligands can effect drastic changes in the molecular crystal packing, which in turn helps determine the crystal conductivity.12 However, a systematic computational means to study the link between ligands and crystal packing/ conductivity has yet to be developed. The recently described, fully nonlocal van der Waals density functional16,17 (vdW-DF) is well tested in numerous systems © 2013 American Chemical Society

where weak van der Waals interactions play an important role,18 including simple dimers,19,20 physi-adsorbed molecules,21−23 DNA,24,25 and molecular crystals.26,27 Indeed, its use has demonstrated that some systems for which van der Waals interactions were thought to be relatively unimportant actually require their inclusion, sometimes to get even qualitative behavior correctly.28−31 The van der Waals density functional approach writes the exchange−correlation energy as (semi)local ExcvdW‐DF = Exc + Excnonlocal 1 (semi)local = Exc + n( r1⃗) Φ( r1⃗ , r2⃗ ) n( r2⃗ ) d r1⃗ d r2⃗ 2

∫∫

(1)

where Φ(r1⃗ ,r2⃗ ) is a nonlocal exchange−correlation kernel that links the charge densities n at r1⃗ and r2⃗ . For the (semi)local exchange and correlation, E(semi)local , in principle a wide variety xc of options exists, but typically ErevPBE + ELDA is used, employing x c standard functionals and thus guaranteeing good performance in many standard cases. The form of the integral in eq 1 was first derived by Dion et al.16 and then efficiently written as a convolution by Román-Pérez and Soler.32 Unlike the DFT-D approach of Grimme,33 the functionals of Truhlar et al.,34,35 and range-separated hybrid functionals,36−38 the vdW-DF approach does not rely on fitting to experimental data, nor does it require the use of orbitals rather than densities. A detailed discussion Received: September 25, 2012 Revised: March 29, 2013 Published: April 1, 2013 3642

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649

The Journal of Physical Chemistry A

Article

dimer bonding orbital with electron density localized between eclipsed α sites.1 This type of binding has a strong angular dependence in addition to the expected distance dependence. Addition of suitable electron withdrawing groups at either the α or β sites has been shown to reduce the binding distance even further.6 This tunable binding distance is inextricably linked to the electron transport capabilities of the dimer system. Replacing the central carbon atom with a nitrogen atom (forming N-Phy) causes the highest-occupied molecular orbital (HOMO) to be doubly occupied and destroys the radical nature of the molecule. Without the ability to form the 2electron/multicenter bond, the binding distance grows and the electron transport properties of the new molecular species likely diminish. Although one might expect the N-Phy2 dimer formed from this species to be governed predominantly by πstacking and van der Waals interactions, we show in this work that there is a considerable rotational barrier set up by the electronic kinetic energy. Both electrostatics and exchange− correlation oppose this barrier but are overpowered by the kinetic energy contribution to the energy. Replacing the central carbon of C-Phy with a boron atom (forming B-Phy) also has the effect of destroying the radical nature of the dimer. In this case the dimer behaves as expected, with only a general preference for the staggered conformation. However, it was found that the local Slater exchange that is a fundamental ingredient in so many modern density functionals yields qualitatively incorrect behavior, causing a spurious apparent rotational barrier considerably weaker than in the N-Phy2 case. In this work we (i) develop an approximation to allow spinpolarized calculations using the recently developed van der Waals density functional (vdW-DF) and demonstrate the applicability and accuracy of the approach in C-Phy2, an extremely complex and delicate physical system, (ii) describe an unexpected rotational barrier governed by electronic kinetic energy in the N-Phy2 system, and (iii) describe the B-Phy2 system, where use of approximate exchange functionals leads to qualitatively wrong behavior for the rotational energy.

and comparison of all these methods can be found in a recent review.39 vdW-DF is becoming increasingly popular and many systems have already been successfully studied with it. This prompted us to consider the more challenging phenalenyl systems as a perfect testing ground to explore the limitations of vdW-DF and we believe that the information gained here will help to improve the functional further. The functional would be well-suited to the problem of studying molecular crystals of C-Phy with an aim toward understanding how side-chain and linker properties determine crystal packing and subsequent conductivity. Unfortunately, two problems exist that preclude this use. First, as will be discussed next, phenalenyl is a neutral radical, meaning spinpolarized DFT must be used to study it. At present, no spinpolarized version of the vdW-DF exists. Second, the radical nature of phenalenyl makes it an extremely difficult molecule to study computationally, requiring special methods and care to obtain properly converged results. The vdW-DF is relatively untested on such difficult computational ground. In this work, an approximate spin-polarized implementation of the vdW-DF was devised and used. To test the validity of the results, secondorder Møller−Plesset perturbation theory (MP2) calculations were also performed for comparison purposes. MP2 is not without its deficiencies and it should be stressed here that the method is not being used as an absolute reference, but rather as a completely independent method with which to compare the vdW-DF results. Coupled cluster calculations would surely serve better here, but the number of atoms involved (44 atoms, 122 valence electrons) makes use of coupled cluster techniques infeasible for these systems. MP2 is the next best choice; its relatively good scaling makes it affordable to consistently apply to all the calculations performed herein. Phenalenyl, shown in Figure 1a, is a neutral molecule that can be visualized as a truncated graphene sheet consisting of

■

METHODS All calculations were carried out with both Møller−Plesset perturbation theory (MP2) and density functional theory utilizing the van der Waals density functional (vdW-DF) of Dion et al.16 All monomer coordinates were fully relaxed within each method. Dissociation and rotation curves were calculated by bringing the relaxed monomers together rigidly in the desired geometry. The rotation curves were calculated at a separation of 3.3 Å, i.e., the experimental distance in the C-Phy2 system,1 and dissociation curves were calculated for each system at its minimum-energy rotation angle. All energies are reported as kcal/mol of dimer. The reference state for all rotational calculations is the fully eclipsed dimer conformation at 0°. MP2 calculations were carried out using the Gaussian 03 quantum-chemistry package.40 Except as otherwise noted, all production calculations utilized the 6-31++G** basis set with an SCF convergence criterion of 1 × 10−6 Hartree (Eh). Use of a larger basis set would certainly have been preferable, but practical considerations limited the present work to the chosen basis set, which is acceptable for systems of this size. Monomer coordinates for each system were fully optimized until all forces were less than 1.5 × 10−5 Eh/Å. All dissociation curves include the counterpoise correction for the basis set superposition error

Figure 1. (a) Phenalenyl molecule. The labels α and β mark two symmetry inequivalent sites that play an important role in the binding. (b) Plot of the calculated singly occupied molecular orbital (SOMO) showing the delocalization of the unpaired spin over the α sites in CPhy.

three aromatic rings terminated by hydrogen. The molecule has 3-fold rotational symmetry with four inequivalent carbon sites. The sites labeled α and β in Figure 1 are particularly important for the binding properties of the dimer. C-Phy has an odd number of electrons making it, necessarily, an open-shell radical with one unpaired electron residing in a singly occupied molecular orbital (SOMO). As shown in Figure 1, this SOMO is delocalized over the α sites of the molecule. When two C-Phy monomers stack face-on, they form a π-dimer notable for its extremely short binding distance. The interactions holding the dimer together go beyond standard π-stacking and dispersion interactions. The close contact between the monomers (3.20− 3.32 Å, less than the sum of van der Waals radii) is a consequence of a strong 2-electron/multicenter bond formed when the two SOMOs overlap to create a doubly occupied 3643

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649

The Journal of Physical Chemistry A

Article

(BSSE). By taking the difference in BSSE between the 0° and 60° rotations, we estimated the BSSE contribution to the rotational energy scans to only be 0.1−0.2 kcal/mol, so it is not included in rotational energy scans. Restricted open-shell Hartree−Fock reference states were used for both B-Phy2 and N-Phy2, whereas a CASSCF(2,2) reference wave function was used for C-Phy2. Triplet states were approximated by requiring two more α electrons than β electrons in a spinunrestricted calculation. The phenalenyl system considered here is difficult to study computationally and application of the methods described above to any open-shell radical has to be interpreted with care. van der Waals density functional calculations were carried out using the Quantum-Espresso package.41 A plane-wave cutoff of 30 Eh, in conjunction with an SCF convergence of 5 × 10−11 Eh and ultrasoft pseudopotentials yields sufficient precision. Monomer relaxations were stopped when each component of every nuclear force was less than 1 × 10−4 Eh/ Å. Because a plane-wave basis presumes periodicity, the choice of box size is critical in calculations of isolated molecules. In all calculations presented herein, a box size of 16.9 Å was used to prevent significant interaction of periodic images. As in the CASSCF(2,2)-MP2 calculations, C-Phy2 requires special care. All C-Phy calculations were treated using spin-polarized density functional theory. A truly spin-polarized version of the vdW-DF does not exist yet. To include spin-polarization and the vdWDF simultaneously, it was assumed that the nonlocal piece of the exchange−correlation functional does not depend explicitly on spin (all local parts are treated properly with LSDA). Under this assumption, the density used for the evaluation of the nonlocal part of vdW-DF, Enl[n(r)⃗ ], is just the total density. That is, n(r)⃗ = n↑(r)⃗ + n↓(r)⃗ , where n↑↓(r)⃗ denotes the up or down spin charge density. The correct nonlocal Kohn−Sham potential for each spin channel is then just the total nonlocal Kohn−Sham potential, so that νnl↑ = νnl↓ = νnl. It should be stressed that, although this approach is not as good as a fully spin-polarized treatmentas, for example, in the functional of Vydrov and Van Voorhis42,43the severity of the approximation is relatively mild. To test the severity of the approximation, spin-unpolarized calculations were carried out for C-Phy2 and the nonlocal piece of the energy was compared to the nonlocal piece from our approximate spin-polarized formulation. The calculations were carried out for a dissociation curve with monomer separations between 2.5 and 6 Å. The maximum difference between the two terms was 0.25 kcal/mol, leading us to believe that the approximation made in our spinpolarized treatment is justified. Both B-Phy2 and N-Phy2 were treated with spin-unpolarized density functional theory. Further calculations were carried out in the NWchem44 package for the B-Phy2 system. These calculations utilized spinrestricted density functional theory with a cc-pVDZ basis set. The calculations were performed with exact exchange, local Slater exchange45 and the gradient corrected functional of Perdew, Burke, and Ernzerhof.46

Figure 2. Rigid rotational energy scan for the C-Phy2 dimer as calculated by CASSCF(2,2)-MP2 and vdW-DF at a separation of 3.3 Å. The two methods are in quantitative agreement and both show a distinct minimum at 60°, in agreement with experiment. The energy curve is periodic in 120° owing to the 3-fold symmetry of the monomers. The reference state is the fully eclipsed dimer conformation at 0°.

0°, the molecules are fully eclipsed, allowing for full overlap of the monomer SOMOs. The resulting dimer bonding orbital is doubly occupied, and the dimer antibonding orbital is fully unoccupied. As the rotation angle is increased, the 2-electron/ multicenter bond is strained, resulting in an energy barrier to rotation. At 60°, the β carbons are staggered, but the α carbons again eclipse. The 2-electron/multicenter bond is restored, but with the β carbons and their attached hydrogen atoms having more space, this geometry is preferred. It is interesting to investigate the effect on the rotation curve of dissallowing electron pairing in the dimer HOMO. This can be accomplished in the post-Hartree−Fock methods by enforcing a spin multiplicity of 3, and in the DFT method by forcing a total spin of 2in each case resulting in 2 unpaired electrons, one in each of the dimer bonding and antibonding orbitals. This effectively breaks the 2-electron/multicenter bond, causing the molecular rotation to be governed primarily by weak van der Waals interactions. The resulting rotation curve for this situation is shown in Figure 3. As can be seen in

Figure 3. Rigid rotation curve of C-Phy2 at a separation of 3.3 Å, made by requiring two unpaired electrons in the system. In this case, the 2electron/multicenter bond is broken and, as expected, the system shows little preference for any particular angle. The fully eclipsed conformation is still unfavored because of steric effects.

the figure, the rotational characteristics are qualitatively changed in this situation. The rotational preference for a 60° rotation angle is broken and, as expected, the molecule shows little preference for any particular rotation angle other than the expected dislike of the eclipsed conformation. The binding energy as a function of monomer separation is shown in Figure 4. In contrast to the rotation curve there is quantitative disagreement between the two methods in this case. The MP2 curve shows stronger binding at a shorter distance than the vdW-DF curve. It should be noted that the CASSCF(2,2)-MP2 binding energy calculated here is in excellent agreement with recent calculations by Mota, Miller, and Novoa48 but disagrees with the earlier calculations of Small

■

RESULTS AND DISCUSSION C-Phy2. Figure 2 shows the energy required to rotate the monomers of C-Phy2 relative to each other. Both CASSCF(2,2)-MP2 and vdW-DF show an energy minimum at a rotation angle of 60°, which is consistent with the D3h symmetry found experimentally for the tri-tert-butyl-substituted dimer. The preference for a rotation of 60° is a consequence of the 2-electron/multicenter bonding present in the system. At 3644

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649

The Journal of Physical Chemistry A

Article

in energy for all relevant geometries, the DFT calculations herein held the total magnetization fixed at zero. However, one may still define an integrated local magnetization as ∫ |n↑(r)⃗ − n↓(r)⃗ | d3r, which is a measure of the total spatial difference between the up and down spin densities. This quantity is plotted as a function of monomer separation in Figure 5a. As

Figure 4. Binding energy as a function of monomer separation for each of the three systems studied in this work. In each case, results are shown for the vdW-DF and MP2 (CASSCF(2,2)-MP2 in the case of C-Phy2) with the monomer rotation set to the calculated optimum value for that system (Table 1). Experimental values for enthalpy47 and binding distance1 depend somewhat on solvent and substituents, respectively. The values presented here are for tert-butyl substituted CPhy2 in hexane (upper value) and dichloromethane (lower value). The binding energies and equilibrium separations for all three systems with both methods are given in Table 1.

Figure 5. Integrated local magnetization in C-Phy2, defined as ∫ |n↑(r⃗) − n↓(r)⃗ | d3r, in Bohr magnetons as a function of (a) monomer separation at a rotation angle of 60° and (b) rotation angle at a separation of 3.3 Å. The vertical dashed line in (a) shows the optimized dimer separation for vdW-DF of 3.31 Å.

can be seen in the figure, spontaneous magnetization develops when the monomers are separated by about 3.3 Å, reaching full saturation at a separation around 4.5 Å. Figure 5b shows the same quantity plotted as a function of monomer rotation at a separation of 3.3 Å. At this separation, magnetization develops when the monomers are rotated beyond about 5°, reaching a peak at 30°. The appearance of spontaneous magnetization is indicative of the point at which spin-polarization effects signifying diradicaloid character become important. When no local magnetization is present, a spin-unpolarized calculation will yield correct results. As magnetization increases, spin up and down densities become more spatially distinct. Experiments suggest that magnetization exists already at the equilibrium separation of C-Phy2.1 Our calculations, as evident from Figure 5a, show that magnetization is developing close to the optimum separation, in line with experiment. These results for C-Phy2 are well-known,1,6,49 and it is encouraging to see that the relatively inexpensive vdW-DF approach captures the physics quite accurately, even in this very difficult chemical situation. Even though the spin-polarization approach used is approximate, it seems to capture the nature of the C-Phy2 system quite well. All the experimentally observed characteristics are well reproduced, and the method agrees qualitatively with MP2 on the behavior of the dissociation curve and nearly quantitatively with MP2 on the characteristics of the rotation curve, both with and without the 2-electron/multicenter bond. N-Phy2. The bonding in C-Phy2 is largely a result of one unpaired electron in each monomer filling the bonding molecular orbital created by SOMO−SOMO overlap to form a 2-electron/multicenter bond. It is interesting, then, to see what happens when an additional electron is added to each monomer to make a closed-shell system. This can be accomplished by replacing the central carbon atom of C-Phy with a nitrogen to form N-Phy. The extra electron in nitrogen as compared to carbon means the HOMO in N-Phy is doubly occupied and the 2-electron/multicenter bond cannot be formed. In this hypothetical situation, one might expect a priori that the N-Phy2 system would show a fairly flat rotation curve, similar to the C-Phy2 case with the 2-electron/

et al.49 The vdW-DF predicts a binding energy of 14.6 kcal/mol at a separation of 3.31 Å. MP2 calculations, by contrast, predict an optimum binding of 30.9 kcal/mol at a separation of 2.94 Å. Note that these numbers do not include the zero-point correction. This disagreement is expected, however, because it is well documented that (i) MP2 tends to overbind, especially in systems with significant π-stacking interactions,50 and (ii) vdW-DF tends to overestimate the binding distance for all systems.19 When these considerations are taken into account, the agreement between the two methods becomes fairly good and both of these values are in reasonable agreement with the experimental result of 7.5−9.5 kcal/mol. It is likely that the MP2 numbers are affected somewhat by an insufficient basis set. Unfortunately, the size of the system (44 atoms) makes using an even larger basis set difficult. It is also worth noting here that the popular semilocal exchange−correlation functional revPBE, the popular functional (U)B3LYP,51 and the Hartree−Fock approach give no binding whatsoever for the dimer, which demonstrates the need for a careful treatment of correlation to correctly account for the binding interactions in this system. It is interesting to investigate the appearance of spontaneous magnetization in C-Phy2 as a function of monomer separation and rotation. The total magnetization of C-Phy2 must be either 2 or 0, corresponding to the ferromagnetic and antiferromagnetic cases, respectively. As the antiferromagnetic case is lower Table 1. Binding Energies ΔE [kcal/mol], Equilibrium Distances d0 [Å], and Minimum-Energy Rotation Angle a0 [deg] (at 3.3 Å) for All Three of the Systems Studied Herea CASSCF(2,2)-MP2

vdW-DF

system

ΔE

d0

a0

ΔE

d0

a0

C-Phy2 N-Phy2 B-Phy2

30.9 15.5 11.8

2.94 3.20 3.42

60° 30° 30°

14.6 11.4 9.4

3.31 3.56 3.74

60° 30° 43.5°

a

Note that the numbers quoted for CASSCF(2,2)-MP2 are those using the 6-31++G** basis. The reference state for all rotational calculations is the fully eclipsed dimer conformation at 0°. 3645

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649

The Journal of Physical Chemistry A

Article

MP2 energy minus the Hartree−Fock energy). Interestingly, electrostatics plus exchange exhibits a local energy maximum at 30° and a local minimum at 60°, opposite the trend in the total energy. The MP2 correlation energy plays little role in the rotational total energy curve but weakly favors a local minimum at 60°. The kinetic energy dominates the rotational energy curve, forcing local maxima at both 0° and 60°. Thus, it is the kinetic energy of the electrons that is responsible for the total energy barrier at 60° in Figure 6. Further analysis shows that the kinetic barrier originates in the antibonding HOMO. This orbital is shown in Figure 8 at several rotational angles. At a

multicenter bond suppressed (Figure 3). However, this is not the case. The rotational curve for N-Phy2 is shown in Figure 6. As can be seen in the figure, the rotational curve is qualitatively

Figure 6. Rigid rotational energy scan for the N-Phy2 dimer as calculated by MP2 and vdW-DF at a separation of 3.3 Å. For a detailed analysis of the barrier at 60°, see Figure 7.

different from the carbon case, which showed an energy minimum at 60°. The N-Phy2 rotation curve has a minimum at 30° and a local maximum at 60°. In fact, the presence of such a strong rotational preference in the N-Phy2 system is somewhat surprising given its lack of a 2-electron/multicenter bond. When the 2-electron/multicenter bond is suppressed in CPhy2, there is only a weak preference for an angle near 40° (Figure 3), but nothing nearly as strong as the almost 8 kcal/ mol rotational barrier at 60° in N-Phy2. Although it would be satisfying to explain the rotational barrier at 60° evident in Figure 6 in terms of a simple picture such as the interaction of frontier orbitals, such a simplistic explanation fails in this case. The full occupation of both the dimer bonding and antibonding orbitals yields an effectively nonbonding combination that does not explain the barrier. The 2-electron/multicenter bond present in C-Phy2 is disrupted here because the two highest-energy electrons in N-Phy2 occupy the antibonding orbital corresponding to that bond. Moreover, although both the electron−electron repulsion and nucleus−nucleus repulsion favor a 60° (staggered) rotation, the electron−nucleus interaction strongly favors an eclipsed conformation, making the overall electrostatic energy rise with torsion away from 0°. The origins of the barrier at 60° is most readily uncovered by splitting the energy into physically meaningful components. Figure 7 shows the total energy split into contributions from electrostatics plus exchange, electron kinetic energy, and the MP2 correlation energy (defined as the

Figure 8. Antibonding highest-occupied molecular orbital (calculated within spin restricted Hartree−Fock) for the N-Phy2 system at (a) 0°, (b) 30°, and (c) 60°.

relative rotation of 0°, the orbital exhibits a nodal plane between the two monomers. The lobes of the orbital are tightly confined on the α-carbons. As the dimer is rotated to 30°, the nodal plane between the monomers vanishes and the orbital lobes are allowed to spread out and connect the monomers. At a relative rotation of 60°, the nodal plane is back, once again tightly confining the lobes on the α carbons. This confinement of the HOMO is the origin of the kinetic energy barrier which manifests itself in the total-energy barrier evident in Figure 6. The binding energy for N-Phy2 as a function of separation (at a rotation angle of 30°) was also calculated and is shown in Figure 4. The curve looks similar to the C-Phy2 case although the binding is not quite as strong. For the DFT calculation, the optimum binding energy is 11.4 kcal/mol at a distance of 3.56 Å. For MP2, the optimum binding energy is 15.5 kcal/mol at a distance of 3.20 Å. Again, these results are consistent with the overbinding of MP2 and the overestimation of the binding distance for vdW-DF. Note that Mota et al. obtained a binding energy of 16.0 kcal/mol for this system at the MP2 level using the 6-31+G* basis set.48 B-Phy2. Another way to remove the radical nature of C-Phy is to replace the central carbon with a boron atom to form BPhy. Boron has one fewer electron than carbon so the SOMO of C-Phy becomes unoccupied in B-Phy. When two B-Phy monomers are brought together, the HOMO−HOMO overlap forms a dimer bonding and a dimer antibonding orbital, each doubly occupied. These are not the same orbitals that are formed from HOMO−HOMO overlap in N-Phy2 because there are four fewer electrons in B-Phy2. The rotational energy curve for B-Phy2 is shown in Figure 9. The curves for MP2 and vdW-DF look nearly identical up to about 30°, at which point they qualitatively diverge. The vdW-DF curve shows a rotational barrier of several kcal/mol at 60°, similar to but smaller than that in N-Phy2, whereas the MP2 curve shows a relatively flat energy curve, similar to that of the triplet C-Phy2 with the 2-electron/multicenter bond suppressed. The most obvious origin of the difference would be the treatment of correlation within the two sets of calculations (vdW-DF vs MP2). However, correlation effects are minimal in the

Figure 7. Breakdown of the computed N-Phy2 total MP2 energy from Figure 6 into a kinetic piece, an electrostatic plus exchange piece, and a (MP2) correlation piece. Both electrostatics plus exchange and correlation favor a local minimum at 60° but the kinetic energy dominates with a local maximum at that angle. 3646

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649

The Journal of Physical Chemistry A

Article

predicting a binding energy of 11.8 kcal/mol at a separation of 3.42 Å.

■

CONCLUSION We have demonstrated the applicability of a simple approximation for the spin-polarization of the vdW-DF in an extremely delicate chemical system. The results calculated here for C-Phy2 agree well both with experiment and, to a large extent, with CASSCF(2,2)-MP2 calculations, the only exception being that the binding distance is predicted to be slightly too large, but this is typical of the vdW-DF. These results pave the way for the use of the vdW-DF and its variants on molecular crystals involving C-Phy with an aim toward a fuller understanding of the crystal packing-conductivity relationship. It is hoped that this functional will be used as an efficient means to predict useful substitutions to effect greater conductivity in both molecular crystals in the family of radical conductors such as those based on C-phy, and others in the much larger class of charge-transfer salts containing π-stacks of cationic or anionic radicals. Further, an unexpected rotational barrier in the nitrogen-substituted N-Phy2 system was discovered. This barrier is attributed to significant changes in the electron kinetic energy upon rotation of the dimer. Finally, it was found that the rotation curve of the closed-shell, boron substituted B-Phy2 system exhibits a spurious rotational barrier at 60°. This effect was found to emanate from the local Slater exchange functional and is not mitigated by addition of either a gradient correction or a self-interaction correction.

Figure 9. Rigid rotational energy curve for B-Phy2 as calculated with both MP2 and vdW-DF at 3.3 Å separation.

rotational energy curves. What is more important is the treatment of exchange. To examine the effects of exchange in more detail, the BPhy2 rotational energy curve was calculated using exact exchange, Slater local exchange, and Slater exchange coupled with a PBE gradient correction. These calculations were performed in the NWchem44 package using a cc-pVDZ basis set and the densest available integration grid. To keep correlation effects from clouding the picture, all correlation was turned off for these calculations. The results are shown in Figure 10. As the figure shows, exact exchange shows a minimal

■

ASSOCIATED CONTENT

S Supporting Information *

Figures like that of Figure 8 for the B-Phy2 and C-Phy2 systems. Curve of the difference between the nonlocal piece of the exchange−correlation energy when calculated using our spinpolarization approximation and a spin-unpolarized calculation for the C-Phy2 binding curve. This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 10. Rotational dependence of the energy (at 3.3 Å separation) of B-Phy2 when exact exchange, local exchange, and gradient-corrected semilocal exchange are used. Local exchange exhibits a spurious energy barrier to rotation at 60°, which is not removed by semilocal gradient correction. No correlation functional was used for the calculations in this figure.

■

rotational barrier at 60°, whereas local exchange gives an almost 4 kcal/mol barrier. Addition of a PBE gradient correction to the exchange does not solve the problem. Thus, the difference evident in Figure 9 does not stem from MP2 vs vdW-DF correlation but from the underlying Hartree−Fock exact exchange vs approximate exchange. In other words, the approximate nature of Slater exchange leads to a slightly wrong density around 60°, compared to the density resulting from exact exchange. To better understand this finding, we analyzed the density dif ference between the dimer calculated with exact exchange and Slater exchange at 60° and find that Slater exchange wrongly tends to move charge to the α carbons, which overlap at 60°, resulting in the observed spurious barrier through steric hindrance. Solutions within DFT thus have to go beyond the standard local and gradient-corrected functionals, such as hybrid functionals. It should be noted that we see no similar difficulty in the exchange part of the functional in the NPhy2 and C-Phy2 cases because there the vdW-DF and MP2 curves qualitatively agree. The dissociation curves for B-Phy2 as calculated with vdWDF and MP2 are shown in Figure 4. Again, the agreement is good between the two approaches with vdW-DF predicting a binding energy of 9.4 kcal/mol at a distance of 3.74 Å and MP2

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported in part by the National Science Foundation under Grant No. CHE-1006702 and PHY1125915.

■

REFERENCES

(1) Goto, K.; Kubo, T.; Yamamoto, K.; Nakasuji, K.; Sato, K.; Shiomi, D.; Takui, T.; Kubota, M.; Kobayashi, T.; Yakusi, K.; et al. A Stable Neutral Hydrocarbon Radical: Synthesis, Crystal Structure, and PHysical Properties of 2,5,8-Tri-tert-butyl-phenalenyl. J. Am. Chem. Soc. 1999, 121, 1619−1620. (2) Itkis, M. E.; Chi, X.; Cordes, A. W.; Haddon, R. C. MagnetoOpto-Electronic Bistability in a Phenalenyl-Based Neutral Radical. Science 2002, 296, 1443−1445. (3) Tagami, K.; Wang, L.; Tsukada, M. Interface Sensitivity in Quantum Transport Through Single Molecules. Nano Lett. 2004, 4, 209−212.

3647

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649

The Journal of Physical Chemistry A

Article

(4) Dutta, P.; Maiti, S. K.; Karmakar, S. Multi-Termial electron Transport Through Single Phenalenyl Molecule: A Theoretical Study. Org. Electron. 2010, 11, 1120−1128. (5) Fan, Z.-Q.; Zhang, Z.-H.; Ming, Q.; Tang, G.-P.; Chen, K.-Q. First-Principles Study of Repeated Current Switching in a Bimolecular Device. Comput. Mater. Sci. 2012, 53, 294−297. (6) Tian, Y. H.; Kertesz, M. Is There a Lower Limit to the CC Bonding Distances in Neutral Radical π-Dimers? The Case of Phenalenyl Derivatives. J. Am. Chem. Soc. 2010, 132, 10648−10649. (7) Chi, X.; Itkis, M. E.; Patrick, B. O.; Barclay, T. M.; Reed, R. W.; Oakley, R. T.; Cordes, A. W.; Haddon, R. C. The First PhenalenylBased Neutral Radical Molecular Conductor. J. Am. Chem. Soc. 1999, 121, 10395−10402. (8) Mandal, S. K.; Samanta, S.; Itkis, M. E.; Jensen, D. W.; Reed, R. W.; Oakley, R. T.; Tham, F. S.; Donnadieu, B.; Haddon, R. C. Resonatng Valence Bond Ground State in Oxygen-Functionalized Phenalenyl-Based Neutral Radical Molecular Conductors. J. Am. Chem. Soc. 2006, 128, 1982−1994. (9) Sarkar, A.; Itkis, M. E.; Tham, F. S.; Haddon, R. C. Synthesis, Structure, and Physical Properties of a Partial π-Stacked PhenalenylBased Neutral Radical Molecular Conductor. Chemistry 2011, 17, 11576−11584. (10) Beer, L.; Mandal, S. K.; Reed, R. W.; Oakley, R. T.; Tham, F. S.; Donnadieu, B.; Haddon, R. C. The First Electronically Stabilized Phenalenyl Radical: Effects of Substituents on Solution Chemistry and Solid-State Structure. Cryst. Growth Des. 2007, 7, 802−809. (11) Wan, X.; Lv, X.; He, G.; Yu, A.; Chen, Y. Synthesis of Neutral Stable Polyradicals and Their Application on Phtovoltaic Devices. Eur. Polym. J. 2011, 47, 1018−1030. (12) Pal, S. K.; Itkis, M. E.; Tham, F. S.; Reed, R. W.; Oakley, R. T.; Donnadieu, B.; Haddon, R. C. Phenalenyl-Based Neutral Radical Molecular Conductors: Substituent Effects on Solid-State Structures and Properties. J. Am. Chem. Soc. 2007, 129, 7163−7174. (13) Liao, P.; Itkis, M. E.; Oakley, R. T.; Tham, F. S.; Haddon, R. C. Light-Mediated C-C σ-Bond Driven Crystallization of a Phenalenyl Radical Dimer. J. Am. Chem. Soc. 2004, 126, 14297−142302. (14) Pal, S. K.; Itkis, M. E.; Tham, F. S.; Reed, R. W.; Oakley, R. T.; Haddon, R. C. Trisphenalenyl-Based Neutral Radical Molecular Conductor. J. Am. Chem. Soc. 2008, 130, 3942−3951. (15) Huang, J.; Kertesz, M. Spin Crossover of Spiro-Biphenalenyl Neutral Radical Molecular Conductors. J. Am. Chem. Soc. 2003, 125, 13334−13335. (16) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D.; Lundqvist, B. van der Waals Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 22−25. (17) Thonhauser, T.; Cooper, V. R.; Li, S.; Puzder, A.; Hyldgaard, P.; Langreth, D. C. van der Waals Density Functional: Self-Consistent Potential and the Nature of the van der Waals Bond. Phys. Rev. B 2007, 76, 125112. (18) Langreth, D. C.; Lundqvist, B. I.; Chakarova-Käck, S. D.; Cooper, V. R.; Dion, M.; Hyldgaard, P.; Kelkkanen, A.; Kleis, J.; Kong, L.; Li, S.; et al. A Density Functional for Sparse Matter. J. Phys.: Condens. Matt. 2009, 21, 084203. (19) Thonhauser, T.; Puzder, A.; Langreth, D. C. Interaction Energies of Monosubstituted Benzene Dimers via Nonlocal Density Functional Theory. J. Chem. Phys. 2006, 124, 164106. (20) Li, S.; Cooper, V. R.; Thonhauser, T.; Puzder, A.; Langreth, D. C. A Density Functional Theory Study of the Benzene-Water Complex. J. Phys. Chem. A 2008, 112, 9031−9036. (21) Björk, J.; Hanke, F.; Palma, C.-A.; Samori, P.; Cecchini, M.; Persson, M. Adsorption of Aromatic and Anti-Aromatic Systems on Graphene through π-π Stacking. J. Phys. Chem. Lett. 2010, 1, 3407− 3412. (22) Mura, M.; Gulans, A.; Thonhauser, T.; Kantorovich, L. Role of van der Waals Interaction in Forming Molecule-Metal Junctions: Flat Organic Molecules on the Au(111) Surface. Phys. Chem. Chem. Phys. 2010, 12, 4759−4767.

(23) Chen, D.-L.; Al-Saidi, W.; Johnson, J. Nobel Gases on Metal Surfaces: INsights on Adsorption Site Preference. Phys. Rev. B 2011, 84, 2−5. (24) Cooper, V. R.; Thonhauser, T.; Puzder, A.; Schröder, E.; Lundqvist, B. I.; Langreth, D. C. Stacking interactions and the Twist of DNA. J. Am. Chem. Soc. 2008, 130, 1304−1308. (25) Cooper, V. R.; Thonhauser, T.; Langreth, D. C. An Application of the van der Waals Density Functional: Hydrogen Bonding and Stacking Interactions Between Nucleobases. J. Chem. Phys. 2008, 128, 204102. (26) Kolb, B.; Thonhauser, T. van der Waals Density Functional Study of Energetic, Structural, and Vibrational Properties of Small Water Clusters and Ice Ih. Phys. Rev. B 2011, 84, 045116. (27) Lin, Y.; Ma, H.; Matthews, C. W.; Kolb, B.; Sinogeikin, S.; Thonhauser, T.; Mao, W. L. Experimental and Theoretical Studies on a High Pressure Monoclinic Phase of Ammonia Borane. J. Phys. Chem. C 2012, 116, 2172−2178. (28) Bil, A.; Kolb, B.; Atkinson, R.; Pettifor, D.; Thonhauser, T.; Kolmogorov, A. van der Waals Interactions in the Ground State of Mg(BH4)2 from Density Functional Theory. Phys. Rev. B 2011, 83, 224103. (29) Klimeš, J.; Bowler, D.; Michaelides, A. van der Waals Density Functionals Applied to Solids. Phys. Rev. B 2011, 83, 195131. (30) Londero, E.; Schröder, E. Role of van der Waals Bonding in the Layered Oxide V2O5: First-Principles Density-Functional Calculations. Phys. Rev. B 2010, 82, 1−8. (31) Hamada, I.; Lee, K.; Morikawa, Y. Interation of Water with a Metal Surface: Importance of van der Waals Forces. Phys. Rev. B 2010, 81, 115452. (32) Román-Pérez, G.; Soler, J. M. Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes. Phys. Rev. Lett. 2009, 103, 096102. (33) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate ab Initio Parameterization of Density functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (34) Zhao, Y.; Truhlar, D. G. Hybrid Meta Density Functional Theory Methods for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions: The MPW1B95 and MPWB1K Models and Comparative Assessments for Hydrogen Bonding and van der Waals Interactions. J. Phys. Chem. A 2004, 108, 6908−6918. (35) Zhao, Y.; Truhlar, D. G. A New Local Density Functional for Main-Group Thermochemistry, Transition Metal Bonding, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Phys. 2006, 125, 194101. (36) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid ExchangeCorrelation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51−57. (37) Vydrov, O. A.; Heyd, J.; Krukau, A. V.; Scuseria, G. E. Importance of Short-Range Versus Long-Range Hartree-Fock Exchange for the Performance of Hybrid Density Functionals. J. Chem. Phys. 2006, 125, 074106−074114. (38) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (39) Kolb, B.; Thonhauser, T. Molecular Biology at the Quantum Level: Can Modern DFT Forge the Path? Nano LIFE 2012, 02, 1230006. (40) 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.; et al. Gaussian 03, Revision D.02; Gaussian, Inc.: Wallingford, CT, 2004. (41) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys.: Condens. Matter 2009, 21, 395502 (19pp). (42) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals Density Functional Made Simple. Phys. Rev. Lett. 2009, 103, 063004. 3648

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649

The Journal of Physical Chemistry A

Article

(43) Vydrov, O. A.; Voorhis, T. V. Nonlocal van der Waals Density Functional: The Simpler the Better. J. Chem. Phys. 2010, 133, 244103. (44) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; et al. NWChem: A Comprehensive and Scalable Open-Source Solution for Large Scale Molecular Simulations. Comput. Phys. Commun. 2010, 181, 1477−1489. (45) Slater, J. C.; Johnson, K. H. Self-Consistent-Field χα Cluster Method for Polyatomic Molecules and Solids. Phys. Rev. B 1972, 5, 844−853. (46) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (47) Suzuki, S.; Morita, Y.; Fukui, K.; Sato, K.; Shiomi, D.; Takui, T.; Nakasuji, K. Aromaticity on the Pancake-Bonded Dimer of Neutral Phenalenyl Radical as Studied by MS and NMR Spectroscopies and NICS Analysis. J. Am. Chem. Soc. 2006, 128, 2530−2531. (48) Mota, F.; Miller, J. S.; Novoa, J. J. Comparative Analysis of the Multicenter, Long Bond in [TCNE].− and Phenalenyl Radical Dimers: A Unified Description of Multicenter, Long Bonds. J. Am. Chem. Soc. 2009, 131, 7699−7707. (49) Small, D.; Zaitsev, V.; Jung, Y.; Rosokha, S. V.; Head-Gordon, M.; Kochi, J. K. Intermolecular π-to-π Bonding Between Stacked Aromatic Dyads. Experimental and Theoretical Binding Energies and Near-IR Optical Transitions for Phenalenyl Radical/Radical versus Radical/Cation Dimerizations. J. Am. Chem. Soc. 2004, 126, 13850− 13858. (50) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (51) Tian, Y.-H.; Huang, J. S.; Kertesz, M. Fluxional σ-Bonds of 2,5,8-Tri-tert-butyl-1,3-diazaphenalenyl Dimers: Stepwise [3,3], [5,5] and [7,7] Sigmatropic Rearrangements via n-Dimer Intermediates. Phys. Chem. Chem. Phys. 2010, 12, 5084−5093.

3649

dx.doi.org/10.1021/jp3095424 | J. Phys. Chem. A 2013, 117, 3642−3649