Role of the Multipolar Electrostatic Interaction Energy Components in

An analysis of the charge sharing between the donor (the π-systems) and the acceptors (the cations) together with the partitioning of total interacti...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/JPCA

Role of the Multipolar Electrostatic Interaction Energy Components in Strong and Weak Cation−π Interactions Pawel Kadlubanski,† Katherine Calderón-Mojica,‡,⊥ Weyshla A. Rodriguez,§,⊥ D. Majumdar,*,∥ Szczepan Roszak,*,† and Jerzy Leszczynski*,∥ †

Institute of Physical and Theoretical Chemistry, Wroclaw University of Technology, Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland ‡ Universidad de Puerto Rico en Humacao, Estacidn Postal CUH 100 Carr.908, Humacao PR00791-4300, Puerto Rico § Universidad Metropolitans, P.O. Box 21150, San Juan 00928-21150, Puerto Rico ∥ Department of Chemistry and Biochemistry, Jackson State University, Jackson, Mississippi 39217, United States S Supporting Information *

ABSTRACT: Density functional and Møller−Plesset secondorder perturbation (MP2) calculations have been carried out on various model cation−π complexes formed through the interactions of Mg2+, Ca2+, and NH4+ cations with benzene, pmethylphenol, and 3-methylindole. Partial hydration of the metal cations was also considered in these model studies to monitor the effect of hydration of cations in cation−π interactions. The binding energies of these complexes were computed from the fully optimized structures using coupled cluster calculations including triple excitations (CCSD(T)) and Gaussian-G4-MP2 (G4MP2) techniques. An analysis of the charge sharing between the donor (the π-systems) and the acceptors (the cations) together with the partitioning of total interaction energies revealed that the strong and weak cation−π interactions have similar electrostatic interaction properties. Further decomposition of such electrostatic terms into their multipolar components showed the importance of the charge− dipole, charge−quadrupole, and charge−octopole terms in shaping the electrostatic forces in such interactions. The computed vibrational spectra of the complexes were analyzed for the specific cation−π interaction modes and have been shown to contain the signature of higher order electrostatic interaction energy components (quadrupole and octopole) in such interactions.

I. INTRODUCTION Cation−π interactions are noncovalent in nature and were first predicted by Kier and coworkers1,2 to be of ion-induced dipole character (from theoretical modeling of acetylcholinesterase inhibition reactions). These interactions were later established as cation−π type by Dougherty and coworkers,3,4 and since then it has emerged as a very common binding force to interpret interactions between cations and aromatics in biomolecular systems such as protein, receptor−ligand complexes, molecular recognition, drug action, and protein folding.3−9 The cation−π interactions are strong enough to compete with the solvation of hydrophilic charged moieties and allow ligand−receptor binding in a hydrophobic pocket constituted of aromatic residues within the core of the receptor.8 Cation−π interactions could be classified into weak and relatively stronger categories. The weak interactions of this class in biological systems usually involve a charged amino/OH2+ group and the localized π-electron cloud of the aromatic amino acids (viz. Phe, Trp, or Tyr).4 Several theoretical calculations as well as experiments have shown that such interactions are mostly on the order of 14−20 kcal/mol (depending on the nature of the cationic part), and although further weaker in © 2013 American Chemical Society

aqueous medium these interactions could still be on the order of 5 to 10 kcal/mol depending on the nature of hydration.5,10 The stronger cation−π interactions involve metal ions instead of charged amino/OH2+ groups, and for monovalent alkali metals these interactions could be as strong as 45 kcal/mol.4,11 In recent times the interactions between H2B2 (3Σg−) and metal ions have been interpreted as cation−π type, and binding energies in such cases have been computed to be several orders of magnitude higher than regular cation−π interactions of metal cations and aromatic systems.12 The cation−π interactions involving metal dications (e.g., Ca2+, Mg2+, Zn2+, Cu2+, Fe2+, Mn2+, etc.) are more recent developments and have been found to be important in the structure-modulation, reactivity, and functions of the metal binding sites of metalloproteins.13,14 Several divalent cations (e.g., Ca2+, Mg2+) have also been found to stabilize the unstacked conformations of DNA and RNA through cation−π interactions with the bases.15 Cation−π interactions are mostly biological events occurring in proteins or nucleic acid systems and are assumed to occur Received: April 29, 2013 Revised: July 25, 2013 Published: July 29, 2013 7989

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

Table 1. Calculated Interacting Distances (r, Å), Binding Energy (ΔEB), Enthalpy (ΔHB298), and Gibbs Free Energy (ΔGB298) of Binding (all in kcal/mol) of Various Cation−π Complexes at the DFT/B3LYP, MP2, CCSD(T), and G4MP2 Levels of Theorya,b cation−π complex benzene-Mg2+

methods

rc,d

ΔEBd

ΔHB298d,e

ΔGB298d,e

DFT/B3LYP

1.95 (1.94) 1.96 (1.98) 1.96 1.92 1.94 1.95 1.95 1.91 1.90 1.91 1.91 1.87 2.46 (2.36) 2.47 (2.46) 2.47 2.35 2.45 2.45 2.45 2.32 2.40 2.40 2.40 2.27 3.02 (3.00) 2.99 (2.89) 2.99 2.92 3.09 2.99 2.99 2.90 2.99 2.90 2.90 2.87

−117.76 (−117.81) −108.75 (−109.66) −109.00 −117.97 −129.68 −119.51 −119.51 −129.14 −145.67 −140.18 −135.26 −145.30 −82.90 (−81.59) −68.60 (−69.43) −68.53 −80.20 −94.41 −74.81 −74.20 −89.07 −103.96 −90.15 −89.53 −103.79 −14.95 (−15.05) −16.95 (−17.46) −16.06 −17.16 −16.84 −18.95 −17.95 −19.26 −21.30 −23.70 −22.50 −23.91

−118.36 (−118.40) −109.71 (−110.32)

−110.32 (−108.90) −100.60 (−100.72)

−118.56 −130.45 −120.63

−109.40 −121.22 −110.88

−129.74 −146.82 −136.72

−119.44 −137.81 −126.83

−145.89 −81.44 (−82.26) −69.37 (−69.93)

−137.26 −73.39 (−73.06) −59.55 (−60.71)

−80.80 −90.98 −75.61

−71.99 −82.28 −66.52

−89.66 −104.39 −91.12

−79.84 −96.23 −81.73

−104.38 −15.71 (−14.78) −16.97 (−17.32)

−96.07 −5.40 (7.47) −8.91 (−11.04)

−17.76 −18.28 −19.03

−10.29 −5.63 −9.98

−19.28 −21.16 −23.90

−9.71 −12.72 −14.60

−24.50

−15.51

MP2

p-methylphenol-Mg2+

3-methylindole-Mg2+

benzene-Ca2+

CCSD(T) G4MP2 DFT/B3LYP MP2 CCSD(T) G4MP2 DFT/B3LYP MP2 CCSD(T) G4MP2 DFT/B3LYP MP2

p-methylphenol-Ca2+

3-methylindole-Ca2+

benzene-NH4+

CCSD(T) G4MP2 DFT/B3LYP MP2 CCSD(T) G4MP2 DFT/B3LYP MP2 CCSD(T) G4MP2 DFT/B3LYP MP2

p-methylphenol-NH4+

3-methylindole-NH4+

CCSD(T) G4MP2 DFT/B3LYP MP2 CCSD(T) G4MP2 DFT/B3LYP MP2 CCSD(T) G4MP2

Calculated energies (i.e., ΔEB, ΔHB298, and ΔGB298) are after BSSE and zero-point energy (ZPE) corrections. bEnergy values at the DFT/B3LYP, MP2, and CCSD(T) levels are computed using aug-cc-pvdz basis set of atoms (see the text for details). cDistance r is defined as the distance between the cation/cation center and the center of mass of hexagonal carbon ring system (excluding hydrogen or any substituent, see Figure 1). dValues for fully optimized structures (using aug-cc-pvdz basis set) are presented within parentheses. eSuperscript 298 indicates the temperature (K) of ΔHB298 and ΔGB298 calculations. a

due to attraction of the quadrupole created by π-electron clouds of the conjugated ring systems and the cations.4,5,16 Various models are available to interpret the nature of electrostatics in cation−π interactions.4 Kim and coworkers16 made a significant contribution in this direction through termby-term computation of the multipolar electrostatic interaction energy components in benzene···ammonium (and its methyl derivatives) systems (up to charge-quadrupolar term) and found that these interaction energy components are important in cation−π interactions. This was also supported through

subsequent studies, although nonelectrostatic energy components are also important.4,5 The total electrostatic interactions could be partitioned as the contribution of charge−charge, charge-dipole, charge−quadrupole, charge−octopole (considering the cationic part to be composed of point charge systems), and higher order terms.17 Although it has been shown that the quadrupole term of the π-systems is responsible for such interactions,4,16 the role of the other electrostatic terms are not properly analyzed. Furthermore, if these electrostatic terms have significant contribution in such interactions, do they 7990

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

the ideal gas, rigid rotator, and harmonic oscillator approximations.26 The computed energy values of the gas-phase clusters were further refined using single-point energy calculations on the optimized structures at the DFT/B3LYP and MP2 levels using the aug-cc-pvdz basis set of atoms. The respective ΔHB298 and ΔGB298 values were approximated from the computed values at the DFT/B3LYP/6-31+g(d) and MP2/6-31+g(d) levels using the following relation 1.

behave similarly in both the strong and weak interaction cases? Could the electrostatic behavior of such cation−π complexes be connected to any observable properties (molecular/spectroscopic)? The aim of the present manuscript is to seek answers to these questions using several model cation−π complexes (both strong and weak types). The estimation of electrostatic contribution in our investigation would require reliable computation of binding energies and vibrational frequencies of the various cation−π complexes. In the present work, we have designed our model π-electron systems as benzene, p-methylphenol, and 3-methylindole. Benzene is the reference molecule because several cation−π interaction calculations have been carried out on this standard π-system and the last two molecules are the standard model aromatic systems for tyrosine and tryptophan.5,18,19 The substituent effect on the π-cloud and its consequence on cation−π interactions could also be visualized through such model molecules. The interactions of these aromatic systems were investigated against NH4+, Mg2+, and Ca2+ ions using density functional (DFT)20 and Møller−Plesset second order perturbation (MP2)21 theories. The strong cation−π interactions were further investigated using the hydrated forms of the metal cations. We have computed the structure, vibrational characteristics, binding energies, and related thermodynamic quantities of binding (viz., Gibbs free energy and enthalpy) at these levels of theory and have further validated the results using G4MP2 and coupled cluster computations. The interactions (both gas-phase and hydrated forms) are primarily analyzed through a hybrid variational-perturbational interaction energy decomposition scheme at the Hartree−Fock level.17,22 The electrostatic interaction energy terms of the gas-phase structures are further analyzed to generate the contribution of the multipolar electrostatic energy terms. In the final part of the manuscript, the vibrational spectra of the cation−π complexes are analyzed, and it would be shown that the modes representing cation−π interactions could be correlated with the nature of electrostatic forces operating to stabilize such complexes.

B B ΔX 298 = ΔX 298 (6‐31+g(d)) + [ΔEB(aug‐cc‐pvdz)

− ΔEB(6‐31+g(d))] (X = H , G)

(1)

The second term in eq 1 reduces the basis set incompleteness error of smaller 6-31+g(d) basis (due to the use of more complete aug-cc-pvdz basis sets) and hence improves the computed ΔX298B. We have previously used this technique to modify the binding energies quite satisfactorily.27 The results derived from eq 1 in the present manuscript have been further calibrated against the fully optimized data of the benzene··· Mg2+, benzene···Ca2+, and benzene···NH4+ complexes using the same aug-cc-pvdz basis sets (see Section III.a, Table 1). The binding energies of the hydrated complexes were modified using similar approximation of eq 1 from the results of DFT/ B3LYP/6-31++g(d,p) calculations. A single-point coupledcluster calculation with triple excitations (CCSD(T))28 using aug-cc-pvdz basis sets was carried out on the MP2-optimized geometries of the complexes (the DFT/B3LYP geometry in the case of hydrated complexes) to improve accuracy of the computed binding energies (ΔEB). These ΔEB values are both ZPE- and BSSE-corrected. The ZPE corrections of ΔEB at the CCSD(T) level were carried out using the ZPE values from MP2 (DFT/B3LYP values in the case of hydrated complexes) results. The reliability of the computed energies was further gauged against the results from G4MP229 technique. A hybrid variational−perturbational interaction energy decomposition scheme21 was used to compute the contributions of various energy terms in such interactions. The SCF interaction energy is partitioned into first-order electrostatic (HL) (ε(10) el ), Heitler−London exchange (εex ), and higher order HF delocalization (ΔEdel ) energy terms.

II. METHODS OF COMPUTATION Structures of the gas-phase cation−π complexes were initially investigated at the DFT and MP2 levels of theory using 631+g(d) basis sets. The B3LYP functional23−25 was used in the DFT calculations. The hydrated metal cation (Mg2+ and Ca2+)−π complexes were investigated at the DFT/B3LYP level using 6-31++g(d,p) basis sets for such studies. The extra diffuse and polarization functions properly take care of the water interactions in such cases. The geometry of the complexes was fully optimized and the minimum energy structures in the respective cases were confirmed through frequency calculations. The results were used to compute the binding energies (ΔEB), enthalpy (ΔHB298), and Gibbs free energy (ΔGB298) of binding (at 298 K) of the respective complexes with counterpoise (CP) and zero-point energy (ZPE) corrections. The raw binding energies (before CP and ZPE corrections) were calculated from the difference between the total energy of the complexes and the sum of total energies of the corresponding isolated cations and the aromatic moieties. All of the species were considered in their respective equilibrium geometries, and in the case of cations Mg2+ and Ca2+ it was made sure that they were in their respective ground electronic states (1S0). The calculations of thermodynamic properties have been carried out by applying

HF ΔEHF = εel(10) + εexHL + ΔEdel

(2)

The electron correlation effects are taken into account by means of the MP2 theory. The ε(2) MP interaction energy term, which includes the dispersion and correlation contributions to the Hartree−Fock components, is calculated in the supermolecular approach as the difference of MP2 energy corrections of the supermolecule and the monomers (eq 3). (2) (2) εMP = EAB − EA(2) − E B(2)

