Intermolecular interactions in solid benzene

THE JOURNAL OF CHEMICAL PHYSICS 124, 044514 共2006兲 Intermolecular interactions in solid benzene G. J. Kearleya兲 Department of Radiation, Radionuclides and Reactors, Faculty of Applied Sciences, Delft University of Technology, Mekelweg 15, 2629 JB Delft, The Netherlands M. R. Johnson Institut Laue-Langevin (ILL), BP 156, 38042 Grenoble, Cedex 9, France J. Tomkinson CCLRC, The ISIS Facility, Rutherford Appleton Laboratory, Chilton, OX fordshire OX11 0QX, United Kingdom 共Received 25 August 2005; accepted 8 November 2005; published online 30 January 2006兲 The lattice dynamics and molecular vibrations of benzene and deuterated benzene crystals are calculated from force constants derived from density-functional theory 共DFT兲 calculations and compared with measured inelastic neutron-scattering spectra. A very small change 共0.5%兲 in lattice parameter is required to obtain real lattice-mode frequencies across the Brillouin zone. There is a strong coupling between wagging and breathing modes away from the zone center. This coupling and sensitivity to cell size arises from two basic interactions. Firstly, comparatively strong interactions that hold the benzene molecules together in layers. These include an intermolecular interaction in which H atoms of one molecule link to the center of the aromatic ring of a neighboring molecule. The layers are held to each other by weaker interactions, which also have components that hold molecules together within a layer. Small changes in the lattice parameters change this second type of interaction and account for the changes to the lattice dynamics. The calculations also reveal a small auxetic effect in that elongation of the crystal along the b axis leads to an increase in internal pressure in the ac plane, that is, elongation in the b direction induces expansion in the a and c directions. © 2006 American Institute of Physics. 关DOI: 10.1063/1.2145926兴 INTRODUCTION Benzene is generally regarded as the prototypical example of an aromatic molecule and in the solid state it provides the simplest real system in which interactions between aromatic molecules can be studied. The aim of the present work is to understand the intermolecular interactions at a microscopic level that lead to molecular packing, lattice dynamics, phase behavior, and ultimately the possibility of new bulk properties in aromatic systems. Over the years there have been a number of studies that are pertinent to this work: crystal structure,1 phases,2,3 and vibrational spectra,4–9 from both the experimental and theoretical standpoints. The present work is mainly concerned with the calculation of intermolecular interactions using density-functional theory 共DFT兲 methods, but in order to connect with experiment we compare our calculations not only with existing crystallographic data and optical vibrational spectroscopies but also with new inelastic neutron-scattering 共INS兲 spectroscopy data. In this way we establish an almost parameter-free model that is capable of reproducing the static and dynamic structure factors. The use of DFT methods for periodic systems for the determination of molecular vibrations and zonecenter 共⌫ point兲 lattice modes, and comparison of these with INS spectra, has become common place in recent years.10,11 Here these methods are extended to full lattice dynamics calculations, taking into account the whole Brillouin zone. a兲 Electronic mail: g.j.kearley@tnw.tudel 0021-9606/2006/124共4兲/044514/9/$23.00 Increasingly, these calculations are being used to simulate not only the coherent INS spectra from single crystals but also the incoherent INS spectra of powdered samples.12–14 The latter is a far more straightforward experimental technique and alleviates the need for large single crystals of deuterated materials in the study of lattice dynamics. Having established a model that reproduces the experimental data we shall exploit it to discern three major intermolecular interactions that hold the crystal together. One of these interactions is between the H atom of one molecule with the center of the aromatic ring of a neighboring molecule. Clearly, this is unique to aromatic systems and it is important to establish the relative strength of this interaction by comparing simulation and experiment, mainly lattice dynamics in this case, and by investigating the effects of uniaxial and isotropic pressures in the simulation. EXPERIMENT Benzene and deuterated benzene were obtained from The Aldrich Chemical Company and used without further purification. Samples were loaded in aluminum sample containers and cooled to 15 K using a standard cryostat. Data were collected using the now defunct TFXA spectrometer 共replaced with TOSCA兲 共Ref. 15兲 at the ISIS pulsed neutron facility in the UK. Raw data were transformed into S共Q␻兲 using standard algorithms. 124, 044514-1 © 2006 American Institute of Physics Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-2 Kearley, Johnson, and Tomkinson J. Chem. Phys. 124, 044514 共2006兲 FIG. 2. Dispersion curves in the low-frequency region for C6H6 using the experimentally determined Ref. 1 unit-cell parameters: a = 7.355 Å, b = 9.371 Å, and c = 6.700 Å 共cell I兲. The negative values are used to indicate the magnitude of the imaginary frequencies. RESULTS AND DISCUSSION Energy minimization and unit-cell optimization FIG. 1. Schematic illustration of the crystal structure of C6D6 from Ref. 1 showing the relative orientation of molecules within the layers. The long b axis is vertical. CRYSTALLOGRAPHIC INFORMATION In a DFT calculation of solid-state structure and dynamics of a molecular crystal, the only input is the measured crystal structure. For benzene, the structure has been determined at 4.2 K to be orthorhombic in the space-group Pbca, the unit cell and its contents being illustrated in Fig. 1. The cell parameters are a = 7.355, b = 9.371, and c = 6.700 Å.1 COMPUTATIONAL METHODS Energy calculations and structural optimization were made using VASP4.5,16 using the PBE exchange-correlation functional and PAW pseudopotentials with an energy cutoff of 600 eV. A single crystallographic unit cell was used for all calculations, with the reciprocal lattice being sampled using eight k points 共=关2 , 2 , 2兴兲. Single-point energy calculations were made for a series of structures in which the crystallographically distinct atoms were displaced by 0.03 Å in positive and negative directions along the x, y, and z directions. These calculations gave the Hellmann-Feynmann 共HF兲 forces acting on each atom and were used as input for the lattice-dynamics program PHONON4.2.4.17 Nonzero force constants were determined using the single unit cell, and it was found that all of these decayed by more than three orders of magnitude in going from the cell center to the nearest cell boundary. Phonon was used to calculate the eigenfrequencies, dispersion curves, and simulated inelastic neutronscattering spectra S共Q␻兲. The same HF forces were used for both C6H6 and C6D6, but the appropriate atomic masses and scattering cross sections were used in the lattice-dynamics calculations and INS calculations. The first step in vibrational analysis is the optimization of the crystal structure so that the total energy is a minimum and the forces acting on the atoms are zero. However, optimization of the unit-cell parameters of weakly bound molecular crystals using DFT is not straightforward because long-range attractive interactions due to mutually induced dipoles cannot, in principle, be built into a theory based on one-electron density such as DFT using local-density approximation 共LDA or GGA兲 exchange-correlation functionals. The fact that the dispersive interactions extend over the spatial range from ⬃3 to ⬃8 Å gives rise to a smoothly varying energy variation within the cell that can be considered in terms of a mean field. This reasoning underlies the correction applied here in which the unit-cell parameters are constrained to experimental 共or other兲 values in order to prevent unphysical cell expansion. For the present type of work, this is the only practical approach to the problem. Where dispersion has been calculated in other cases8,9 the approach has been found to work reasonably well. A more direct test of this correction is the calculation of weak rotational potentials for methyl groups, which depend significantly on van der Waals 共vdW兲 interactions and are obtained with a precision of ⬃90%, see, for example, Ref. 18, which does not include nonlocal, longrange correlation effects such as dispersive interactions. In the present work this is crucial since we will show that rather small changes to the unit-cell parameters have important effects on the lattice dynamics. At worst, this can be conceived as three adjustable parameters: the pressure along each of the crystallographic directions. This means that the calculated variation of properties as a function of pressure will be incorrect by a constant factor that will be rather close to unity, but that the trends will be correct. The starting model was taken from the most recent crystal structure determination 共4.2 K兲,1 the atomic positions being relaxed, but with the unit-cell parameters being held constant. A lattice dynamics calculation using the optimized structure resulted in dispersion curves illustrated in Fig. 2. It Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-3 Interactions in solid benzene FIG. 3. Dispersion curves in the low-frequency region for C6H6 using the bigger unit-cell parameters: a = 7.397 Å, b = 9.422 Å, and c = 6.737 Å 共cell II兲. Notice the absence of negative 共imaginary兲 values compared with Fig. 2. is immediately clear from this figure that the acoustic modes and one of the optic modes become imaginary around Y 共0 , 21 , 0兲, S共 21 , 21 , 0兲, and T共0 , 21 , 21 兲. While this may at first seem to be consistent with the proposal that this structure of benzene is only entropically stable,2 as we will show later, the lattice-mode INS spectrum calculated with this structure is in poor agreement with that measured. As alluded to above, there is a small uncertainty in the unit-cell parameters due to the shortcoming of the current DFT method, and consequently we investigated the effects of slight scaling of the unit-cell parameters. A number of calculations were made using larger and smaller unit cells and it transpired that reducing the unit-cell parameters 共0.5%2.0%兲 had little effect on the dispersion curves, while increasing the cell parameters by as little as 0.5% 共⬃0.04 Å兲 removed all imaginary frequencies with the exception of some very small values for the acoustic modes at the ⌫ point, k = 0. The dispersion curves for the low-energy region of the smaller unit cell are illustrated in Fig. 3, and the observed and calculated INS spectra are shown in Fig. 4. FIG. 4. Observed 共upper兲 and calculated 共lower兲 INS profiles for C6H6 using the model cell II. Calculation includes multiphonon modes up to five and convolution with an analytical instrumental resolution function. The inset shows the lattice-mode region: observed 共upper兲, cell II 共middle兲, and cell I 共lower, dashed兲. J. Chem. Phys. 124, 044514 共2006兲 FIG. 5. Observed 共upper兲 and calculated 共lower兲 INS profiles for deuterated benzene, C6D6, using the model cell II. Calculation includes multiphonon modes up to five and convolution with an analytical instrumental resolution function. Molecular and lattice vibrations We will denote the crystallographically determined unit cell as I and the 0.5% enlarged cell as II. A comparison of Fig. 2 and 3 reveals considerably different lattice dynamics for such a small change in unit-cell size, without change of symmetry. The inset in Fig. 4 compares the observed and calculated INS spectra in the lattice-mode region for cells I and II, which clearly reveals lost spectral density of the acoustic modes in the experimental cell, with an almost complete absence of intensity in the region around 40 cm−1. In addition the cell-I calculation does not show a Debye-type spectrum in the limit ␻ − ⬎ 0. While the agreement between the observed spectrum and the spectrum from cell II is not perfect in this region, it is a vast improvement for such a small change in lattice parameters. Agreement between the observed and calculated spectra of the internal modes is rather good 共Fig. 4兲 and is similar for either unit cell because these modes are much less sensitive to weak intermolecular interactions. The observed and calculated spectra for C6D6 are shown in Fig. 5, and again it was found that cell II was required to avoid imaginary frequencies and to give good agreement in the low-energy part of the spectrum. Formal assignments of all modes for both isotopomers are given in Table I. Assignments for the calculated frequencies in this table are based on the eigenvectors, and comparison with the experimental values is based on symmetry species where possible, otherwise, best match. Dispersion for the internal modes is generally less than about 15 cm−1, with notable exception of the 1000 cm−1 region. This spectral region, between 980 and 1010 cm−1, is rather complicated because 12 crystal modes exist, arising from the three molecular modes: ring breathing ␯1, H wagging ␯5, and in-plane ring deformation ␯12, these modes 共and the proximate ␯17兲 being illustrated in Fig. 6. Dispersion of these modes is illustrated in Figs. 7共a兲 and 7共b兲 and 8 for cells II and I, respectively. We will first consider the larger cell, II. A comparison of Figs. 7共a兲 and 7共b兲 shows that the ␯17 wagging modes, between 960 and 980 cm−1 at the zone center, are essentially H-atom displacements over the whole zone. This consistent behavior across the zone is found for Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-4 Kearley, Johnson, and Tomkinson J. Chem. Phys. 124, 044514 共2006兲 TABLE I. Vibrational frequencies of crystalline benzene and deuterated benzene at the ⌫ point 共k = 0兲. Observed frequencies are given in parentheses. Calculated assignments and their symmetry species are based on the atomic displacements. R and I denote Raman and infrared active modes, respectively. Agreement with the INS data is based on comparison of the observed and calculated spectral profiles 共Figs. 4 and 5兲. C6H6共cm−1兲 C6D6共cm−1兲 Assignment 共crystal兲 Molecular mode −1 −1 0 55 共57兲a 56 57 共57兲a 59 共61兲a 共53兲b 64 共70兲b 66 72 共79兲a 74 共53兲b 78 共79兲a 80 共100兲a 87 共90兲a 89 共92兲a 共94兲b 95 共92兲a 98 102 共85兲b 104 共94兲b 126 共128兲a 127 共128兲a 397 399 共401兲c 399 共404兲c 402 共405兲c 412 412 597 598 598 600 共600兲d 601 602 603 共606兲d 664 666 共685兲c 691 共674兲c 692 共690兲c 703 703 707 708 838 839 842 847 850 853 共852兲d 866 867 961 962 963 共974兲c 966 967 共983兲c 977 978 共977兲c −1 −1 0 50 54 53 53 57 62 64 66 71 71 73 79 81 86 86 95 93 95 101 115 117 345 347 共353兲c 348 351 361 361 共364兲c 568 570 571 571 572 571 573 576 489 490 共515兲c 507 共506兲c 508 共525兲c 599 600 609 611 653 653 656 660 662 665 671 672 784 784 785 共795兲c 788 786 788 792 792 共799兲c B2u共I兲 B3u共I兲 B1u共I兲 Ag共R兲 Au B1g共R兲 B3g共R兲 + B2u共I兲 B1u共I兲 Au Ag共R兲 B3u共I兲 B2g共R兲 B1g共R兲 B2g共R兲 Ag共R兲 + B3u共I兲 B3g共R兲 Au B2g共R兲 + B1u共I兲 B2u共I兲 B3g共R兲 B1g共R兲 B2u共I兲 + B3u共I兲 Au Au + B1u共I兲 B2u共I兲 B3u共I兲 B1u共I兲 B2g共R兲 Ag共R兲 B3g共R兲 Ag共R兲 + B1g共R兲 B1g共R兲 B2g共R兲 B3g共R兲 B3u共I兲 B1u共I兲 Au B2u共I兲 B3g共R兲 B1g共R兲 B2g共R兲 Ag共R兲 Ag共R兲 B1g共R兲 Ag共R兲 B2g共R兲 B3g共R兲 B2g共R兲 B3g共R兲 B1g共R兲 B3u共I兲 B2u共I兲 Au B1u共I兲 Au + B2u共I兲 B3u共I兲 B1u共I兲 Lattice modes ␯16E2u ␯6E2g ␯11A2u ␯4B2g ␯10E1g ␯17E2u Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-5 Interactions in solid benzene J. Chem. Phys. 124, 044514 共2006兲 TABLE I. 共Continued.兲 C6H6共cm−1兲 d 988 共978兲 989 共983兲d 992 992 共1010兲c 993 995 共991兲d 996 共997c,1006d兲 997 1003 1006 共1006兲d 1031 1034 1035 共1035兲c 1036 1039 共1038兲c 1040 1041 1132 1134 共1142兲c 1138 1142 1154 1159 1159 1160 1162 共1169兲d 1163 共1174兲d 1164 共1177兲d 1169 共1181兲d 1328 1330 1332 1334 1346 1348 1348 共1312兲c 1349 1459 1461 1463 共1470兲c 1463 1463 共1478兲c 1465 共1475兲c 1467 1468 1588 共1585兲d 1589 1590 1590 1592 3100 3102 3102 3106 共3044兲d 3107 共3048兲d 3107 3111 3111 3112 3113 C6D6共cm−1兲 Assignment 共crystal兲 Molecular mode 823 823 953 953 共970兲c 954 827 950 826 957 共967兲c 952 951 952 801 803 803 共810兲c 807 808 共814兲c 808 810 809 810 共823兲c 813 816 844 849 849 849 850 853 852 857 1036 1037 1039 1040 1324 1322 1325 共1286兲c 1324 1325 1325 1327 共1326兲c 1327 1344 1343 共1329兲c 1343 1344 1551 1552 1551 1552 1553 1555 1553 1555 2291 2292 2293 2298 2298 2299 2299 2302 2302 2303 2304 B3g共R兲␯5 B1g共R兲␯5 B3u共I兲␯12 B1u共I兲␯12 B2u共I兲␯12 Ag共R兲␯1 + B1g共R兲␯1 B2g共R兲␯1 + Au␯12 B3g共R兲␯1 B2g共R兲␯5 Ag共R兲␯5 B2u共I兲 B1u共I兲 + B3u共I兲 Au Au B3u共I兲 B1u共I兲 B2u共I兲 B3u共I兲 Au B1u共I兲 B2u共I兲 Ag共R兲 B1g共R兲 B2g共R兲 B3g共R兲 Ag共R兲 B2g共R兲 B3g共R兲 B1g共R兲 B3g共R兲 Ag共R兲 B1g共R兲 B2g共R兲 B3u共I兲 B1u共I兲 Au B2u共I兲 B1u共I兲 B2u共I兲 Au B3u共I兲 B2u共I兲 B1u共I兲 B3u共I兲 Au Ag共R兲 + B2g共R兲 + B3g共R兲 B1g共R兲 + B3g共R兲 Ag共R兲 B2g共R兲 B1g共R兲 B1u共I兲 + B3u共I兲 Au B2u共I兲 B2g共R兲 Ag共R兲 + B1g共R兲 B3g共R兲 B2g共R兲 Ag共R兲 B3g共R兲 B1g共R兲 ␯1A1g ␯5B2g ␯12B1u ␯18E1u ␯15B2u ␯9E2g ␯3A2g ␯14B2u ␯19E1u ␯8E2g ␯13B1u ␯7E2g Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-6 Kearley, Johnson, and Tomkinson J. Chem. Phys. 124, 044514 共2006兲 TABLE I. 共Continued.兲 C6H6共cm−1兲 C6D6共cm−1兲 Assignment 共crystal兲 Molecular mode 3117 3118 共3033兲c 3118 共3038兲c 3119 共3088兲c 3125 共3069兲c 3125 共3092兲c 3126 3128 共3061兲d 3130 3131 3131 2315 2315 共2267兲c 2315 共2278兲c 2316 2319 共2282兲c 2320 2320 2327 2328 2329 2329 B2u共I兲 Au B1u共I兲 B3u共I兲 B2u共I兲 Au B1u共I兲 + B3u共I兲 B3g共R兲 B1g共R兲 B2g共R兲 Ag共R兲 ␯20E1u ␯2A1g a Reference 4. Reference 5. c Reference 7. d Reference 6. b all modes except ␯1, ␯5, and ␯12, which mix strongly with each other at different wave vectors as is evident from Figs. 7共a兲 and 7共b兲 between 980 and 1010 cm−1. Away from the zone center, the displacements in these modes are clearly a mixture of the formal molecular modes: ring breathing 共symmetric and antisymmetric兲 with the out-of-plane wagging modes. The situation for the smaller cell, I, is markedly different 关Figs. 8共a兲 and 8共b兲兴. Not only is the dispersion in this region more pronounced but here the higher-frequency components of ␯17 also mix with ␯1, ␯5, and ␯12 away from the zone center. This mixing of out-of-plane H-wagging modes with inplane ring breathing modes away from the zone center suggests a significant intermolecular interaction between the H atoms of one molecule and the aromatic core of the neighboring molecule 共see below兲. The molecular center-of-mass displacement of these modes also varies across the zone due to the coupling of these internal modes with the lattice modes. FIG. 6. Schematic illustration of the atomic displacements in the molecular modes that arise in the 960–1010 cm−1 spectral region 共⌫ point兲. Intermolecular interactions The crystal structure of benzene illustrated in Fig. 1 is conveniently regarded as composed of layers of molecules stacked along the long b axis, the molecules in each ac layer are tilted by about 38° to b. Inspection of the crystal structure reveals the three types of interaction that are illustrated in Fig. 9. The interactions labeled A and B are between the layers, while the layers themselves are held together by in- FIG. 7. Mixing of formal molecular modes with each other and lattice modes across the zone. The behavior is not seen in other spectral regions. 共a兲 Dispersion curves for cell II in the regions of ␯17 共between 961 and 978 cm−1兲, ␯5 共⬃990 and ⬃1005 cm−1兲, ␯12 共⬃992 and 996 cm−1兲, and ␯1共⬃996 cm−1兲, these frequencies being at the ⌫ point, k = 0. The density of the lines 共gray scale兲 reflects the relative amplitude of the C-atom displacements. 共b兲 Same as 共a兲, but the intensity of the lines 共gray scale兲 represents the relative amplitude of the H-atom displacements. Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-7 Interactions in solid benzene J. Chem. Phys. 124, 044514 共2006兲 FIG. 10. Development of the electron density for interactions within the ac plane 共lower molecules兲 and between planes 共upper molecule兲, as the isosurface is decreased from 0.045e to 0.025e. Notice the large change in the H-aromatic overlap between 0.035e and 0.025e. FIG. 8. 共a兲 Dispersion curves for cell II in the regions of ␯17, ␯5, ␯12, and ␯1, to be compared with Fig. 7共a兲. The intensity of the lines 共gray scale兲 reflects the relative amplitude of the C-atom displacements. 共b兲 Same as 共a兲, but the intensity of the lines 共gray scale兲 represents the relative amplitude of the H-atom displacements. teraction D, and a more unusual type of interaction between the H atoms and the aromatic rings of neighboring molecules, labeled C. In this interaction, the distance from the H atom to each C atom of the neighboring ring is almost the same, varying from 3.015 to 3.061 Å, the distance from the H atom to the center of the ring 共labeled E兲 being only 2.701 Å. These interactions are examined in more detail in Fig. 10, where the electron-density isosurface has been calculated at three different values. At the isosurface 0.045e the beginning of interaction type C can be seen between the two lower molecules, but the H-aromatic link has not actually been made. Looking at lower electron density, 0.035, the interaction type D is established but there is no actual link between the upper and lower molecules, and the H-aromatic link is still missing. Finally, at the isosurface 0.025e, all links are established, but it is interesting to note that the greatest overlap is by a combination of interactions C and D, which effectively merge the densities of the two lower molecules. The picture that emerges from this electron-density figure is of molecules that are held together in layers 共lower molecules兲 by a quite strong interaction plus the H-aromatic interactions. The layers are linked to each other by weaker interactions. This picture is broadly consistent with the pattern of acoustic modes in Fig. 2 and 3, where the lowestenergy mode is translation along Y, while the corresponding modes for X and Z are generally about 50% higher. The anisotropy of the mean-square displacements, collected in Table II, is also in agreement with this picture. Considering the H atom involved in the H-aromatic interaction, C, displacement perpendicular to the approximate bond direction is noticeably larger than that in the x-z plane. Similarly, H-atom interactions involved in connecting neighboring planes of molecules, along y, show less overall displacement in this direction. The overall isotropic values are in reasonable agreement with those measured crystallographically.1 TABLE II. Calculated mean-square displacements for the crystallographically distinct atoms. The interaction types, a − d, are illustrated in Fig. 9. Figures in parentheses are for the H atom after subtraction of the C atom displacement. Atom interaction Atom interaction MSD-x共Å2兲 MSD-y共Å2兲 MSD-z共Å2兲 Isotropic 共Å2兲 FIG. 9. Illustration of the three major intermolecular interactions and their distances. The oval represents the center of the aromatic ring, the intermolecular distance to the nearest neighbor being E, 2.701 Å. The orientation of this fragment is similar to that of the unit cell in Fig. 1. Ha Hc Hd Ca Cc Cd 0.74共0.52兲 0.70共0.51兲 0.74共0.52兲 0.22 0.19 0.22 0.71共0.54兲 0.78共0.57兲 0.61共0.50兲 0.17 0.21 0.11 0.73共0.53兲 0.75共0.54兲 0.79共0.55兲 0.20 0.21 0.24 0.73共0.53兲 0.74共0.54兲 0.71共0.52兲 0.20 0.20 0.19 Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-8 Kearley, Johnson, and Tomkinson FIG. 11. Plot of calculated pressure in cell II along each direction as a function of the fractional change in the lattice parameter b 共for 1.000 b = 6.737 Å兲. The atomic positions were relaxed to the energy minimum for each value of c. Auxetic effect of pressure along b direction Intuitively we would expect the stronger interactions in the ac plane to be more sensitive to changes in the unit-cell size, but this is clearly not the case. In order to proceed we have investigated the structural consequences on progressively changing the unit-cell b parameter. Figure 11 shows the results of these calculations starting from the parameters of cell II, and it can be seen that within the range on the b parameters of cells I and II, there is a marked hardening of the cell along this direction. Perhaps more surprising is that in the range from 1.002 to 0.996, the pressure along b increases by 3.1 kbars, while the pressure along both a and c reduces by 0.7 kbars. In terms of cell parameters, shortening b leads to compression along a and c, giving a negative Poisson’s ratio or auxetic behavior. This behavior can be understood by considering Figs. 1 and 9. The structural differences at the points 0.99 and 1.01 in Fig. 11 correspond to changes in C and D 共Fig. 9兲, of only 0.002 Å this being consistent with the greater overlap of electron density shown on the right of Fig. 10. In contrast, distances A and B, between neighboring layers, decrease by 0.013 and 0.023 Å for compressing the b direction from 1.01 to 1.00, but then increase by 0.040 and 0.016 Å, respectively, when compressing the b direction further to 0.99. In order to establish the consistency of this effect we have also calculated the consequences of a 1% elongation of the cell along a. This leads to decreases in pressure along a and c of 1.9 and 1.2 kbars, respectively, but an increase of 0.9 kbars along b, this being entirely consistent with the results obtained above. Without going into the fine detail of the rather small molecular reorientations, the basic scheme is as follows. It is important to notice that interactions A and B also have components in the ac plane so that initially, as the lattice is compressed along b, interactions A and B increase pulling molecules in the ac plane together, reducing the pressure in this plane. Further compression 共below ⬃0.998兲 forces A and B interactions beyond their optimum, allowing relaxation in the xz plane. J. Chem. Phys. 124, 044514 共2006兲 This pattern of interactions accounts for the high sensitivity of the acoustic modes at points Y 共0 , 21 , 0兲, S共 21 , 21 , 0兲, and T共0 , 21 , 21 兲 in the dispersion curves shown in Figs. 2 and 3. The b parameter of cell I corresponds to 0.995 in Fig. 11, where A and B interactions are not optimal, and hence at some points in the Brillouin zone y displacements lead to an overall drop in the energy, and hence imaginary frequencies. When A and B interactions are near optimal, all frequencies are real 共apart from a very small error at k = 0 for the acoustic modes兲, as seen in Fig. 3. This change in A and B interactions not only accounts for the changes in the dispersion of the lattice modes between cells I and II 共Figs. 2 and 3兲 but also for the surprisingly large changes in the internal-mode dispersion of the H-wagging modes ␯5 and ␯17. Inspection of Figs. 1 and 6 reveals that all modes that wag the interlayer hydrogens will strongly modulate A and B 共Fig. 9兲. CONCLUSIONS DFT calculations are now sufficiently rapid and accurate to allow the vibrational dynamics of crystals, such as benzene, and to verify these calculations against experimental spectroscopies. The principle difficulty is the uncertainty in dispersion energy, but this can be overcome by using pressure to constrain the lattice to the experimental values, effectively preventing physically unrealistic lattice expansion. In the case of benzene, however, it is found that very small changes in unit-cell dimensions have a dramatic effect on the lattice dynamics and an expansion of only 0.5% above the experimentally determined values takes the cell from an unstable to a stable state. Although this leaves a small unknown scalar in the pressure, it is clear that the phonon dispersion and the dispersion of some of the internal modes depend crucially on small changes to the lattice parameters. The net interactions holding the molecules together in layers are stronger than those holding neighboring layers together. Because some interactions play both roles, forcing the layers together can increase the net interaction within the layers leading to a “contraction” of the layer. Changes in these interactions are entirely consistent with the sensitivity of the lattice modes and molecular vibrations to small changes in the unit-cell size. Constraining the unit-cell parameters to values close to those experimentally determined is effectively a correction of the DFT method to take account of dispersive interactions. This has the consequence of introducing an offset in the pressure of about 10 kbars, as seen in Fig. 11. This will have some small effect on the relative values at which the auxetic effect occurs, and it would be interesting to see if the predicted auxetic effect could be observed experimentally. Dispersion of internal modes is normally only important where strong hydrogen bonding interactions are involved. In the present case there is significant intermolecular coupling of wagging and breathing modes of the aromatic ring that can be seen as Davidov splitting at the zone center, but which couples strongly to optic and acoustic phonons away from the zone center causing extensive mixing. This interaction is considerably stronger than would be suggested by an analysis of the optical spectra alone. Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 044514-9 1 Interactions in solid benzene W. I. F. David, R. M. Ibberson, G. A. Jeffrey, and J. R. Ruble, Physica B 180 & 181, 597 共1992兲. 2 P. Raiteri, R. Martoňák, and M. Parrinello, Angew. Chem., Int. Ed. 44, 2 共2005兲. 3 R. Martoňák, A. Laio, M. Bernasconi, C. Ceriani, P. Raiteri, F. Zipoli, and M. Parrinello, Z. Kristallogr. 220, 489 共2005兲. 4 H. Bonadeo, M. P. Marzocchi, E. Castellucci, and S. Califano, J. Chem. Phys. 57, 4299 共1972兲. 5 G. Taddei, H. Bonadeo, M. P. Marzocchi, and S. Califano, J. Chem. Phys. 58, 966 共1973兲. 6 M. M. Thiéry and J. M. Léger, J. Chem. Phys. 89, 4255 共1988兲. 7 M. M. Thiéry, J. M. Besson, and J. L. Bribes, J. Chem. Phys. 96, 2633 共1992兲. 8 M. R. Johnson and H. P. Trommsdorff, Chem. Phys. Lett. 364, 34–38 共2002兲. 9 J. Stride, M. Adams, and M. R. Johnson, Chem. Phys. 317, 143 共2005兲. 10 L. van Eijck, M. R. Johnson, and G. J. Kearley, J. Phys. Chem. A 107, J. Chem. Phys. 124, 044514 共2006兲 8980 共2003兲. M. Montejo, A. Navarro, G. J. Kearley, J. Vazquez, and J. J. LopezGonzalez, J. Am. Chem. Soc. 126, 15087 共2004兲. 12 M. R. Johnson, K. Parlinski, I. Natkaniec, and B. Hudson, Chem. Phys. 53, 291 共2003兲. 13 P. Hermet, J.-L. Bantignies, A. Rahmani, J.-L. Sauvajol, M. R. Johnson, and F. Serein, J. Phys.: Condens. Matter 16, 7385 共2004兲. 14 P. Hermet, J.-L. Bantignies, A. Rahmani, J.-L. Sauvajol, and M. R. Johnson, J. Phys. Chem. A 109, 4202 共2005兲. 15 P. C. H. Mitchell, S. F. Parker, A. J. Ramirez-Cuesta, and J. Tomkinson, Vibrational Spectroscopy with Neutrons 共World Scientific, NJ, 2005兲. 16 G. Kresse and J. Furthmüller, software VASP, Vienna, 1999; G. Kresse, Phys. Rev. B 54, 169 共1996兲; Comput. Mater. Sci. 6, 15 共1996兲. 17 K. Parlinski, AIP Conf. Proc. 479, 121 共1998兲. 18 M. R. Johnson, M. Prager, H. Grimm, M. A. Neumann, G. J. Kearley, and C. C. Wilson, Chem. Phys. 244, 49 共1999兲. 11 Downloaded 10 Sep 2010 to 131.180.130.114. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions