Experimental and Theoretical Investigation of the First Overtone

Jul 13, 2011 - Aberdeen Proving Ground, Maryland 21005-5066, United States ... modeling study of the CАH stretch first overtone of bulk 1,3,5-trinitr...
0 downloads 0 Views 3MB Size
ARTICLE pubs.acs.org/JPCA

Experimental and Theoretical Investigation of the First Overtone Spectrum of 1,3,5-Trinitrotoluene Jerry B. Cabalo* and Rosario Sausa Aberdeen Proving Ground, Maryland 21005-5066, United States ABSTRACT: Vibrational overtone spectroscopy is a powerful tool for studying intramolecular and intermolecular interactions. We report on a combined experimental and modeling study of the CH stretch first overtone of bulk 1,3,5-trinitrotoluene (TNT) and TNT on fumed-silica powder. We recorded the overtone spectra by laser photoacoustic spectroscopy and compared them with those predicted with the harmonically coupled anharmonic oscillator model in the 56006600 cm1 region. The model systems included single molecules and hybrid quantum and molecular mechanical (QM:MM) clusters to account for the effects of intermolecular interactions on the observed spectra. We performed the hybrid QM:MM calculations at the HF/6-31+G(d,p), B3LYP/6-31+G(d, p), and MP2/6-31+G(d,p) levels of theory and with the universal force field (UFF) to account for van der Waals and electrostatic effects from surrounding molecules. Overtone spectra calculated from the MP2 level of theory, using a HF/3-21+G* calculation to assign molecular charges in the MM layer, and the MerzSinghKollman population analysis for assigning partial charge in the QM layer and determining the transition dipole moment agreed best with the experimental data.

’ INTRODUCTION Intermolecular interactions play an integral part in physical and chemical processes, such as adsorption, corrosion, and heterogeneous catalysis.16 An understanding of these processes has profound implications for the properties, functions, and reactivities of molecules. Vibrational overtones are uniquely sensitive to a molecule’s local environment and are potentially an effective tool in examining intermolecular and moleculesurface interactions because they involve higher vibrational energy levels than the ground state levels. Knowledge of how these interactions are reflected in the overtone spectra permits the refinement and validation of computational models that are in turn used to predict other important properties, such as molecular orientation, molecular motion, and binding energies. When the surrounding molecules perturb a molecule’s vibrational potential energy surface, the higher energy levels experience more change than the lower levels. The molecule’s average bond lengths increase with vibrational level, and thus, the molecule experiences greater mechanical interactions with the surroundings. This intermolecular interaction also affects the nuclear wave functions that in turn modify the transition dipole moment and the absorption intensities. As a result, the magnitude of the perturbation from the surroundings to the vibrational spectrum increases as higher vibrational levels are involved. The vibrational self consistent field (VSCF) method710 and the harmonically coupled anharmonic oscillators (HCAO) method1114 are two theoretical approaches that can be used to predict overtone spectra. The VSCF approach makes the approximation that the nuclear wave functions along each normal coordinate are r 2011 American Chemical Society

independent of the motions along the other coordinates. A consequence of this assumption is that the vibrational Hamiltonian can be broken into a one-dimensional vibrational Hamiltonian for each normal coordinate and that each single coordinate Hamiltonian can be solved with a self-consistent field approach. This is analogous to the HartreeFock approach in which the electronic wave function is treated as the product of single electron wave functions. HCAO is a much simpler, local mode model that applies only to bond stretches between a hydrogen atom and a heavier atom, such as carbon, and treats them as a diatomic Morse oscillator. HCAO also corrects for coupling between oscillators by treating that coupling as a harmonic, second-order interaction. Although VSCF is more general than HCAO, HCAO was implemented in this study because of its simplicity and applicability to the CH stretches that were the focus of this study. In this paper, we report our combined experimental and modeling study on the CH stretch first overtone of solid 1,3,5trinitrotoluene (TNT) because it is relevant to energetic materials, it can be an important environmental contaminant, and its crystal structure is well-known.15 We measured the overtone spectra of solid TNT and TNT on fumed silica by laser photoacoustic spectroscopy. The results were then compared with spectra calculated with the HCAO model at the HF/6-31+G(d,p), B3LYP/ 6-31+G(d,p), and MP2/6-31+G(d,p) levels of theory using single

Received: March 18, 2011 Revised: July 11, 2011 Published: July 13, 2011 9139

dx.doi.org/10.1021/jp2025895 | J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Molecular arrangement of the TNT models. Single molecule models are shown in part A. Positions in the cluster that were treated quantum mechanically (QM) and undergo HCAO analysis are shown as tube models, and molecules treated classically (MM) are shown as wireframes. Labeled atomic positions correspond to distances and bond angles shown in Table 4. Part B is a model of the TNT crystal unit cell that was optimized with the UFF and the partial charge assigned with the QEq algorithm. Part C is a large cluster model containing eight unit cells of the TNT crystal. The numbered positions indicate positions in the cluster that were modeled with HCAO.

molecule, 8 member clusters, and 64 member clusters. Our experimental data were used to test and refine our model. The objective is to obtain an understanding of how molecular interactions are reflected in the overtone spectra and to use this understanding for predicting overtone spectra from models of intermolecular interactions. The capability to predict changes in the overtone spectrum due to moleculebulk and moleculesurface interactions permits the refinement and validation of the computational models of intermolecular interactions. The validated models in turn allow other important properties, such as binding energies, to be predicted.

’ THEORETICAL METHODS Theoretical calculations were performed using Gaussian 2003 version B.04 on a Sunstation V as well as a LAN Parallel cluster utilizing 48 Intel processors and Gaussian 2009 version A.02.16 Overtone vibrational frequencies of TNT molecules were calculated for single molecules and for TNT molecules at various positions within 8 member and 64 member clusters of TNT molecules, as shown in Figure 1. The molecular cluster models were used to simulate the effect of intermolecular interactions within the TNT crystal as well as between TNT and a silica surface on the overtone spectrum. Two methods of assigning partial atomic charge were used on a large 64 member cluster model (shown in Figure 1C) that was arranged according to the crystal structure determined by X-ray crystallography.15 In the

first method, the geometry of the cluster was first optimized with the universal force field (UFF) before any subsequent calculations, in which the QEq algorithm was used to assign partial charges to all the atoms in the cluster. The cluster was then treated as a hybrid quantum mechanical/classical mechanical (QM:MM) model in which the molecule undergoing vibrational analysis was treated with ab initio/density functional methods (i.e., second order MoellerPlesset (MP2)), and the surrounding molecules were treated with UFF. The hybrid QM:MM calculations were implemented in the Gaussian code with the Our own N-layered Integrated molecular Orbital molecular Mechanics (ONIOM) method.16 The effect of intermolecular interactions was treated as mechanical and electrostatic perturbations from the MM layer to the QM layer. The second method of assigning partial charge used the HF/321+G* level of theory and the MerzSinghKollman (MK) population analysis. A second geometry optimization was performed with UFF and the ab initio assigned partial charges. Several positions were selected within the cluster to be completely surrounded by other molecules to avoid edge effects. For each selected position, the core molecule was treated with the MP2/6-31+G(d,p) level of theory. It was assumed that only interactions between nearest neighbors were important, and all molecules except the nearest neighbors to the core molecule were removed from the model for subsequent steps of the calculation. The core molecule’s geometry was then optimized while leaving the neighboring molecules frozen. The partial 9140

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

atomic charges were left on the molecules in the MM layer and were allowed to interact with the QM layer. The structure of the unit cell of a TNT crystal served as the initial geometry of the molecular cluster in Figure 1B. The molecules in Figure 1B,C shown as tube models indicate the positions of the molecules treated quantum mechanically. Only one molecule was treated quantum mechanically in a single calculation so that a series of calculations were performed to cover the positions shown. For the cluster in Figure 1B, position 1 was chosen so that the methyl hydrogens were directed into the cluster, and positions 2 and 3 were chosen so that the aryl hydrogens interacted with the interior of the cluster. The hybrid cluster model geometries were further optimized, with the molecule targeted for vibrational analysis treated quantum mechanically. Figure 1C shows a 64 member cluster model that had its partial charges assigned with an ab initio calculation and its geometry subsequently optimized. In the model of a single TNT molecule interacting with a silica surface (not pictured), the silica surface was modeled as amorphous, fully hydrated silica, and partial charge was assigned and the geometry optimized in the same manner as above. All the silicon atoms were in the tetrahedral geometry with no long-range order. Two possible approaches were utilized here for predicting overtone frequencies. The first approach utilized the local mode model,1720 of which HCAO21 is a variant and was used to simulate overtone transitions from the TNT molecule interacting with its surroundings. The HCAO model allowed calculation of both absorption frequencies and intensities of the overtone transitions. Because HCAO was shown in the literature to be successful at describing the absorptions of several molecules, including chemical agents,1113 only the details particular to this study of TNT are described below. The second approach applied anharmonic corrections to the frequencies calculated for the normal modes, which is well described in the literature and is available in the Gaussian code.2225 This approach was not generally used in this study because absorption intensities were not calculated and because of the computational expense of obtaining the third-order force constants. For TNT, there are a total of five CH oscillators: three oscillators on the methyl group and two oscillators that are part of the aromatic ring. The three methyl oscillators were treated as kinetically and potentially coupled, and the two aryl oscillators were treated as coupled potentially with the kinetic energy coupling neglected due to the physical separation of the oscillators. The Hamiltonian used for the methyl CH oscillators is shown in eq 1: ^ i j0æj j0æk ^ i j0æj j0æk  Hj0æ Hjυæ ¼ υωeðiÞ jυæi j0æj j0æk hc  ωe xeðiÞ ðυ2 þ υÞ jυæi j0æj j0æk   0 þ  þ  jij a a þ a a jυæi j0æj j0æk i j j i



i6¼ j



∑ jik i6¼ j 0

 þ  þ jυæi j0æj j0æk a i ak þ ak ai

ð1Þ

where υ is the excited state vibrational quantum number induced by optical absorption; |υæi, |0æj, and |0æk denote the Morse oscillator wave function of the ith, jth, and kth methyl CH oscillator with the ith oscillator absorbing the optical radiation, and the overall wave function is the product of the three individual wave functions; ωe(i), ωe(j), ωe(k) denote the Morse

oscillator harmonic constants, and ωexe(i), ωexe(j), and ωexe(k), the anharmonicity constants; ϕ0ij and ϕ0ik are the pairwise oscillator coupling parameters between the ith and jth, and ith and kth coupled oscillators; and a and a+ are the harmonic oscillator raising and lower operators. It is assumed they act on the Morse oscillator wave functions in the same way as the Harmonic wave functions. Insertion of the appropriate vibrational quantum numbers and Morse oscillator parameters gives the transition frequencies in cm1. For the coupled oscillators in the TNT molecule, the nuclear wave function for the coupled oscillators is written as products of the individual wave functions along each stretching coordinate. Because the aryl CH oscillators were physically isolated from the methyl CH oscillators, the Hamiltonian allows only for interaction between the aryl oscillators, and it is shown in eq 2, ^ i j0æj ^ i j0æj  Hj0æ Hjυæ ¼ υωeðiÞ jυæi j0æj hc   0 þ  þ  ωe xeðiÞ ðυ2 þ υÞ jυæi j0æj  jij a jυæi j0æj i aj þ aj ai



i6¼ j