(3)

The energy terms on the right-hand side of eq 3 represent the difference between the MP2 and Hartree−Fock energies of the supermolecule (AB) and the monomers (A and B). All of the interaction energy terms are calculated consistently in the dimer-centered basis set and are therefore free from the basis-set superposition error due to the full CP correction. The contribution of the multipolar electrostatic interactions in the cation−π complexes has been carried out using the distributed atomic multipolar expansion of the charge distributions of monomers.22 The multipolar expansion technique is based on the numerically equivalent spherical harmonic formulations of 7991

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

Stone et al.30,31 The optimized structures at the MP2 level (optimized DFT/B3LYP structures in the case of hydrated metal ion clusters) were used to generate the electrostatic interaction energy surfaces (using aug-cc-pvdz basis sets). It is to be noted in this connection that the results of energy decomposition calculations, which would be presented here, are not rigorous because there is an element of arbitrariness in any of such decomposition schemes. The geometry optimization and frequency calculations at the DFT/B3LYP and MP2 levels were carried out using Gaussian 09 code.32 The interaction energy decomposition scheme implemented in the GAMESS code33 was used for energy partitioning analyses and computation of the multipolar components of the total electrostatic interaction energies.34 The molecular graphics have been generated using GaussView 05 software.35

(strong and weak) of interactions studied here (Tables 1 and 2). The computed binding energies are dependent on the levels of theory used. The DFT/B3LYP and G4MP2 results agree quite well as the basic geometry optimization procedures are similar (use B3LYP functional). The MP2 and CCSD(T) results generate a lower estimate of binding energies, and they agree in most of the cases (Table 1). We have tested the reliability of the computed ΔEB at different levels against the experimental results of the benzene···Na+ complex. The experimental ΔEB of benzene···Na+ complex is available from two different sources. They are −92.6 ± 5.8 kJ/mol (−22.2 ± 1.4 kcal/mol)36 and −88.3 ± 4.3 kJ/mol (−21.1 ± 1.03 kcal/ mol),37 respectively. We have carried out DFT/B3LYP, MP2, and G4MP2 calculations on this complex at the same basis-set levels of the corresponding Mg2+ and Ca2+ complexes (Table 1), and the results revealed that the respective ΔEB values at the three levels are −23.5, −20.0, and −25.14 kcal/mol. The results show that B3LYP and G4MP2 generate the upper limit, while MP2 computations generate the lower limit of ΔE B (considering the error bar of the experiment). The computed ΔGB298 values (DFT/B3LYP: −14.3 kcal/mol; MP2: −12.3 kcal/mol; G4MP2: −17.1 kcal/mol) also show a similar trend with respect to the estimated ΔGB298 (−15.3 ± 1.96 kcal/mol) from the experimental data.36 The previously computed38 ΔEB values of benzene···Mg2+ (CCSD(T): −109.2 kcal/mol, MP2: −114.4 kcal/mol) and benzene···Ca2+ (CCSD(T): −71.5 kcal/ mol, MP2: −80.6 kcal/mol) complexes were used for further calibration. These computed ΔEB values are comparable to our present results at the CCSD(T) level (Table 1), while our results at the MP2 level are underestimated (Table 1). The computed ΔEB values at the G4MP2 level (Table 1) are again comparable to the previously computed38 MP2 results. The interactions of NH4+ ion with all three of aromatics have been previously studied (both experimentally and theoretically).4,5,10,16 The structure calculations in the present cases have been carried out using the bidentate orientation of NH4+ ion over the aromatic moieties (Figure S1 of the Supporting Information) because the previous studies showed this orientation to be the most favorable one for such cases.5,10,18 In the case of benzene···NH4+ complex, the results at different levels show that MP2 (−17.5 kcal/mol), CCSD(T) (−16.1 kcal/mol), and G4MP2 (−17.2 kcal/mol) computations produce similar ΔEB values (Table1) comparable to previous theoretical results.16,39,40 The computed ΔHB values are also similar (MP2: −17.32 kcal/mol; G4MP2: −17.8 kcal/mol) and comparable to the experiment (ΔHB298: −19.3).41 The DFT/ B3LYP generates a somewhat lower estimate of both ΔEB and ΔHB (−15.1 and −15.7 kcal/mol, respectively, Table 1) values. It was also previously observed that higher basis set systematically underestimates the computed binding energy values of benzene···NH4+ complex at the DFT/B3LYP and CPMD computation levels, and the results of CPMD computations with 100 Ry cutoff match with our present DFT/B3LYP results (r: 3.05; Å; ΔEB: −15.4 kcal/mol).42 The interaction of NH4+ ion with p-methylphenol and 3methylindole is stronger than benzene. The computed binding energies at the MP2 and CCSD(T) levels are similar to the previous results.5 The G4MP2 technique predicts somewhat stronger binding. In all of these computations, the NH4+ ion prefers an orthogonal orientation above the aromatic nucleus. Such orientation of the cation with respect to the aromatic ring leads to the overlap of the π and σ*N−H orbitals. In a geometric

III. RESULTS AND DISCUSSION III.a. Structures and Energetics of the Cation−π Complexes. The results of geometry optimization of various gas-phase cation−π complexes are shown in Table 1 together with their respective binding energy values (ΔEB, ΔHB298, and ΔGB298) and optimal distances (r) between the cation center and the aromatic nucleus. The distance is measured between the cation center and the center of mass (CM) of the aromatic nucleus of the various π-systems. Only the carbon centers of the aromatic part are considered for CM calculations, and the strategy is shown schematically in Figure 1. The optimized

Figure 1. Schematic representation of the cation−π complexes for (a) benzene, (b) p-methylphenol, and (c) 3-methylindole. The green sphere represents the position of the cations (Mg2+, Ca2+, NH4+, Mg(H2O)2+2, and Ca(H2O)2+2) and the red sphere represents the position of the center of mass (CM). The CM is computed using the carbon atoms of the six-member aromatic π-system (as indicated by the gray spheres). The distance r is measured between the CM and the cation center. The optimized structures of the individual complexes are available in Figures S1 and S2 of the Supporting Information.

structures are shown in Figure S1 of the Supporting Information. The results show that the Mg2+ and Ca2+ ions form very strong cation−π complexes and Mg2+ ion−complexes are more strongly bound than their Ca2+ ion counterparts. The interactions are also dependent on the effect of substitution on the aromatic moieties. The cation−π interactions are the strongest with 3-methylindole, and they follow the order 3methylindole > p-methylphenol > benzene in both the classes 7992

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

Table 2. Calculated Interacting Distances (r, Å), Binding Energy (ΔEB), Enthalpy (ΔHB298), and Gibbs Free Energy (ΔGB298) of Binding (all in kcal/mol) of Various Dihydrated Cation−π Complexes at the DFT/B3LYP, CCSD(T), and G4MP2 Levelsa,b dihydrated complexc benzene-Mg

2+

p-methylphenol-Mg2+

3-methylindole-Mg2+

benzene-Ca2+

p-methylphenol-Ca2+

3-methylindole-Ca2+

method

rd

ΔEB

ΔHB298e

ΔGB298e

DFT/B3LYP G4MP2 CCSD(T) DFT/B3LYP G4MP2 CCSD(T) DFT/B3LYP G4MP2 CCSD(T) DFT/B3LYP G4MP2 CCSD(T) DFT/B3LYP G4MP2 CCSD(T) DFT/B3LYP G4MP2 CCSD(T)

2.09 2.07 2.09 2.07 2.04 2.07 2.05 2.02 2.05 2.55 2.46 2.55 2.55 2.44 2.55 2.48 2.38 2.48

−65.97 −71.68 −66.87 −73.75 −80.59 −74.85 −85.83 −92.66 −87.16 −56.30 −55.98 −60.10 −60.55 −64.86 −65.20 −70.93 −73.62 −78.60

−65.80 −72.27

−56.82 −62.39

−73.88 −81.19

−61.25 −68.54

−85.76 −93.25

−74.04 −81.53

−53.31 −56.57

−45.02 −47.63

−60.75 −65.45

−47.45 −51.76

−70.57 −74.21

−59.76 −63.74

Calculated energies (i.e., ΔEB, ΔHB298, and ΔGB298) are after BSSE and zero-point energy (ZPE) corrections. bCalculated energy values at the DFT/B3LYP and CCSD(T) levels are using aug-cc-pvdz basis set of atoms. (See the text for details.) cRefer to Figure S2 in the Supporting Informationfor the structures of the dihydrated cation−π complexes. dSee Table 1 and the text for the definition of r. eSuperscript 298 indicates the temperature (in K) of ΔHB298 and ΔGB298 calculations. a

interacting with the π-aromatic systems. The results of geometry optimization of such complexes are shown in Table 2 together with their computed binding energies and optimized distances (r, as defined in Figure 1). The optimized structures are shown in Figure S2 (Supporting Information). The results show that the hydration of the metal ions substantially reduces the strength of the cation−π interactions. Moreover, the hydrated p-methylphenol···M(H2O)22+ (M: Mg, Ca) complexes show the formation of weak hydrogen bonds between phenolicoxygen and one of the hydrogen atoms of the water of hydration (Figure S2 of the Supporting Information). Thus the interactions are not fully cation−π-type in these two cases. The computed binding energy (ΔEB, ΔHB298, and ΔGB298) trend shows that the CCSD(T) computations slightly improve them over the DFT/B3LYP results. This is understandable, as the CCSD(T) results are obtained using the optimized structures from the DFT/B3LYP analysis. The binding energy improvement is due to the higher correlation effect of CCSD(T) technique. The G4MP2 computations, on the other hand, estimate these values on the higher side, as the calculated r values in various complexes are shorter than those computed at the DFT/B3LYP level. The distance (r) between the metal cation and the aromatic systems is presented in the Tables 1 and 2 and is much shorter and stronger with respect to those in proteins and nucleic acid system because in such cases the metal ions are further coordinated to the other ligands.13,45,46 For example, in the case of 3-methylindole···Ca2+, the computed distance between the metal ion and the aromatic moiety is ∼2.40 Å in the optimized structure (Table 1), while it is 3.37 Å in the metalloprotein SBA,46 with tilted orientation of the metal ions due to coordination with other moieties. In the case of NH 4+ complexes (Tables 1 and 2), although the computed distances and interaction energies are similar to the previous studies,5,10,16,18 they are not pure cation−π types in the real biological systems. For example, our previous investigations on

sense this would lead to dissimilar N−H bond lengths and H− N−H angles. These are depicted in Figure S3 (Supporting Information) for the geometries computed at the DFT/B3LYP and G4MP2 levels. The slight elongation of the N−H length along the direction of the π-system and distortion of the corresponding H−N−H angle (with respect to the other H− N−H angle − not pointing to the π-system) show the overlap of the aromatic π and N−H antibonding orbitals. In all of these NH4+···π interactions, the nitrogen center remains negatively charged. The positively charged hydrogen atoms of NH4+ interact with the π-cloud of the aromatic counterpart. This is considered to be a general characteristic of the electrostatic model of cation−π interactions.43 The hydration of the cationic part influences the binding strength as well as the geometric pattern of the cation−π complexes. The data analysis from protein data bank has revealed that the cation−π interactions are usually stabilized through stacked interactions.13,14 In the present paper our interest is, of course, different. The focus here is to investigate the influence of cation hydration on the gas-phase binding energies of the model cation−π complexes, and we intend to model it using the complexes involving Ca2+ and Mg2+ cations. Previously, similar cases of weak complexes involving NH4+ and N(CH3)4+ cations were investigated,10,18,19 and substantial reduction in the strength of cation−π interactions was observed due to hydration of the cations. The Mg2+ and Ca2+ ions are known to form octahedral Mg(H2O)62+, and Ca(H2O)62+ complexes in their fully hydrated (first hydration shell) forms.44 Although hydrated metal ion clusters are known to stabilize the cation−π interactions through stacking effects, the individual binding energies in such cases are found to be quite low in comparison with the gas-phase complexes.14 Because our present goal is to investigate the effect of hydration on the gas-phase structures of strong cation−π complexes (as discussed here), we restrict our hydrated model systems to Mg(H2O)22+ and Ca(H2O)22+ 7993

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

Table 3. Two-Body Interaction Energy Decomposition Components (aug-cc-pvdz, kcal/mol) of Various Gas-Phase Cation−π Complexes in Their Equilibrium Geometry interaction energy componentsa cation−π complexes

b

2+

benzene-Mg p-methylphenol-Mg2+ 3-methylindole-Mg2+ benzene-Ca2+ p-methylphenol-Ca2+ 3-methylindole-Ca2+ benzene-NH4+ p-methylphenol-NH4+ 3-methylindole-NH4+

εel

(10)

−41.99 −48.46 −54.07 −36.66 −40.99 −47.85 −14.82 −15.86 −19.11

εexHL

ΔEdelHF

ΔEHF

εMP(2)c

28.44 30.87 33.40 26.75 28.44 32.03 12.85 14.70 14.99

−105.87 −113.92 −124.51 −67.76 −73.68 −82.43 −11.61 −14.07 −13.85

−119.42 −131.50 −145.18 −77.66 −86.23 −98.26 −13.58 −15.24 −17.97

6.63 4.60 3.49 −0.36 −2.12 −3.51 −4.44 −5.29 −6.13

a

See the text for the definition of the energy components. bSee Table 1 and Figure S1 in the Supporting Information for the definition of the complexes. cDispersion and correlation contribution as defined in eq 3 of the text.

Table 4. Two-Body Interaction Energy Decomposition Components (aug-cc-pvdz, kcal/mol) of Various Dihydrated Cation−π Complexes in Their Equilibrium Geometry interaction energy componentsa cation−π complexesb 2+

benzene-Mg p-methylphenol-Mg2+ 3-methylindole-Mg2+ benzene-Ca2+ p-methylphenol-Ca2+ 3-methylindole-Ca2+

εel(10)

εexHL

ΔEdelHF

ΔEHF

εMP(2)c

−32.82 −39.95 −44.77 −29.49 −38.61 −42.63

23.26 26.97 28.15 20.36 26.22 25.48

−67.53 −71.84 −79.12 −46.28 −50.32 −52.52

−77.09 −84.83 −95.74 −55.41 −62.71 −69.66