ð2Þ where the variables are the same as in eq 1, except that the overall wave function is the product of only the two individual aryl oscillator Morse wave functions, and only the coupling parameter ϕ0ij between the two aryl oscillators is considered. Calculation of the Morse parameters and the coupling parameters were performed as described in the literature. The diagonal and offdiagonal elements of the Wilson F matrix26,27 and force constants were determined from the coefficients of a polynomial fit to the results of a potential energy surface scan along the CH stretching coordinates. Tabulated expressions26,27 were used for the elements of the Wilson G matrix. To extract force constants from the computational output, the geometry of the target molecule and its environment were energy-optimized with Gaussian 2009 at the HF/6-31+G(d,p), B3LYP/6-31+G(d,p), and MP2/6-31+G(d,p) levels of theory, then the potential energy surface (PES) was calculated as the bond lengths were varied (0.30 Å in steps of 0.05 Å around the equilibrium position. The potential energy scans were performed along two stretch coordinates simultaneously to capture any interaction between two coordinates. For the methyl CH oscillators, each oscillator was scanned pairwise so that interaction among all three CH oscillators was captured. The aryl CH coordinates and one of the aryl and methyl CH coordinates were scanned together to check for coupling. For determination of the cubic and square force constants, the potential energy as a function of motion along one coordinate was fitted to a fourth-order polynomial, and all other coordinates were held at the equilibrium position. The coefficients of the fit generated by Matlab R2006 were substituted into the equations described in Kjaergaard at al.,11,12,28,29 and the output of these equations yielded the Morse oscillator parameters. Similarly, the F matrix elements were determined from the fitting coefficients of a two-variable polynomial function for the potential energy surface with only second-order terms considered in the fitting function. To calculate predicted intensities, the matrix elements of the electric dipole moment operator, |Æψ1|μ|ψ0æ|,2 were calculated for each oscillator. The dipole moment function μ was approximated with a Taylor expansion. Terms up to the fifth order were 9141

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

used for the three coupled methyl CH oscillators (eq 3), and terms up to the fourth order were used for the aryl oscillators (eq 4), i i þ j þ k¼5  1 ∂i μ ∂j μ ∂k μ  μ q1 , q2 , q3 ¼ j k q1 þ r0, ðiÞ ∑ i i, j, k ¼ 1 i!j!k! ∂q 0 ∂q 0 ∂q 0  j  k q2 þ r0, ðjÞ q3 þ r0, ðkÞ ð3Þ  μ q4 , q5 ¼

i þ j¼4



i, j ¼ 1

i  j 1 ∂i μ ∂j μ  q þ r q þ r 4 5 0, ðiÞ 0, ðjÞ i!j! ∂qi ∂qj 0

0

ð4Þ where q1, q2, and q3 are the internuclear displacements along methyl CH bond stretches and q4 and q5 are the displacements for the aryl CH bonds. Fourth- and fifth-order terms were included because the higher-order terms can contribute significantly to the absorption intensity. The dipole moment derivatives were determined on the same grid points as the molecular potential energy in the quantum mechanical calculations by specifying an electron population analysis at each grid point. Several schemes for the population analysis were tried when the Mulliken, the MSK, and the Charges from Electrostatic Potentials (CHelp) analyses were used. The derivatives were found by fitting to polynomials as in eqs 3 and 4. The resulting dipole moment function and the Morse oscillator wave functions, as reported by Brabson30 but without the phase factors, were used in the calculation of the transition dipole moment |μe|2. Equation 5 was used for the aryl oscillators, and eq 6, for the methyl CH oscillators,  ð5Þ μe ¼ Æυj1 Æ0j2 μ q1 , q2 j0æ1 j0æ2  μe ¼ Æυj1 Æ0j2 Æ0j3 μ q1 , q2 , q3 j0æ1 j0æ2 j0æ3

ð6Þ

where the state with the quantum number υ indicates the oscillator that absorbs the IR radiation, and the υ = 0 states indicate ground state wave functions. The oscillator strength, f, of a given absorption line is then given by eq 7 and is combined with the Lorentzian line shape function to simulate the contribution of a given absorption line to the IR spectrum, 2 4πme ð7Þ f ¼ 2 ω0 μe 3e p where ω0 is the transition frequency, me is the mass of an electron, e is the elementary charge, and p is the reduced Planck’s constant. Calculated spectra were constructed by summing all the Lorentzian line functions at each absorption frequency.

’ EXPERIMENTAL METHOD To provide a basis for comparison with the synthetic spectra, overtone and fundamental vibrational spectra were collected using photoacoustic spectroscopy in a commercial MTEC 300 cell that was flushed with Argon gas to remove water vapor. Photoacoustic spectroscopy was used because it is suitable for powdered samples. Tunable near-IR laser light was obtained from the idler output of a 10 Hz Nd:YAG pumped Continuum

Sunlite EX OPO laser system in the range of 56006600 cm1. The pulse energy varied from 600 μJ/pulse at 6400 cm1, to 150 μJ/pulse at 5700 cm1. A Schott glass filter removed residual visible light from the near-IR beam. Laser pulse energy was measured using a Coherent Field-Max II meter with a J-10 energy sensor. The MTEC 300 had an onboard preamplifier with an amplification setting of 20. The output was directed to a Tektronix 5104B digital storage oscilloscope and a Stanford Research boxcar integrator, with a maximum gate setting of 3 μs and a delay of 10 ms. The output of the boxcar was digitized and recorded on a computer. Power spectra were recorded by scanning over the tuning range and monitoring the laser energy with a photodiode directed into the boxcar. To calibrate the power curve, the average pulse energy of the IR idler beam was measured with the pulse energy meter at 5 nm intervals of the visible output wavelength. All spectra were normalized with the power curve. Background spectra were recorded with an empty sample cup and subtracted from the sample data. Near-IR overtone spectra were collected from ∼50 mg pure, recrystallized TNT obtained from the Weapons Materials Research Directorate at the Army Research Laboratory and from ∼100 mg of munitions grade TNT obtained from the Pyrotechnics Branch at the Edgewood Chemical Biological Center. For preparation of TNT on silica, 12 mg of the munitions grade TNT was dissolved in 10 mL of acetone. Nine milliliters of this solution was added to 35 mg of Cabosil 5 M fumed silica, and the acetone was evaporated, leaving a coating of TNT on the fumed silica. According to the manufacturer, the fumed silica had a surface area of 200 m2/g of material, which means that there was 7.0 m2 of amorphous silica surface. Assuming TNT has a diameter of 7 Å and an area of 3.8  1019 m2, the 35 mg of TNT covered approximately 35 m2 in a monolayer. Assuming the TNT spread evenly, the layer of TNT on the silica surface was approximately five molecules deep. Because benzoic acid is also an aromatic solid, its overtone spectrum was recorded for the sake of comparison. The TNT overtone spectra were compared with the predicted spectra from HCAO. Saturation of the sample did not occur because plots of the integrated acoustic signal as a function of input laser power were measured to be linear.

’ RESULTS AND DISCUSSION As shown in Figure 2, both pure and munitions grade TNT spectra show the same features at 5811, 5865, and 5956 cm1 and the shoulders at 6010 and 6095 cm1. There is a small shift in the strongest feature from 6051 to 6063 cm1 and a reduction in intensity for the 5880 cm1 peak when going from pure TNT to the munitions grade TNT. For the TNT adsorbed onto silica, there are two strong features centered at 6039 and 5860 cm1. Most notably, there is the absence of the peak at 5956 cm1, which is present in the munitions grade and pure TNT. The spectrum of benzoic acid is distinguishable from TNT, although there is a similar dominant feature that is red-shifted to 6016 cm1, with no additional distinct features in the spectrum. The simulations of overtone spectra using single molecules, an 8 molecule cluster, and 64 molecule clusters (eight unit crystal cells) have elucidated three key findings on the effect of intermolecular interactions on the overtone spectrum. First, the simulations showed that the crystal structure in the classical layer of the model was critical. Second, electron correlation must be accounted for to correctly capture third-order force constants for 9142

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

Figure 4. Predicted spectra from single molecule models of the A and B forms reported from X-ray crystallography measurements. Figure 2. Near-IR laser photoacoustic spectra of pure, recrystallized TNT, munition grade TNT, TNT deposited with solvent on fumed silica powder, and benzoic acid in the frequency range of 57006700 cm1.

Figure 5. Overtone spectra produced from eight member TNT molecular clusters.

Figure 3. TNT overtone spectrum compared with predicted spectra from the local mode model for (1) HF/6-31+G(d,p), (2) DFT B3LYP/ 6-31+G(d,p), and (3) MP2/6-31+G(d,p). HCAO corrections to the predicted local mode frequencies are not shown in the predicted spectra.

the CH stretches. Finally, intermolecular forces must be accounted for by correct assignment of partial electric charge in both the quantum and classical layers of the model. The results of the single, 8, and 64 molecule cluster models illustrate how sensitive the simulated overtone spectra are to molecular placement. Single molecule models that do not include intermolecular interactions do not account for several features. When the spectra produced from single molecule models in Figure 3 are compared with the TNT data in Figure 2, it can be seen that the single molecule models for TNT do not account for the peak at 5956 cm1 or the shoulders around the main aryl CH stretch peaks at 6104 or 6017 cm1. Although the spectra from the single molecule models have the same general shape as the laser PAS data and in some cases come within 20 cm1 of the measured aryl CH stretches, there also is poor agreement for the methyl CH stretching frequencies and intensities. The absence of several features in the predicted spectra shown in

Figure 3 can be explained by degeneracy in the highly symmetric single molecule models. Application of anharmonic corrections to the normal-mode frequencies did not reproduce the additional features that appear in the experimental data. To investigate the influence of the two crystal polymorphs of TNT as identified by X-ray crystallography, a single molecule model of each polymorph (denoted as “A” and “B”) was optimized, and the HCAO analysis was applied. The nitro functional group dihedral angles of each respective model was set to match those measured with X-ray crystallography. The results are shown in Figure 4, where the predicted overtone spectra are insensitive to the intramolecular geometry. The intramolecular OH separations are comparable to the intermolecular OH separations for both simulation and the X-ray measurements. The insensitivity to the positions of the intramolecular nitro functional groups makes sense because the electrostatic interactions were directed almost orthogonally to the axes of the CH stretches, whereas the OH intermolecular interactions were directed more closely along the CH stretching axis, thus applying a greater perturbation. This result reinforces the fact that a multimolecular model was, indeed, necessary to interpret the experimental result. The results of the predicted overtone spectra from the eight member cluster model are shown in Figure 5. Although use of this simple cluster model eliminated the issue with degeneracy, 9143

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

Table 1. Predicted First Overtone Vibrational Frequencies CH stretching frequency(cm1) aryl stretches single molecule

methyl stretches

mode 1

mode 2

mode 3

mode 4

mode 5

HF/6-31+G(d,p) anharmonic correction

6580.3

6570.0

6349.3

6349.3

6362.5

B3LYP/6-31+G(d,p) anharmonic correction

6172.0

6161.5

5977.1

5960.6

5878.0

HF/6-31+G(d,p) local mode

6415.2

6415.2

6172.8

6171.6

6109.7

HCAO correction

6415.0

6415.0

6075.0

6068.8

6006.7

B3LYP/6-31+G(d,p) local mode

6039.7

6039.7

5825.6

5824.9

5750.5

HCAO correction

6039.4

6039.4

5728.1

5727.3

5660.1

MP2/6-31+G(d,p) local mode

6136.1

6135.3

5949.7

5949.0

5927.7

HCAO correction

6135.6

6134.8

5841.5

5835.9

5816.2

HF/6-31+G(d,p) anharmonic correction

6550.4

6538.5

6530.4

6371.5

6286.4

B3LYP/6-31+G(d,p) anharmonic correction

6138.4

6125.9

6123.0

5959.5

5954.2

B3LYP/6-31+G(d,p) molecule 1

6012.7

6005.5

5939.0

5762.5

5760.6

HCAO correction

6012.7

6005.4

5842.1

5661.8

5659.9

MP2/6-31+G(d,p) molecule 1

6144.0

6127.8

5997.6

5938.3

5923.6

HCAO correction

6143.6

6127.3

5884.2

5830.3

5809.7

MP2/6-31+G(d,p) molecule 2

6138.5

6145.9

6129.6

5953.3

5984.6

HCAO correction MP2/6-31+G(d,p) molecule 3

6137.9 6181.1

6145.3 6083.2

6017.0 5939.3

5838.1 5946.9

5870.4 6030.1

HCAO correction

6180.4

6082.5

5831.3

5833.7

5914.9

MP2/6-31+G(d,p) molecule 1

6147.9

6211.2

6010.9

5927.6

6175.5

HCAO correction

6147.3

6210.7

5900.8

5811.4

6057.1

MP2/6-31+G(d,p) molecule 2

6112.4

6138.9

6010.9

5997.2

6092.7

a

1 Unit Cell Cluster

8 Unit Cell Cluster

a

HCAO correction

6111.7

6138.2

5902.6

5882.5

5974.7

MP2/6-31+G(d,p) molecule 3

6189.4

6161.1

6128.1

6011.4

5926.4

HCAO correction MP2/6-31+G(d,p) molecule 4

6188.6 6085.4

6160.4 6146.3

6014.6 5942.7

5895.5 5993.5

5810.4 6095.1

HCAO correction

6084.6

6145.5

5834.8

5879.3

5977.1

TNT on silica

6147.3

6131.7

5936.6

5960.7

5961.5

HCAO correction

6147.3

6131.1

5828.9

5846.8

5847.6

aryl CH modes with anharmonic corrections to the normal modes refer to antisymmetric and symmetric stretches.

remaining differences between the predicted spectra and experiment show that the model inadequately captures the intermolecular interactions. The reason is that there was no position in the model that was completely caged by surrounding molecules, as occurs in the bulk crystal. It became clear that much improved results are obtained from the larger models that better simulate the TNT crystal. However, the eight molecule model yielded a key result. The results in Figure 5 from the eight TNT molecule model show that the handling of electron correlation in the two-layer QM:MM models is critical. Figure 5 is a comparison of predicted overtone spectra from the MP2/6-31+G(d,p) and B3LYP/6-31 +G(d,p) levels of theory. When B3LYP/6-31+G(d,p) was used in the hybrid QM:MM eight molecule cluster model, it predicted a large drop in overtone intensity compared with the single molecule spectrum predicted in Figure 3 such that the DFTpredicted spectrum was attenuated by a factor of ∼33. In contrast, the absorption intensities of the MP2/6-31+G(d,p)predicted spectrum in Figure 5 is in reasonable agreement with

the predicted intensities of all three levels of theory in Figure 3. The fact that the MP2-predicted spectra are not attenuated in the hybrid QM/MM models is quantitatively shown in Table 3. The oscillator strengths predicted from the single unit cell model and the MP2 method agree with the single molecule oscillator strengths to within an order of magnitude. Comparison of the HCAO predicted Morse parameters in Table 2 provide additional insight into the reason the MP2 method gave better results. For the single molecule models, the Morse parameters were roughly equivalent among the three levels of theory, where anharmonicities of 8387 cm1 were typical. The results obtained from a model cluster with the MP2 method were about the same as from single molecule models. The case is different for the B3LYP/6-31+G(d,p) level of theory, for which the anharmonicities reduced to 0.010.02 cm1 for the aryl CH stretches, and 15 cm1 for the methyl CH stretches. As expected, the overtone absorption intensities became correspondingly smaller, and the shape of the absorption spectrum became radically different. A plot of the PES predicted by the 9144

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

Table 2. Local Mode Morse Parameters aryl stretches mode 1 1

B3LYP/6-31+G(d,p) ωe (cm )

3229.2

anharmonicity ωexe (cm1)

methyl stretches mode 2

mode 3

mode 4

mode 5

3229.2

3094.5

3131.1

3131.3

83.73

MP2/6-31+G(d,p) ωe (cm1)

83.73

3277.8

anharmonicity ωexe (cm1)

87.72

3277.5

83.89

3189.7

83.93

87.45 3189.7

85.94

85.95

87.41 3178.8 85.95

1 Unit Cell Molecular Cluster B3LYP/6-31+G(d,p) molecule 1, ωe (cm1)

3006.4

anharmonicity ωexe (cm1)

3002.8

0.01

2971.7

0.02

2893.2

0.87

4.77

2892.1 4.72

MP2/6-31+G(d,p) molecule 1, ωe (cm1)

3281.3

3274.0

3213.9

3186.2

3182.1

anharmonicity ωexe (cm1) MP2/6-31+G(d,p) molecule 2, ωe (cm1)

83.74 3282.9

84.03 3284.1

86.03 3275.8

86.82 3199.7

88.13 3211.6

anharmonicity ωexe (cm1)

85.47

MP2/6-31+G(d,p) molecule 3, ωe (cm1)

84.46

3316.5

anharmonicity ωexe (cm1)

84.39

3253.7

90.39

3185.1

84.84

89.24 3187.8

86.20

85.73

87.71 3226 84.38

8 Unit Cell Molecular Cluster MP2/6-31+G(d,p) molecule 1, ωe (cm1)

3295.9

anharmonicity ωexe (cm1)

3384.5

3234.7

3183.8

3299.7

88.78

111.56

91.71

87.99

MP2/6-31+G(d,p) molecule 2, ωe (cm1) anharmonicity ωexe (cm1)

3288.7 92.99

3293.7 89.71

3234.7 91.71

3241.5 97.17

3274.3 91.16

MP2/6-31+G(d,p) molecule 3, ωe (cm1)

3305.9

3325.5

3273.7

3245

3187.6

anharmonicity ωexe (cm1)

84.492

MP2/6-31+G(d,p) molecule 4, ωe (cm1)

97.99

3277.5

anharmonicity ωexe (cm1)

83.87

3296.9

93.9

3212.2

89.493

95.72 3240.1

96.32

97.36

84.8

89.77 3275.1 91.03

TNT on Silica MP2/6-31+G(d,p) ωe (cm1)

3283.8

anharmonicity ωexe (cm1)

3274.4

83.94

3184.1

83.42

3201.4

86.34

88.41

3199.1 87.35

Table 3. Predicted CH Stretching Overtone Oscillator Strengths aryl CH stretches

methyl CH stretches

mode 1

mode 2

mode 3

mode 4

mode 5

4.58  108

4.58  108

8

2.13  108

4.83  108

2.13  108

6.74  10

8

9.08  10

8

8

2.19  108

5.15  10

8

1.72  10

8

8

6.47  10

1.72  108

Single Molecule HF/6-31+G(d,p) B3LYP/6-31+G(d,p)

6.74  10

8

2.19  10

MP2/6-31+G(d,p)

5.17  10

B3LYP/6-31+G(d,p) molecule 1

5.13  109

5.10  109

7.21  108

3.51  108

3.13  108

7

8.82  10

7

8

7

4.54  109

3.03  10

8

7

7.17  107

5.31  10

8

7

2.08  10

4.21  109

1 Unit Cell Molecular Cluster MP2/6-31+G(d,p) molecule 1 MP2/6-31+G(d,p) molecule 2

8.91  10

8

6.77  10

8

7.81  10 1.69  10

7

1.83  10

7

5.81  10

4.19  10

MP2/6-31+G(d,p) molecule 3

6.35  10

MP2/6-31+G(d,p) molecule 1

4.76  108

3.60  108

1.10  108

6.88  108

2.51  108

8

1.13  10

7

7

1

4.50  109

1.45  10

8

8

3.99  108

7

1

3.43  10 0

1.32  108

2.53  108

2.42  108

8 Unit Cell Molecular Cluster MP2/6-31+G(d,p) molecule 2 MP2/6-31+G(d,p) molecule 3

4.27  10

8

7.15  10

9

MP2/6-31+G(d,p) molecule 4

3.59  10

1.25  10

MP2/6-31+G(d,p)

4.69  108

3.41  108

1.17  10 2.23  10

8

5.13  10

8

6.37  10 0 2.87  10

TNT on Silica

different levels of theory for the single and cluster models shows what happened.

3.78  108

Figure 6 shows the HF/6-31+G(d,p)-, B3LYP/6-31+G(d,p)-, and MP2/6-31+G(d,p)-predicted aryl CH stretch PES for 9145

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

Figure 6. Plot of potential energy surface along the aryl CH stretch coordinate for TNT within an eight molecule cluster. Neither the B3LYP or X3LYP functionals nor HF methods reproduced the anharmonic potential energy surface, whereas the surface calculated with MP2 retained anharmonic character.

both single molecule and cluster models. A PES from an additional calculation using the X3LYP/6-31+G(d,p) is also presented because this functional is designed to better handle correlation effects.31 For all levels of theory on a single molecule model, there is a clear third- or higher-ordered component to the PES that contributes to the cubic force constant used to calculate the anharmonicity. For most of the cluster models, the PES is primarily parabolic, and only the PES from the MP2/6-31+G(d,p) level of theory keeps the higher order component that resembles the single molecule models. A consequence of the loss of the cubic component is that the anharmonicity is greatly reduced, and overtone intensity is lost. The results of both B3LYP and HF differ from MP2-predicted results because when these two approaches solve for single electron wave functions, they treat the other electrons as an average charge density, that is, a mean field approach. For the hybrid classical/quantum mechanical models used in this study, the embedding of static point charges in the classical layer interferes with the mean field approach. In other words, when the single electron wave functions are calculated, the mean field includes contributions from charges in the MM layer of the model, and it no longer describes the behavior of just the electrons within the QM layer. The charges in the MM layer differ from the electrons in the QM layer because they are not mobile, but by allowing them to contribute to the mean field, they are treated the same as the electrons in QM layer, and this is not physical. This conclusion is further supported by the observed trend that the anharmonicity reduced further as the MM layer was expanded with additional charges. In contrast, the MP2 method produces predicted spectra that are in much better agreement with experiment, and the method maintains reasonable values for anharmonicities. The primary difference among the MP2, the HF, and B3LYP methods is that only the MP2 method provides a correction for electron correlation energies. This result is consistent with published results3133 that show that accounting for electron correlation is necessary to model weak forces, such as van der Waals or other intermolecular interactions. It is expected that multiconfigurational methods will produce better agreement with experiment because they account for electron correlation more accurately.

Figure 7. The laser photoacoustic overtone spectrum compared with predicted spectra from HCAO and the MP2/6-31+G(d,p) level of theory for charge assigned by the QEq/Mulliken, QEq/MSK, and HF/ 3-21+G* calculation/MSK. The model is optimized with the universal force field.

To further investigate the role of electron correlation effects, the cluster model was optimized, and a potential energy surface scan was performed with the X3LYP functional. The X3LYP functional is an extended version of B3LYP that is designed to take into account electron correlation effects such as longer range interactions, and it agrees with experimental measurements of binding energies of noble gas dimers. However, when the X3LYP calculation was permitted to interact with the partial charges in the molecular mechanics layer of the cluster, the resulting PES had negligible third-order character and was similar to the B3LYP predicted result. A primary difference between the two layer study here and the systems used to develop the X3LYP functional is that those model systems are treated entirely with DFT. Thus, in those systems the intermolecular interactions are not described by an interaction between a quantum mechanical system and point charges. As a result, the X3LYP functional could not correctly handle the correlation effects imposed by the point charges in the MM layer. Calculations that were attempted using X3LYP to describe the target molecule and its nearest neighbors did not converge within a reasonable time. The effect of intermolecular forces on the model spectra is reflected in the sensitivity to charge assignment in both the QM and MM layers of the model. The HCAO spectra (Figure 7) that were calculated from the large cluster model shown in Figure 1C illustrate this. The predicted spectra from different combinations of methods for assigning electric charge are shown where the scheme that uses an ab initio calculation to assign charge in the low layer and the MerzSinghKollman population analysis to assign charge in the high layer yielded the best results. The first predicted spectrum in Figure 7 is representative of models using the Mulliken population analysis in the determination 9146

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

Figure 8. Molecular dipole moment as a function of CH bond length for the Mulliken and the MerzSinghKollman schemes, from both a single molecule model and cluster models.

of molecular dipole moment, and it has poor agreement with the experimental spectrum. Neither performing the HCAO analysis on different positions within the cluster nor extending the cluster size improved agreement with the data, except for the position where the aryl hydrogens were directed outside of the cluster. As the model cluster size increased, the intensity of the low-frequency peaks generally increased relative to the high-frequency peaks. Although the spectra from single molecule models matched for both the Mulliken34 and the MerzSinghKollman (MSK)35 schemes, only the results based on fitting partial charges to the molecular electronic potential, such as the MSK or CHelp36 methods, resembled the experimental data. This remained true, even as the size of the MM layer was increased. The plot of molecular dipole moments as a function of CH bond length in Figure 8 helps explain this trend. For the single molecule model of TNT, the slopes of the dipole moment versus aryl CH bond length predicted by the Mulliken and MSK schemes were close in value. For the HCAO calculation of transition dipole matrix elements, the derivatives of the dipole moment function contributed a major portion to the absorption intensity along with the overlap of nuclear wave functions; thus, similar oscillator strengths were obtained for both population analyses. For the QM/MM cluster models, the dipole moment as a function of CH bond length for the MSK scheme had a slope similar to that of the single molecule models. This was true although the static dipole moment was higher. In contrast, the slope of the dipole moment as predicted from the Mulliken population analysis was reduced. Clearly, the Mulliken analysis did not handle the point charges in the MM layer well. The Mulliken population analysis assigns partial charge based on the weighting constants of the atomcentered basis sets as determined from the self-consistent-field (SCF) calculations. Although the SCF-determined wave functions responded to the surrounding charges in the model, partial charges calculated from these wave functions did not sufficiently change in response to the MM layer. In contrast, methods based on fitting the atom centered partial charges to match the molecular electric potential, such as MSK or CHelp, responded better to the

ARTICLE

surrounding point charges. This makes sense in view of the fact that these calculations must directly take into account surrounding point charges when fitting to match the potential. It should be noted that the two schemes predicted the same dipole moment at the equilibrium geometry for both the single molecule and cluster models, and the permanent dipole moments predicted by these two approaches alone did not reveal the reason for the different results in the HCAO predicted spectra. The assignment of charge in the MM layer of the model is almost equally important to correctly represent the intermolecular forces acting on the QM layer. The empirically based QEq algorithm is used for the previously discussed spectra to assign charge in the MM level of the model. For the last spectrum in Figure 7, an MSK population analysis calculation at the HF/3-21 +G* level of theory is performed to assign partial charge in the MM layer, and the HCAO analysis is applied to the resulting system. As can be seen, the resulting spectrum bears a strong resemblance to the data, indicating that the model has captured more of the physics of the experiment. Because the QEqgenerated charges are used for the optimization step in this model, the model is refined further to improve the molecular placement and bring agreement with experiment. The model used to generate the spectra in Figure 9 was arranged to match the intermolecular distances reported in the X-ray crystallography study. Charges were then assigned with the MSK population analysis of the result of an HF/3-21+G* calculation, and a subsequent optimization with the UFF was performed. Spectra were obtained from four positions within the model. As shown in the Figure 9, the four individual spectra have significant differences. The intensities and frequencies of the aryl CH stretches are similar, but the locations of the most prominent methyl CH stretch absorptions are separated by ∼90 cm1. The differences between position 2 and the others are even more pronounced, where the aryl CH absorptions occur at higher frequencies and the methyl CH have an intensity comparable to the aryl CH stretches. For each position sampled within the cluster, there is an almost distinct spectrum. However, when the spectra from all positions are averaged together, the resulting spectrum resembles the experimental spectrum, including many small features, such as the shoulders at 6100 and 6130 cm1 or the peak at 5930 cm1. Two observations explain the variation in predicted spectra from each position in the model. As stated previously, the intermolecular interactions break the molecular symmetry and lead to the variety of features in the experimental spectrum. The A and B crystal polymorphs do not contribute to the variety of features in the spectrum. Second, because the crystal unit cell has eight distinct positions, it is expected that the experimental spectrum is an average of all the spectra from each position. Although the general form of the intensities and frequencies of the spectrum in Figure 9 resembled the data, perfect agreement was not obtained. Only some of the intermolecular distances were reported from X-ray data, so the TNT cluster model could not completely replicate the physical crystal. Improved agreement with experiment is expected from a TNT geometry that is optimized either with a MM or QM calculation with periodic boundary conditions. To test the ability of this approach to predict overtone spectra from molecules on surfaces, the output of a single molecule model on a fully hydrated silica surface was compared with the measured spectrum of a thin film of TNT adsorbed onto fumed silica shown in Figure 10. The difference between the experimental 9147

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

Figure 9. Spectra from the different positions in Figure 1 C, including the sum of all the spectra and the experimental spectrum from pure TNT.

Figure 10. Single molecule spectrum of TNT on deposited on silica.

spectrum on silica and the pure TNT and munitions grade TNT in Figure 2 showed a change in the intermolecular placement of the TNT molecules. The spectrum became greatly simplified so that only two features were discernible. The features in the predicted overtone spectrum in Figure 10 resemble the experimental features. This indicated that the local environment experienced by the TNT molecule in the model was similar to that in the experimental measurement. However, it is likely that the real system, as with the pure TNT results, was an average over an ensemble of molecules. Thus, although the predicted and observed spectra resembled each other, to conclude whether the modeling approach was consistent with experiment, an orientation-averaged ensemble of molecules on a silica surface is necessary.

Finally, with respect to the HCAO corrections to the local mode frequencies, the inclusion of harmonic coupling between anharmonic oscillators did not bring additional agreement with experiment. For the simulated spectra in Figures 3, 4, 5, 7, 9, and 10, only the pure local mode predicted frequencies are represented. As shown throughout Table 1, the HCAO corrections for the CH methyl oscillators were on the order of ∼100 cm1. Whenever these corrections were applied, differences with the experimental spectra were increased. This indicated that the effect of the surrounding molecules in the model far outweighed the effect of interaction between intramolecular CH oscillators. This was consistent with the observation that intramolecular interactions between the hydrogen atoms and the nitro functional groups were insufficient to explain all the features in the experimental spectrum. Thus, a pure local mode treatment appeared to be appropriate for molecules interacting with a matrix of surrounding molecules. It is difficult to draw definitive conclusions on how both intraand intermolecular geometries were reflected in the predicted overtone spectra. A comparison of bond lengths and intermolecular distances from the cluster model in Figure 1C (see Table 4), which yielded the best agreement between predicted spectra and the experiment, and the X-ray data, did not reveal great differences. For most bond lengths within the skeleton of the molecule (e.g., CC and CN bonds, not shown), the models came to within 1% of both the A and B crystal polymorphs, and there was no discernible pattern between a given bond length and features in the overtone spectrum. For the nitro functional group dihedral angles in Table 4 (the atom labels for Table 4 are in Figure 1A), there was deviation from the pattern that can be correlated to the overtone spectrum. In terms of dihedral angles, model positions 1 and 3 tended to be similar, and positions 2 and 4 were also similar to each other. Positions 1 and 3 had a large negative angle for the N1 nitro functional group, which was adjacent to the methyl group, a small positive angle for the N3 nitro functional group adjacent to the methyl group, and a large negative angle for the N2 nitro group. For model positions 2 and 4, the N1 and N3 nitro groups had large positive dihedral angles and a small 9148

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

Table 4. Intramolecular Distances between Hydrogens and Oxygens, And Dihedral Angles for the Nitro-Functional Groups, Compared to Distances Reported for A and B Forms from X-ray Dataa atomic distances (Å)

nitro group dihedral angles (deg)

H1O1

O2H4

H4O3

O4H5

C4C3N1-O2

C6C5N2-O4

C6C7N3-O5

X-ray TNT-A

2.39

2.52

2.57

2.43

43.0

24.0

51.0

X-ray TNT-B

2.62

2.52

2.57

2.63

45.0

30.0

60.0

HF/6-31+G(d, p)

2.42

2.46

2.40

2.40

37.8

0.1

37.8

B3LYP/6-31+G(d, p)

2.41

2.50

2.42

2.42

37.9

0.2

37.9

Single Molecule Models

MP2/6-31+G(d, p)

2.56

2.69

2.42

2.42

52.6

0.3

52.7

constrained TNT A MP2/6-31+G(d, p)

2.11

2.65

2.50

2.50

43.0

24.0

51.0

constrained TNT B MP2/6-31+G(d, p)

2.52

2.59

2.54

2.54

45.0

30.0

60.0

position 1

2.45

2.52

2.42

2.42

38.4

2.1

56.7

position 2

2.34

2.77

2.64

2.62

61.4

39.3

25.0

position 3

2.27

2.62

2.48

2.50

52.2

23.3

40.1

position 1

2.44

2.57

2.45

2.40

33.3

0.3

32.1

position 2

2.21

2.33

2.52

2.70

51.2

38.3

14.5

position 3

2.33

2.63

2.50

2.46

43.6

12.3

21.9

position 4

2.33

2.58

2.56

2.52

52.6

42.6

13.8

Eight-Membered Cluster

64 Membered Cluster

a

Atom labels are from Figure 1A.

dihedral angle on the N2 nitro group. In the overtone spectra, we see that molecular positions 1 and 3 contributed significantly in intensity to the peak at ∼5920 cm1 and also to the shoulder at ∼6180 cm1. Positions 2 and 4 contributed strongly to the peak at ∼6140 cm1. Breaking the grouping pattern, positions 1, 2, and 3 contributed to the feature at ∼6010 cm1. A comparison of intermolecular distances in the model and distances from the published X-ray crystal showed no discernible correlation between the predicted spectra and these distances, and all the distances were close to those measured by X-ray crystallography. Most likely, interactions between a CH stretch and multiple nearest neighbors is necessary to directly demonstrate the cause and effect.

’ CONCLUSIONS We investigated the CH overtone spectrum of neat TNT and TNT deposited on fumed silica by laser photoacoustic spectroscopy and modeling. The local mode model was applied to the CH oscillators in the TNT molecule, including HCAO corrections to the local mode frequencies. Single TNT molecules, 8 member TNT clusters, and 64 molecule TNT clusters yielded reasonable results; the large cluster models yielded the best results. Comparison of the overtone spectra calculated from the single molecule models and 8 molecule cluster models showed agreement with some features in the experimentally measured spectra, specifically, the aryl CH stretches; however, they were inadequate because they did not account for many features in the experimental spectra. The 64 molecule model that incorporated interaction between the QM layer and point charges in the MM layer yielded the best results between model and experiment. For the two-layer QM:MM approach, both the treatment of electron correlation and assignment of charge are critical.

Predicted spectra showed that the overtone spectrum was very sensitive to the intermolecular forces acting on the molecule undergoing HCAO analysis. A consequence of this was that frequencies, intensities, and general appearance of the spectra could change drastically in response to intermolecular distances in the model as well as to how partial electrostatic charges were assigned. Clearly, agreement between theoretical and experimental spectra indicates correct modeling of the intermolecular forces. In conclusion, the HCAO/local mode approach in combination with hybrid QM:MM cluster models of condensed phase TNT predicts reasonable spectra. Furthermore, it is shown that modeling the intermolecular interactions as electrostatic interactions accounts for most of the significant features in the experimental measurement. Thus, the approach described here is reasonable for investigating intermolecular interactions.

’ AUTHOR INFORMATION Corresponding Author

*Address: 5183 Blackhawk R., E5951, Aberdeen Proving Ground, MD 21010. Phone: 410-436-5634. Fax: 410-4362742. E-mail: [email protected].

’ ACKNOWLEDGMENT The authors thank Dr. Barry Williams and Dr. Ronald Miles for their assistance and Dr. Rose Pesce-Rodriguez of the Army Research Laboratory for providing purified TNT samples. We also acknowledge support from the Edgewood Chemical Biological Center’s In House Independent Research program for supporting this work. ’ REFERENCES (1) Roddick-Lanzilotta, A. D.; Connor, P. A.; McQuillan, A. J. Langmuir 1998, 14 (22), 6479–6484. 9149

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150

The Journal of Physical Chemistry A

ARTICLE

(2) Johnson, W. G.; Buch, V.; Trenary, M. J. Chem. Phys. 1990, 93 (12), 9167–9179. (3) Sung, C. Y.; Broadbelt, L. J.; Snurr, R. Q. J. Phys. Chem. C 2009, 113 (35), 15643–15651. (4) Bermudez, V. M. J. Phys. Chem. C 2009, 113 (5), 1917–1930. (5) Tao, Y. T.; Lin, W. L.; Hietpas, G. D.; Allara, D. L. J. Phys. Chem. B 1997, 101 (47), 9732–9740. (6) Fujita, S.; Sato, H.; Kakegawa, N.; Yamagishi, A. J. Phys. Chem. B 2006, 110 (6), 2533–2540. (7) Lundell, J.; Chaban, G. M.; Gerber, R. B. J. Phys. Chem. A 2000, 104 (34), 7944–7949. (8) Gerber, R. B.; Ratner, M. A.; Bacic, Z. Tunneling; Jortner, J., Pullman, B., Eds.; D. Reidel: Dordrecht, The Netherlands, 1986; pp 2538. (9) Gerber, R. B.; Horn, T. R.; Williams, C. J.; Ratner, M. A. Dynamics of Polyatomic van der Waals Complexes; Halberstadt, N., Janda, K. C., Eds.; Plenum Press: New York, 1990; pp 343354. (10) Ratner, M. A.; Gerber, R. B. J. Phys. Chem. 1986, 90 (1), 20–30. (11) Kjaergaard, H. G.; Henry, B. R.; Wei, H.; Lefebvre, S.; Carrington, T., Jr.; Mortensen, O. S.; Sage, M. L. J. Chem. Phys. 1994, 100 (9), 6228–6239. (12) Kjaergaard, H. G.; Henry, B. R. J. Chem. Phys. 1992, 96 (7), 4841–4851. (13) Petryk, M. W. P.; Henry, B. R. J. Phys. Chem. A 2005, 109 (18), 4081–4091. (14) Petryk, M. W. P. Prog. Biomed. Opt. Imaging 2006, 7 (43), 637818.1–637818.14. (15) Carper, W. R.; Davis, L. P.; Extine, M. W. J. Phys. Chem. 1982, 86 (4), 459–462. (16) Gaussian 2009, Revision A.02, Gaussian, Inc.: Wallingford CT, 2009. (17) Child, M. S.; Halonen, L. Adv. Chem. Phys. 1984, 57, 1–58. (18) Henry, B. R.; Tarr, A. W.; Mortensen, O. S.; Murphy, W. F.; Compton, D. A. C. J. Chem. Phys. 1983, 79 (6), 2583–2589. (19) Henry, B. R. Acc. Chem. Res. 1977, 10, 207–212. (20) Sowa, M. G.; Henry, B. R.; Mizugai, Y. J. Phys. Chem. 1991, 95, 7659–7664. (21) Child, M. S.; Lawton, R. T. Faraday Discuss. Chem. Soc. 1981, 71, 273–285. (22) Miller, W. H.; Hernandez, R.; Handy, N. C.; Jayatilaka, D.; Willetts, A. Chem. Phys. Lett. 1990, 172, 62–68. (23) Clabo, D. A.; Allen, W. D.; Remington, R. B.; Yamaguchi, Y.; Schaefer, H. F., III Chem. Phys. 1988, 123, 187–239. (24) Barone, V. J. Chem. Phys. 2005, 122, 014108. (25) Barone, V. J. Chem. Phys. 2005, 120, 3059–3065. (26) Wilson, E. B. J. Chem. Phys. 1939, 7, 1047–1052. (27) Decius, J. C. J. Chem. Phys. 1948, 16 (11), 1025–1034. (28) Howard, D. L.; Jorgensen, P.; Kjaergaard, H. G. J. Am. Chem. Soc. 2005, 127 (48), 17096–17103. (29) Kjaergaard, H. G.; Henry, B. R. J. Phys. Chem. 1995, 99, 899–904. (30) Brabson, G. D. J. Chem. Educ. 1973, 50 (6), 397. (31) Xu, X.; Goddard, W. A., III Proc. Natl. Acad. Sci. 2004, 101 (9), 2673–2677. (32) Wodrich, M. D.; Cormeinboeuf, C.; Schleyer, P. v. R. Org. Lett. 2006, 8 (17), 3631–3634. (33) Kone, M.; Illien, B.; Graton, J.; Laruence, C. J. Phys. Chem. A 2005, 109 (51), 11907–11913. (34) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833–1831. (35) Besler, B. H.; Merz, K. M., Jr.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431–439. (36) Chirlian, L. E.; Francl, M. M. J. Comput. Chem. 1987, 8, 894–905.

9150

dx.doi.org/10.1021/jp2025895 |J. Phys. Chem. A 2011, 115, 9139–9150