1.39 −0.77 −1.50 −2.02 −4.93 −3.34

a See the text for the definition of the energy components. bSee Table 2 and Figure S2 in the Supporting Inforamtion for the definition of the complexes. cDispersion and correlation contribution as defined in eq 3 of the text.

aromatic systems to the cation (Table 5). This Δq is much less in the NH4+ complexes but not insignificant. A plot of Δq against ΔEB (Figure 2) shows that they are linearly correlated and ΔEB increases with higher charge sharing. Such correlations show that forces involving strong and weak cation−π complexes have similar electrostatic origin, hence the need to analyze the importance of various electrostatic interaction components in such complexes. The contribution of the multipolar components of electrostatic terms in cation−π interactions was investigated in some detail in previous studies.6 The strong electron sharing in the complexes shows that they might have considerable influence in shaping the contribution of total electrostatic interactions between the donors and the acceptors. Figures 3 and 4 represent the natures of various multipolar components of the electrostatic interaction energies in Mg2+, Ca2+, and NH4+ complexes as a function of r (defined in Figure 1). The total electrostatic interaction energies (T-El) together with charge− charge (C−C), dipole-charge (D−C), quadrupole-charge (Q− C), and octopole−charge (O−C) components were found to be important and are presented in the Figures. The influences of the other multipolar components were found to be negligible. The results show that the total interaction energy curves, with few exceptions, are either repulsive or dispersive in nature. The T-El curve of Mg2+ complexes is bound (i.e., having a minimum), while the similar curves for the Ca2+ complexes are dispersive for benzene and 3-methylindole (Figure 3a,c). In this respect, the interactions of p-methylphenol···Ca2+ show a weak minimum (Figure 3c). The nature of the curves representing the Q−C and O−C components of both the complexes as represented in Figure 3b,d shows that they have either repulsive or attractive contributions to the total

the nature of interactions of acetylcholine (Ach) within the acetylcholinesterase (AchE) active cavity showed that the quaternary ammonium interactions of Ach to the Tyr-121 residue of the AchE cavity are of cation−π-type, albeit other interactions are also involved there (due to the surrounding amino acid residues).47 Thus the further interaction energy analyses of the present model systems (which are continued in the next section) should generally be valid in ideal cases, but the results could be extended for the basic understanding of cation−π interactions. III.b. Nature of Cation−π Interactions in Strong and Weak Complexes. Table 3 presents the interaction energy components of the strong (involving Mg2+ and Ca2+) and weak (NH4+) cation−π interactions at the HF level. Similar analyses for the dihydrated metal ion complexes are shown in Table 4. The Tables also contain the contribution of the εMP(2) term, as represented in eq 2. It could be seen that in both gas-phase and hydrated metal ion complexes the contribution of εMP(2) term is not important for too-strong cation−π interactions because they have repulsive contribution in such cases. For the relatively weaker interactions, this term has some contribution to the total interactions (indicated by the computed negative value of εMP(2)) (Tables 3 and 4). In the case of NH4+ complexes (Table 3), of course, this term contributes significantly to the total interaction energies. The main attractive contributors to the total interaction energies are electrostatic (εel(10)) and delocalization (ΔEdelHF) energies. Some part of these attractive interaction energies is offset by the exchange-repulsion term (εexHL), but still the electrostatic contribution plays a significant role in the stabilization of the cation−π complexes. Moreover, the optimized structures of the cation−π complexes of Mg2+ and Ca2+ show substantial charge sharing (Δq) from the 7994

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

(Figure 3a−d) to mostly offset the attractive components to generate overall shape of the T-El curves. The curves representing the electrostatic components of the cation−π complexes of NH4+ ion show an altogether different nature (Figure 4). The D−C components are bound around the equilibrium r for all of the three complexes. The rest of the components around the equilibrium r behave in such a way that the T-El are dispersive (contributing attractive T-El around the equilibrium position, Table 3), although such components sometimes show bound character at larger r (not of much significance here). The curves representing the electrostatic interaction energy components of the dihydrated metal ions (Mg2+ and Ca2+) with various aromatic systems (studied here) are more or less similar to the nonhydrated cases (shown in Figure 3). These curves are shown in Figure S4 of the Supporting Information. Only the C−C electrostatic interaction energy component of the complexes of dihydrated metal ions (both Mg2+ and Ca2+) and p-methylphenol show somewhat deeper potential minima with respect to the corresponding nonhydrated complexes. This is probably due the possibility of hydrogen bond formation in such structures (Figure S2 of the Supporting Information). The differences of the nature of the curves of NH4+···aromatic complexes (Figure 4) with respect to those of the metal ion complexes (Figure 3) could be due to the anisotropic potential expansion of the NH4+ ion. Although dihydrated metal ions are expected to have some anisotropic character in this respect, the curves show they are not strong enough to make too much difference with respect to nonhydrated cases (where the metal ions are strictly isotropic). Thus it is evident from these curves that the multipolar components of electrostatic interactions play a major role in the formation of the cation−π complexes and the nature of interactions in such complexes (as is evident from the shape of the curves in Figures 3 and 4 and Figure S4 of the Supporting Information) depends on the substituents of the aromatic moiety. These substituents are responsible for the charge sharing (or electron flux) from the aromatic moiety to the cation, thus contributing to the stabilizing forces of the cation−π complexes. III.c. Computed Cation−π Vibrational Modes and Their Relation to Electrostatic Interaction Components. The nature of the stabilizing forces in cation−π complexes, as analyzed in the previous section, does not have any direct experimental evidence. The signature of such electrostatic effects could, of course, be predicted from analyses of the proper cation−π vibrational modes. Two different modes could be identified to contain the signature of cation−π interactions in lower frequency range, viz. intramolecular stretching frequency (Sz: up and down motion of the cation combined with the out-of-plane twist of the aromatic moiety) and the outof-plane CH bend (νopCHb) of the π-system.47 The latter one is calculated to be highly infrared (IR)-active and exhibits characteristic blue shift (ΔνopCHb) with respect to the similar mode of the bare π-system. The calculated frequencies involving these two modes (together with their respective IR intensities) for the various cation−π systems are presented in Table 5. A plot of the ΔEB of these complexes against their respective Sz modes (Table 5) show that they are linearly correlated (Figure 5a). The ΔνopCHb does not have such a general correlation with ΔEB values (Figure 5b). The linear correlation of ΔνopCHb with the ΔEB is dependent on the πsystems (but independent of the cations), although it clearly

Table 5. Calculated Charge Shift to the Cation (Δq), Intermolecular Stretching Frequency (Sz, cm−1), and Blue Shift of the Out-of-Plane CH Bend (ΔνopCHb, cm−1) at the G4MP2 Levela,b Δqc

Sz

benzene···Mg p-methylphenol···Mg2+ 3-methylindole···Mg2+ benzene···Ca2+ p-methylphenol···Ca2+ 3-methylindole···Ca2+ benzene···Mg(H2O)22+ p-methylphenol···Mg(H2O)22+ 3-methylindole···Mg(H2O)22+ benzene···Ca(H2O)22+ p-methylphenol···Ca(H2O)22+ 3-methylindole···Ca(H2O)22+ benzene···NH4+

1.02 1.09 1.13 0.63 0.70 0.74 0.74 0.78 0.82 0.40 0.55 0.58 0.14

p-methylphenol···NH4+

0.15

3-methylindole···NH4+

0.16

364 (68.0) 379 (69.0) 390 (72.2) 267 (47.3) 286 (45.6) 298 (66.7) 248 (17.4) 264 (14.1) 273 (16.2) 204 (14.5) 218 (12.1) 228 (15.0) 174 (3.5) 202 (77.0)e 207 (11.0) 222 (63.9)e 221 (22.0) 202 (38.7)e

cation−π complex 2+

ΔνopCHbd 99.0 75.6 85.3 49.0 68.0 77.0 83.1 70.2 79.5 65.9 62.0 63.9 34.8

(76.1) (68.0) (52.0) (81.3) (49.9) (31.7) (68.8) (60.1) (52.3) (53.3) (45.1) (35.8) (27.3)

31.4 (30.8) 34.4 (18.1)

Respective IRI and intensity enhancements for Sz and νopCHb are also included in the Table (within parentheses). bThese values at the G4MP2 level are extracted after the optimization step. Obviously results are from DFT/B3LYP optimized geometries using high-level basis set (as prescribed in the G4MP2 method). cΔq values are calculated from the Mulliken population analysis results. The Mulliken charges considered here could be expected to some error bars, as in many cases elements of different rows are involved. dΔνopCHb and the corresponding intensity increments are computed using the results of the respective π-systems of the cation−π complexes. eThese alternate frequencies also resemble Sz modes of the respective cation−π complexes, but the respective intensities do not reflect the ΔEB complexes nor do they show any linear correlation with respect to the other ΔEb values (as discussed in Figure 2). Thus they are not further discussed in the text. a

Figure 2. Correlation of Δq (see Table 5) with the respective ΔEB of various cation−π complexes.

electrostatic interactions in a significant way around the equilibrium r values, thus revealing the importance of such interactions in strong cation−π complexes. The D−C and (C− C) components of these metal complexes also behave similarly 7995

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

Figure 3. Plots of total electrostatic interaction energies (T-El) and its other multipolar components as a function of r (defined in Figure 1) for the complexes of Mg2+ and Ca2+ ions with benzene (Bz), p-methylphenol (PMP), and 3-methylindole (3MI). Panels a and b represent the curves for the Mg2+ complexes, while panels c and d are for the Ca2+ ion complexes. In the Figures, C−C, D−C, Q−C, and O−C represent the multipolar electrostatic interaction energy components (C−C: charge−charge; D−C: dipole−charge; Q−C: quadrupole−charge; O−C: octopole−charge).

1 1 1 αijVV β VV γ VV i j − i jVk − i jVkVl 2 3! ijk 4! ijkl 1 1 ··· − Ai , jk Vi ∇Vjk··· − Bij , kl VV i j∇Vkl··· 3 6 1 Ci , jklVi ∇2 Vjkl − (4) 15

distinguishes between strong and weak interactions in a specific π-system. The IR intensities (IRI) of the respective modes also show linear changes with respect to the binding energies of the cation−π complexes (Figure 5c,d). The linearity is not always as good as those in the corresponding frequency cases, but the change of IRI for the respective modes could be seen to be more or less regular with respect to the corresponding ΔEB of the complexes. It is well known that IRI is proportional to the square of the magnitude of the dipole derivative, and it has been recently shown that intensity enhancement of a particular mode showing the interactions of constituent components in a weak complex is proportional to the charge flux involving the counterpart molecule/molecules in such a complex.46 In the present complexes, the enhancement of IRI of the νopCHb modes could be related to the measure of intermolecular charge flux, like the previously studied hydrogen bonded complexes.48 These IRI could also be shown to bear the signature of higher order moment contributions through a simple analysis of the of the total energy of the complex. Following the treatment of Qian and Krimm for hydrogen-bonded systems,49 we consider the present cation−π complexes in simple Onsagar reaction field. Such a constant external field is being created by the acceptor atoms of the complexes. The total energy (E) of the molecular system in such a homogeneous electric field (V) could be written as50,51

E = E0 − μi0 Vi −

E0 denotes the energy without any external perturbation. The terms μ0 and α denote dipole moment and polarizability, β and γ are the higher order polarizabilities, Ai,kj is the dipole−dipolequadrupole hyperpolarizability, and Bij,kl and Ci,jkl are dipole− quadrupole and dipole−octopole polarizabilities, respectively. The subscripts i, j, k, and l are the ranks of these tensorial terms. The dipole moment μi generated in the electric field V is defined as μi = −

∂E ∂Vi

(5)

Thus, using eq 4, μi could be written as50,51 1 1 1 βijk VjVk + γijklVjVkVl ··· + Ai , jk ∇Vkj··· 2 6 3 1 1 2 + Bij , kl Vj∇Vkl··· + Ci , jkl∇ Vjkl··· (6) 3 15

μ = μ0 + αijVj +

7996

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

Figure 4. Plots of total electrostatic interaction energies (T-El) and their multipolar components as a function of r (defined in Figure 1) for the complexes of NH4+ ion with benzene (Bz) (a), p-methylphenol (PMP) (b), and 3-methylindole (3MI) (c). See Figure 3 for the definition of the curves representing the cases of C−C, D−C, Q−C, and O−C.

combination of eqs 6−8 reveals that this charge flux during the cation−π complex formation has contributions from the dipole moment, polarizability, and higher order polarizabilites and hence from the higher order moments. The combined properties of the eqs 7 and 8 further prove that the multipolar components of electrostatic interactions play an important role in the stabilization of cation−π complexes because the ∂(∂μ/ ∂Q) term is related to the ΔEB of such complexes. It is to be noted that although Q−C and O−C interactions are present in the cation−π complexes the possibility of detecting such contributions from experiment (except indirect IR intensity evidence) is quite slim. The quadrupolar effects could be observed in some cases from the polar plots of the hyper Rayleigh scattering intensities52 (if the complexes satisfy the necessary conditions for such observation), but the presence of O−C contribution is conceivable from theoretical viewpoint only.

Similar higher order induced moment terms are neglected because they are quite unlikely to contribute in such a weak field. If Q represents the relevant normal coordinates, then the derivative of μl with respect to Q is given by ⎛ ∂μ 0 ⎞ ∂αij 1 ∂βijk = ⎜⎜ i + Vj + VjVk···⎟⎟ ∂Q ∂Q 2 ∂Q ⎝ ∂Q ⎠ ∂μi

⎛ ∂V ⎞ ∂V + ⎜αij i + βijk i Vk + ···⎟ ∂Q ⎝ ∂Q ⎠

(7)

Here Q could be defined as representing the normal coordinates of the cation−π mode. The enhancement of IRI of a particular mode in cation−π complex could be measured from the change in dipole derivative defined as: ⎛ ∂μ ⎞ ⎛ ∂μ ⎞ ⎛ ∂μ ⎞ ∂⎜ ⎟ = ⎜ ⎟ −⎜ ⎟ ⎝ ∂Q ⎠ ⎝ ∂Q ⎠complex ⎝ ∂Q ⎠isolated



(8)

CONCLUSIONS Theoretical investigation has been carried out on several strong and weak cation−π complexes to monitor the nature of electrostatic forces, which contribute to the stabilization of such complexes in the gas phase. Three model π-systems, viz. benzene, p-methylphenol, and 3-methylindole, are chosen to

The enhancement of IRI is directly measurable from the intensities of νopCHb modes of the complex and the isolated species, and they bear linear relation with the ΔEB of the complexes (Figure 5d). It has already been demonstrated that (∂μ/∂Q) is directly related to the charge flux.48 Thus the 7997

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

Figure 5. (a,b) Correlation of ΔEB with the intermolecular stretching (Sz, panel a) and blue-shifted out-of-plane CH bending mode of aromatic moiety (ΔνopCHb, panel b) modes of various cation−π complexes. (b) Red, blue, and green lines respectively, represent the complexes of the cations with benzene, p-methylphenol, and 3-methylindole. (c,d) Correlation of ΔEB with IRI for the respective Sz (panel c) and the intensity enhancement for the respective ΔνopCHb (panel d) modes of the various cation−π complexes. The red line in panel c represents the correlation with all the complexes, while the blue line represents correlation excluding the hydrated metal ions. The red, blue, and green lines in panel d represent correlations for the similar cases in panel b.



bind Mg2+, Ca2+, and NH42+ ions. These investigations also included partially hydrated (dihydrated) metal ions to monitor the effect of hydration on the strong cation−π interactions. The structure of the complexes has been optimized initially at the DFT/B3LYP and MP2 levels with further modification of the final binding energies (ΔEB, ΔHB298, and ΔGB298) using CCSD(T) and G4MP2 techniques. It is observed that the cation shares electronic charges (Δq) from the π-systems and the computed Δq of various complexes is linearly correlated with their ΔEB. This result shows that both the strong and weak cation−π interactions have similar electrostatic origin. The decomposition of interaction energies of the complexes revealed that the electrostatic component of interactions is important in such complexes, and further decomposition of electrostatic interaction contribution into its multipolar components revealed the importance of higher order electrostatic energy components (charge−dipole, charge−quadrupole, and charge−octopole) in shaping the total electrostatic forces in such interactions. The cation−π interaction modes (Sz and ΔνopCHb) of such complexes are analyzed from their computed vibrational spectra, and it has been demonstrated that the intensity enhancement of these modes could be correlated with the multipolar electrostatic behavior of such interactions.

ASSOCIATED CONTENT

S Supporting Information *

Complete reference of Gaussian 09 code and Figures S1−S4. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (J.L.), [email protected] (D.M.), [email protected] (S.R). Notes

The authors declare no competing financial interest. ⊥ K.C.-M. and W.A.E. worked in the project as CRESTsupported ICN summer school students.



ACKNOWLEDGMENTS We acknowledge the support of NSF CREST funding for this work (No. HRD-0833178). This research was further supported by the Extreme Science and Engineering Discovery Environment (XSEDE) by National Science Foundation grantnumber OCI-1053575 and XSEDE award allocation number DMR110088. S.R. acknowledges the financial support 7998

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

(19) Duffy, E. M.; Kowalczyk, P. J.; Jorgensen, W. L. Do Denaturants Interact with Aromatic Hydrocarbons in Water? J. Am. Chem. Soc. 1993, 115, 9271−9275. (20) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford: New York, 1989. (21) Møller, C.; Plesset, M. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618−622. (22) Gora, R. W.; Sokalski, W. A.; Leszczynski, J.; Pett, V. B. The Nature of Interactions in the Ionic Crystal of 3-Pentenenitrile, 2-Nitro5-oxo, Ion (−1) Sodium. J. Phys. Chem. B 2005, 109, 2007−2033. (23) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (24) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58, 1200−1211. (25) Lee, C.; Wang, W.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (26) McQuarrie, D. A.; Simon, D. A. Physical Chemistry A Molecular Approach; University Science Books: Sausalito, CA, 1997. (27) Cato, M. A.; Majumdar, D.; Roszak, S.; Leszczynski, J. Exploring Relative Thermodynamic Stabilities of Formic Acid and Formamide Dimers − Role of Low Frequency Hydrogen-Bond Vibrations. J. Chem. Theor. Comput. 2013, 9, 1016−1026. (28) Pople, J. A.; Head-Gordon, M.; Raghavachari, K. Quadratic Configuration Interaction - A General Technique for Determining Electron Correlation Energies. J. Chem. Phys. 1987, 87, 5968−5975. (29) Curtiss, L. A.; Redfem, P. C.; Raghavachari, K. Gaussian-4Theory using Reduced Order Perturbation Theory. J. Chem. Phys. 2007, 127, 124105-1−124105-8. (30) Stone, A. J. Distributed Multipole Analysis, or How to Describe a Molecular Charge Distribution. Chem. Phys. Lett. 1983, 83, 233−239. (31) Stone, A. J.; Alderton, M. Distributed Multipole Analysis  Methods and Applications. Mol. Phys. 1985, 56, 1047−1064. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B. et al.GAUSSIAN 09, revision C.01; Gaussian, Inc.: Wallingford, CT, 2009. (33) Schmidt, M. S.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (34) Gora, R. W. EDS Package, version 2.8.3; Wroclaw, Poland, 1999−2009. (35) Dennington, R.; Keith, T.; Millam, J. GaussView, version 5; Semichem, Inc.: Shawnee Mission, KS, 2009. (36) Amicangelo, J. C.; Armentrout, P. B. Absolute Binding Energies of Alkali-Metal Cation Complexes with Benzene Determined by Threshold Collision-Induced Dissociation Experiments and Ab Initio Theory. J. Phys. Chem. A 2000, 104, 11420−11432. (37) Armentrout, P. B.; Rodgers, M. T. An Absolute Sodium Cation Affinity Scale: Threshold Collision-Induced Dissociation Experiments and Ab Initio Theory. J. Phys. Chem. A 2000, 104, 2238. (38) Tsuzuki, S.; Uchimaru, T.; Mikami, M. Is the Cation/pi Interaction in Alkaline-Earth- Metal Dication/Benzene Complexes a Covalent Interaction? J. Phys. Chem. A 2003, 107, 10414−10418. (39) Caldwell, J. W.; Kollman, P. A. Cation−π Interactions: Nonadditive Effects Are Critical in Their Accurate Representation. J. Am. Chem. Soc. 1995, 117, 4177−4178. (40) Mavri, J.; Koller, J.; Hadzi, D. Ab Initio and AM1 Calculations on Model Systems of Acctylcholine Binding: Complexes of Tetramethylammonium with Aromatics, Neutral and Ionic Formic Acid. J. Mol. Struct.: THEOCHEM 1993, 283, 305−312. (41) Deakyne, C. A.; Meot-Ner (Mautner), M. Unconventional Ionic Hydrogen Bonds. 2. NH+···π. Complexes of Onium Ions With Olefins and Benzene Derivatives. J. Am. Chem. Soc. 1985, 107, 474−479. (42) Sa, R.; Zhu, W.; Shen, J.; Gong, Z.; Cheng, J.; Chen, K.; Jiang, H. How Does Ammonium Dynamically Interact with Benzene in Aqueous Media? A First Principle Study Using the Car-Parrinello Molecular Dynamics Method. J. Phys. Chem. B 2006, 110, 5094−5098.

by a statutory activity subsidy from Polish Ministry of Science and Technology of Higher Education for the Faculty of Chemistry of Wroclaw University of Technology. We thank the Mississippi Center for Supercomputing Research and Wroclaw Supercomputing and Networking Centers for a generous allotment of computer time.



REFERENCES

(1) Kier, L. B; Aldrich, H. S. A Theoretical Study of Receptor Site Models for Trimethylammonium Group Interaction. J. Theor. Biol. 1974, 46, 529−541. (2) Höltje, H.-D.; Kier, L. B. Nature of Anionic or α-Site of Cholinesterase. J. Pharm. Sci. 1975, 64, 418−420. (3) Dougherty, D. A. Cation-p Interactions in Chemistry and Biology: A New View of Benzene, Phe, Tyr, and Trp. Science 1996, 271, 163−168. (4) Ma, J. C.; Dougherty, D. A. The Cation−π Interactions. Chem. Rev. 1997, 97, 1303−1324. (5) Minuox, H.; Chipot, C. Cation-π Interactions in Proteins: Can Simple Models Provide an Accurate Description? J. Am. Chem. Soc. 1999, 121, 10366−10372. (6) Kim, K. S.; Tarakeswar, P.; Lee, J. Y. Molecular Clusters of πSystems: Theoretical Studies of Structures, Spectra, and Origin of Interaction Energies. Chem. Rev. 2000, 100, 4145−4186. (7) Dougherty, D. A.; Stauffer, D. A. Acetylcholine Binding by a Synthetic Receptor: Implications for Biological Recognition. Science 1990, 250, 1558−1560. (8) Kumpf, R. A.; Dougherty, D. A. A Mechanism for Ion Selectivity in Potassium Channels: Computational Studies of Cation−π Interactions. Science 1993, 261, 1708−1710. (9) Gallivan, J. P.; Dougherty, D. A. Cation−π Interactions in Structural biology. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 9459−9461. (10) Majumdar, D.; Leszczynski, J. Theoretical Modeling of Cation−π interactions in Various Environments: Case Study Using Benzene···Ammonium and Benzene···Tetramethyl-Ammonium Ion Interactions as Model Systems. Comput. Lett. 2007, 3, 257−265. (11) Nicholas, J. B.; Hay, B. P.; Dixon, D. A. Ab Initio Molecular Orbital Study of Cation−π Binding between the Alkali-Metal Cations and Benzene. J. Phys. Chem. A 1999, 103, 1394−1400. (12) Xu, W-z.; Ren, F-d.; Ren, J.; Liu, S-n.; Yue, Y.; Wang., W-l.; Chen, S-s. A UB3LYP and UMP2 Theoretical Investigation on Unusual Cation-p Interaction between the Triplet State HB=BH (3Σg−) and H+, Li+, Na+, Be+, or Mg+. J. Mol. Model. 2010, 16, 615− 627. (13) Zarić, S. D.; Povović, D. M.; Knapp, E.-W. Metal Ligand Aromatic Cation−π Interactions in Metalloproteins: Ligand Coordinated to Metal Interact with Aromatic Residues. Chem.Eur. J. 2000, 6, 3935−3942. (14) Yang, C. M.; Li, X.; Wei, W.; Li, Y.; Duan, Z.; Zheng, J.; Huan, T. Dissecting the General Physicochemical Properties of Noncovalent Interactions Involving Tyrosine Side Chain as a Second-Shell Ligand in Biomolecular Metal-Binding Site Mimetics: An Experimental Study Combining Fluorescence, 13C NMR Spectroscopy and ESI Mass Spectrometry. Chem.Eur. J. 2007, 13, 3120−3130. (15) Lori, M.-I.; Xiuqi, S.; Loren, D. W. Divalent Cations Stabilize Unstacked Conformations of DNA and RNA by Interacting with Base Π Systems. Biochemistry 1998, 37, 17105−17111. (16) Kim, K. S.; Lee, J. Y.; Lee, S. J.; Ha, T.-K.; Kim, D. H. On Binding Forces between Aromatic Ring and Quaternary Ammonium Compound. J. Am. Chem. Soc. 1994, 116, 7399−7400. (17) Gora, R. W.; Barkowiak, W.; Roszak, S.; Leszczynski, J. A New Theoretical Insight into the Nature of Intermolecular Interactions in the Molecular Crystal of Urea. J. Chem. Phys. 2002, 117, 1031−1039. (18) Chipot, C.; Maigret, B.; Pearlman, D. A. Kollman, Molecular Dynamics Potential of Mean Force Calculations: A Study of the Toluene-Ammonium π-Cation Interactions. P. A. J. Am. Chem. Soc. 1996, 118, 2998−3005. 7999

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000

The Journal of Physical Chemistry A

Article

(43) Dougherty, D. A. The Cation−π Interactions. Acc. Chem. Res. 2013, 46, 885−893. (44) Asthagiri, D.; Pratt, L. R.; Paulaitis, M. E.; Rempe, S. B. Hydration Structure and Free Energy of Biomolecularly Specific Aqueous Dications, Including Zn2+ and First Transition Row Metals. J. Am. Chem. Soc. 2004, 126, 1285−1289. (45) Bellsolell, L.; Prieto, J.; Serrano, L.; Coll, M. Magnesium Binding to the Bacterial Chemotaxis Protein CheY Results in Large Conformational Changes Involving its Functional Surface. J. Mol. Biol. 1994, 238, 489−495. (46) Olsen, R. L.; Dessen, A.; Gupta, D.; Sabesan, S.; Sacchettini, J C.; Brewer, C. F. X-ray Crystallographic Studies of Unique CrossLinked Lattices between Four Isomeric Biantennary Oligosaccharides and Soybean Agglutinin. Biochemistry 1997, 36, 15073−15080. (47) Majumdar, D.; Roszak, S.; Leszczynski, J. Probing the Acetylcholinesterase Inhibition of Sarin: A Comparative Interaction Study of the Inhibitor and Acetylcholine with a Model Enzyme Cavity. J. Phys. Chem. B 2006, 110, 13597−13607. (48) Tori, H. Intermolecular Charge Flux as the Origin of Infrared Intensity Enhancement upon Halogen-bond Formation of the Peptide Group. J. Chem. Phys. 2010, 133, 034504-1−034504-10. (49) Qian, W.; Krimm, S. Vibrational Spectroscopy of Hydrogen Bonding: Origin of the Different Behavior of the C-H···O Hydrogen Bond. J. Phys. Chem. A 2002, 106, 6628−6636. (50) Zyss, J.; Ledoux, I. Nonlinear Optics in Multipolar Media: Theory and Experiments. Chem. Rev. 1994, 94, 77−106. (51) McLean, A. D.; Yoshimine, M. Theory of Molecular Polarizabilities. J. Chem. Phys. 1967, 47, 1927−1935. (52) Ray, P. C. Size and Shape Dependent Second Order Nonlinear Optical Properties of Nanomaterials and Their Application in Biological and Chemical Sensing. Chem. Rev. 2010, 110, 5332−5365.

8000

dx.doi.org/10.1021/jp404245q | J. Phys. Chem. A 2013, 117, 7989−8000