Subtle Balance between Interionic Interactions and Ionic Solvation Effect

May 1, 2014 - Shinji Aono,. ‡. Hiroshi Nakano,. § and Takeshi Yamamoto*. ,†. †. Department of Chemistry, Graduate School of Science, Kyoto Univ...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Like-Charge Attraction of Molecular Cations in Water: Subtle Balance between Interionic Interactions and Ionic Solvation Effect Taichi Inagaki,† Shinji Aono,‡ Hiroshi Nakano,§ and Takeshi Yamamoto*,† †

Department of Chemistry, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan Fukui Institute for Fundamental Chemistry, Kyoto University, Kyoto 606-8103, Japan § Department of Molecular Engineering, Graduate School of Engineering, Kyoto University, Kyoto 615-8510, Japan ‡

S Supporting Information *

ABSTRACT: Despite strong electrostatic repulsion, likecharged ions in aqueous solution can effectively attract each other via ion−water interactions. In this paper we investigate such an effective interaction of like-charged ions in water by using the 3D-RISM-SCF method (i.e., electronic structure theory combined with three-dimensional integral equation theory for molecular solvents). Free energy profiles are calculated at the CCSD(T) level for a series of molecular ions including guanidinium (Gdm+), alkyl-substituted ammonium, and aromatic amine cations. Polarizable continuum model (PCM) and mean-field QM/MM free energy calculations are also performed for comparison. The results show that the stability of like-charged ion pairs in aqueous solution is determined by a very subtle balance between interionic interactions (including dispersion and π-stacking interactions) and ionic solvation/hydrophobic effects and that the Gdm+ ion has a rather favorable character for like-charge association among all the cations studied. Furthermore, we investigate the like-charge pairing in Arg-Ala-Arg and Lys-Ala-Lys tripeptides in water and show that the Arg-Arg pair has a contact free-energy minimum of about −6 kcal/mol. This result indicates that arginine pairing observed on protein surfaces and interfaces is stabilized considerably by solvation effects.

1. INTRODUCTION Electrostatic interactions between charged species play important roles in many chemical and biological processes. In the gas phase, like-charged ions repel each other, while oppositely charged ions attract each other. In aqueous solution, the strength of these interactions is weakened significantly by the dielectric shielding of solvent, and more importantly, they can be qualitatively different from what is expected in the gas phase. For example, it is known1 that fluoride ions in aqueous solution tend to form contact ion pairs despite their electrostatic repulsion. The effective attraction between fluoride ions has been attributed to the favorable interaction energy between the ion pair and water molecules that compensates for the ion−ion repulsion energy. Another example of like-charge attraction is that between macroions in solution, which is considered to arise from the accumulation of multivalent counterions at the macroion’s surface.2−5 Like-charge attraction has also received considerable attention regarding the pair formation of positively charged arginine (Arg) residues in proteins. The ionic groups in proteins mainly form positive−negative pairs, about 40% of which involve guanidinium−carboxylate salt bridges. On the other hand, several statistical surveys of the PDB database6−11 demonstrated that a rather large number of protein structures contain like-charged Arg-Arg pairs in close contact. An interesting observation here is that most of the Arg-Arg pairs © 2014 American Chemical Society

are found in the vicinity of protein surfaces, suggesting that solvation contributes significantly to the Arg-Arg pairing. A recent analysis of 70 000 protein structures and complexes also revealed that clusters of four to eight arginine residues (organized as rings, stacks, or strings of stacked arginines) exist at the interfaces of oligomeric proteins.12 These results suggest that the unusual Arg-Arg pairs on protein surfaces may be among the factors responsible for protein−protein interactions. A number of theoretical7,13−22 and experimental23−27 studies have been performed to obtain a deeper understanding of likecharge attraction in water with particular focus on the Arg-Arg pairing. In many of these studies the guanidinium (Gdm+) cation has been used as a simple model of charged arginine residues. Boudon et al.13 showed via Monte Carlo calculations that the Gdm+ ions exhibit an attractive behavior in water with a deep and broad free energy minimum of 9.5 kcal/mol. Similar conclusions were obtained by No et al.14 and Sortens et al.,7 although the well depth of the free energy minimum ranged from ∼3 to 10 kcal/mol depending on the water models used. The free energy profile for the pairing of charged arginine side chains was also obtained from classical molecular dynamics Received: February 3, 2014 Revised: April 30, 2014 Published: May 1, 2014 5499

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

(MD) simulations.15−19 Vondrasek et al.20 performed energy decomposition analysis for the Arg-Arg pair in water and indicated the importance of dispersion and solvent-exclusion effects. More recently, Vazdar et al.21 performed an ab initio MD simulation on the Gdm+ ions in water and observed the clear tendency for like-charge pairing. Experimentally, a variety of techniques have been applied to detect like-charge pairing in water, including neutron diffraction with isotopic substitution,23 capillary electrophoresis,24 dielectric relaxation spectroscopy,25 conductivity measurement,26 and X-ray absorption spectroscopy.27 While several studies23,24,27 support the Gdm+ pairing in water, the results also depend on the experimental techniques, and the conclusion remains controversial.25,26 In many of these studies the amphiphilic character of the Gdm+ ion has been proposed to be crucial for their peculiar association properties. The planarity, aromaticity, and charge distribution of the Gdm+ ion result in weak hydration of Gdm+ faces, with strong hydrogen bonding with water molecules only in the molecular plane. In contrast, it is known that organic cations such as tetramethylammonium do not favorably form a contact ion pair, whereas tetraalkyl ammonium with larger alkyl groups tends to associate in water.28−30 That is, the propensity for likecharge pairing is strongly dependent on specific molecular properties such as charge distribution and hydrophobicity. Therefore, our purpose in this paper is to obtain more physical insight into the like-charge pairing in water by calculating the free energy profiles for relevant molecular cations and making a systematic comparison of energy components that contribute to the pair formation. Specifically, we calculate the free energy profile for eight kinds of molecular cations, i.e., guanidinium, ammonium, methylammonium, tetramethylammonium, tert-butyl ammonium, anilinium, pyridinium, and N-methylpyridinium. These cations have been selected because they are of comparable size to Gdm+, while possessing qualitatively different characters such as charge distribution and hydrophobicity. In addition to these small cations, we consider arginine-alanine-arginine (Arg-Ala-Arg) and lysine-alanine-lysine (Lys-Ala-Lys) tripeptides in water, which is a minimal model for charged side-chain pairs in proteins.19 The free energy profile for these ion pairs is calculated using the three-dimensional reference interaction site model self-consistent field (3D-RISM-SCF) method, which is an electronic structure theory combined with integral equation theory for molecular solvents. By making an energy decomposition analysis we make a systematic comparison of different energy terms that contribute to the pair formation. Furthermore, we perform polarizable continuum model (PCM) and mean-field QM/MM calculations for selected ion pairs for examining the dependence of free energy profiles on solvation models. On the basis of the obtained results, we discuss the competition among interionic interactions and ionic solvation effects that determines the stability of ion pairs in solution.

evaluated using the 3D-RISM method.31,33,34 This method represents the solvent in terms of spatial distribution functions of each solvent site, ρα(R) (α = O,H), defined on a dense grid in a cubic simulation box. The solvent distribution functions on the grid are determined by solving the 3D-RISM integral equation together with the Kovalenko−Hirata closure relation35 (see Supporting Information for details). The 3DRISM method is capable of describing specific solute−solvent interactions such as hydrogen bonds, and it has been applied previously to various systems such as alkali-halide solutions in water36,37 and molecular recognition in solvated proteins.34 The 3D-RISM-SCF method combines electronic structure theory for the solute molecule and the 3D-RISM theory for the solvent to evaluate G(R) in eq 1 in a self-consistent manner.32 This method is conceptually similar to the traditional quantummechanical solvation models such as PCM,38 the difference being that the solvent is represented in terms of spatial distribution functions rather than a dielectric continuum. Note that the free energy G(R) in eq 1 is calculated for fixed solute geometry R, and hence the solute entropic effect is not included in the free energy profile (see Section 3.4 for related discussions). More details on the 3D-RISM-SCF calculations are given in the Supporting Information. In the following sections, we denote the solute electronic energy in eq 1 as Eslt(R), namely Eslt(R) = ⟨Ψ|H0|Ψ⟩

The PMF profile for each cation pair was calculated as follows. First, we located the equilibrium geometry of the contact ion pair by the MP2/3D-RISM method with the 6311G(d,p) basis set. The obtained geometry corresponds to a local minimum on the free energy surface G(R) in eq 1. A reaction path was then obtained by separating the cation pair in a given direction (see Section 3). If a local minimum was not found on the free energy surface (i.e., in the case of a purely repulsive profile), the reaction path was obtained by first performing a constrained geometry optimization with a fixed reaction coordinate and then separating the obtained ion pair in a given direction. The reaction coordinate was chosen as the C−C distance for the Gdm+ pair and the N−N distance for other species. The PMF was then calculated along the reaction path by using the CCSD(T)/3D-RISM method with the 6311++G(2d,2p) basis set. Because the CCSD(T) method is computationally demanding, we calculated the PMFs for tripeptide molecules using the M06-2X/3D-RISM method.39,40 For comparison, we also calculated the PMF using two different solvation models. The first is the continuum solvation model (specifically, the C-PCM41 and SMD42 methods implemented in the GAMESS program).43 The C-PCM method accounts for cavity formation energy in addition to solute−solvent electrostatic energy, whereas the SMD method further includes solute−solvent dispersion and solvent structure terms. The second is a mean-field QM/MM method44−46 which calculates Δμ(R) in eq 1 by explicitly sampling the solvent molecules via molecular dynamics (MD) method. The system consisted of one solute molecule and 860 water molecules described by the SPC/E model.47 All the electrostatic interactions were evaluated with the Ewald method under periodic boundary conditions. We performed the QM calculation for the solute molecule and the MD sampling for water molecules several times until self-consistency was reached. Each MD cycle consisted of 3 ns runs, which is sufficiently long for converging the mean force acting on the

2. COMPUTATIONAL METHODS In this paper we utilized the 3D-RISM-SCF method31−33 for calculating the potential of mean force (PMF) for like-charge association. This method calculates the Helmholtz free energy of a solute molecule with geometry R dissolved in water as follows G(R) = ⟨Ψ|H0|Ψ⟩ + Δμ(R)

(2)

(1)

where the first term is the solute electronic energy and the second term is the solvation free energy. The Δμ(R) is 5500

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

solute atoms. The PMF was then obtained by integrating the mean force acting on the solute atoms along the reaction path. We confirmed that the PMF thus obtained was statistically well converged within an error of 0.5 kcal/mol. The other details of the mean-field QM/MM calculation are the same as described in ref 46. The basis set superposition error (BSSE) was corrected for each value of the reaction coordinate by using the counterpoise (CP) method.48 The CP method reduced the well depth of a free energy minimum by ∼1 kcal/mol for the Gdm+ pair, ∼0.5 kcal/mol for alkyl-substituted amine pairs, and ∼2 kcal/mol for aromatic amine pairs. The CP method was applied to all the PMFs obtained from the 3D-RISM, PCM, and mean-field QM/ MM calculations, except for the tripeptide molecules for which the CP method was not applied for simplicity. All the 3DRISM-SCF calculations were performed using a locally modified version of the GAMESS program.43

3. RESULTS AND DISCUSSION 3.1. Guanidinium and Ammonium Ion Pairs. We first performed the geometry optimization of the Gdm+ pair in water using the MP2/3D-RISM method. Figure 1(a) displays

Figure 2. PMF profiles for the guanidinium (Gdm+) and ammonium (Am+) pairs in water obtained from CCSD(T)/3D-RISM calculation. The x-axis represents the reaction coordinate R, which is chosen as R = R(C−C) and R = R(N−N) for the Gdm+ and Am+ pairs, respectively. The variations in the Helmholtz free energy G(R), solute electronic energy Eslt(R), and solvation free energy Δμ(R) are plotted in red, blue, and green, respectively. Circles represent the point-charge Coulomb potential (i.e., 1/R). All the profiles are anchored at R = 20 Å. Figure 1. Optimized structures of the guanidinium (Gdm+) and ammonium (Am+) ion pairs in water obtained from the MP2/3DRISM calculation: (a) Gdm+ pair and (b) Am+ pair. The reaction path is obtained by separating the ion pairs in the direction shown by the arrows.

to the stacked structure. Indeed, a previous ab initio MD simulation suggests that the conversion from the T-shape to stacked structure occurs very rapidly within 20 ps.21 A similar calculation was performed for the ammonium (Am+) cation pair. Because a free energy minimum was not found for the Am+ pair, the reaction path was determined by first performing a constrained geometry optimization with R(N−N) = 2.5 Å and then separating the cations as shown in Figure 1. The corresponding PMF is displayed in Figure 2(b), which turns out to be purely repulsive in qualitative agreement with previous studies. To obtain more insights, it is useful to decompose the free energy into the solute electronic energy and solvation free energy as follows

the optimized geometry thus obtained, which takes a stacked structure with the NH 2 groups assuming a staggered configuration. This geometrical feature is very similar to that observed in several previous studies.7,13,14,20,21,23 To calculate the free energy profile, a reaction path was defined by separating the ion pairs in the direction shown in Figure 1(a). Here, the intramolecular geometry and the relative orientation of the ion pair were held fixed with only the C−C distance varied. Figure 2 (red line) displays the free energy profile thus obtained at the CCSD(T)/3D-RISM level as a function of R = R(C−C). It should be noted that this figure displays the variation of G(R) in eq 1 as measured from R = R∞, i.e. ΔG(R ) = G(R ) − G(R ∞)

ΔG(R ) = ΔEslt(R ) + ΔΔμ(R )

(4)

where ΔEslt(R) = Eslt(R) − Eslt(R∞) and ΔΔμ(R) = Δμ(R) − Δμ(R∞). The ΔEslt(R) profile in Figure 2 (blue line) represents the variation in the ion−ion interaction energy. For likecharged ions in solution, the ΔEslt(R) is dominated by electrostatic repulsion, and it behaves as ∼1/R at large distances. This is confirmed by comparing the ΔEslt(R) curve with 1/R (depicted by circles in Figure 2). The ΔEslt(R) converges to 1/R much faster for the Am+ pair than for the Gdm+ pair by reflecting the localized charge distribution of the Am+ ion. At short distances, the value of ΔEslt(R) for the Gdm+ pair is much smaller than 1/R because of the flat dendritic

(3)

with R∞ = 20 Å. With this definition the well depth of the free energy minimum is found to be −2.3 kcal/mol, which is roughly comparable to that reported in previous studies.7,14,18,20 In addition to the stacked structure, we also obtained a T-shape structure as a local minimum on the free energy surface. However, the well depth of the T-shape structure was found to be very small (0.2 kcal/mol), indicating that it readily converts 5501

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

shape of Gdm+. That is, the delocalized charge distribution of Gdm+ serves to reduce the interionic electrostatic repulsion, as pointed out in the previous PCM calculation.20 At even shorter distances, the exchange−repulsion interaction becomes operative, and the ΔEslt(R) increases rapidly with decreasing R. The variation of solvation free energy ΔΔμ(R) is also depicted in Figure 2 (green line). In contrast to ΔEslt(R), the ΔΔμ(R) takes on negatively large values by reflecting the attractive interaction between the ion and surrounding waters. It should be noted that ΔΔμ(R) increases (in the negative direction) with decreasing R. This is because the cation pair with a smaller interionic distance exerts a stronger electric field on the surrounding waters, which increases the ion−water interactions (see ref 1 for relevant discussions). According to the continuum model, it can be expected that ΔG(R) behaves as ∼1/(εR) at long distances, where ε is the dielectric constant of water. Because ΔG(R) is the sum of ΔEslt(R) and ΔΔμ(R), it follows that ΔΔμ(R) behaves as −(1 − 1/ε)/R ≃ −1/R at long distances. On the other hand, the molecular nature of solute−solvent interactions becomes important at short distances, and the ΔΔμ(R) deviates from the dielectric behavior. Here, it should be noted that the rather small values of ΔG(R) in Figure 2 result from the cancellation of very large values of ΔEslt(R) and ΔΔμ(R) in opposite signs. This means that the electrostatic repulsion between like-charge ions is largely compensated by the favorable ion−water interactions. The competition between these two factors determines the stability of the contact ion pair in water. The solute electronic energy ΔEslt(R) includes the dispersion interactions between the two ions. To separate the latter contribution, we compare in Figure 3 the PMF profiles

calculated at the Hartree−Fock (HF) and CCSD(T) levels. For the Gdm+ pair, the contact minimum obtained at the CCSD(T) level is much deeper than that at the HF level, whereas the PMF for the Am+ pair is almost identical between the two levels. This result clearly suggests that the dispersion interaction plays a crucial role in stabilizing the contact Gdm+ pair, as noted in previous calculations.20,21 In Table 1 we summarize various properties associated with the contact ion pairs studied in this paper. The Edisp in this table represents the dispersion energy calculated as the difference of ΔG(R) between the CCSD(T) and HF levels. (The other properties in Table 1 will be described later.) Figure 4 depicts the spatial distribution functions (SDFs) of water oxygen atoms around the Gdm+ ions obtained from the 3D-RISM calculation. The SDF for the isolated Gdm+ ion [panel (a)] shows that the solvent density is higher in the molecular plane while much lower above and below the molecular plane. This is because water molecules form hydrogen bonds with the NH2 groups of Gdm+ exclusively in the molecular plane, whereas they tend to be repelled from the hydrophobic faces of Gdm+, thus yielding water-deficient regions above and below the molecular plane. Figure 4(b) depicts the SDF of the contact Gdm+ pair. Because the contact pair assumes a staggered geometry, the SDF exhibits an alternating hydration pattern around the ions. The above features of SDF agree qualitatively with the previous observation that the in-plane hydration is more favorable than the facial hydration for the contact Gdm+ pair.13,19,20,23 While our focus in this paper is on the comparison of PMFs among different species, it is useful to see to what extent the PMF depends on the solvation models used. For this purpose we have calculated the PMF using the PCM and mean-field QM/MM methods (Figure 5). The difference of these solvation models lies in how to evaluate the solvation free energy Δμ in eq 1. Figure 5(a) compares the PMF profiles obtained from the PCM model with and without nonelectrostatic terms. For the Gdm+ pair, all the PCM models yield a free energy minimum of ∼1 kcal/mol, which is somewhat smaller than that obtained with the 3D-RISM method (2.3 kcal/mol). A closer look at this figure also reveals that the cavitation energy slightly increases the well depth for the Gdm+ pair, whereas all the PCM models provide essentially the same results for the Am+ pair. Figure 5(b) presents the PMFs obtained from the mean-field QM/MM method. As seen, the PMF for the Gdm+ pair exhibits a free-energy minimum of about 3 kcal/mol, which is slightly greater than that obtained from the 3D-RISM-SCF method. In addition to Gdm+, we have also calculated the SDF and PMF for the nitrate (NO3−) anion pair in water (Figure S1 in the Supporting Information). The nitrate ion is interesting because of its similarity to Gdm+ regarding the flat dendritic shape and so-called Y-aromaticity. Panels (a) and (b) in this figure depict the SDF of water hydrogen atoms around an isolated nitrate ion and a nitrate ion pair, respectively, obtained from the 3D-RISM calculation. It is seen that the hydration pattern of these ions is qualitatively similar to that of the Gdm+ ions. The solute−solvent hydrogen bonds are formed preferentially in the molecular plane, whereas the solute is much less hydrated above and below the hydrophobic face. Panel (c) displays the PMF for the nitrate ion pair obtained from the CCSD(T)/3D-RISM calculation. The value of ΔG at the contact minimum is −0.8 kcal/mol, which is somewhat smaller than that for the Gdm+ pair at the same level of theory

Figure 3. Comparison of PMF profiles obtained from CCSD(T)/3DRISM and HF/3D-RISM calculations: (a) guanidinium (Gdm+) pair and (b) ammonium (Am+) pair. 5502

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

Table 1. Properties of Various Cation−Cation Pairs in Water Obtained from 3D-RISM-SCF Calculationsa cation +

Gdm Am+ MeAm+ Me4N+ BuAm+ Ani+ Py+ MePy+

Rmin b

3.34 (3.3)c 6.39 (6.5)c 7.93 7.90 5.16 5.48

ΔG

ΔEslt

ΔΔμ

Edispd

ΔEcave

Δ⟨Uuv⟩f

Qg

grouph

−2.3 5.4 1.4 2.0 0.1 −2.8 −0.4 0.9

50.8 73.7 37.8 34.2 28.8 25.0 43.9 39.0

−53.1 −68.3 −36.4 −32.2 −28.7 −27.8 −44.3 −38.1

−4.1 −0.4 −0.9 −1.1 −3.7 −6.9 −5.3 −5.2

−1.7 −1.0 −1.1 0.1 −3.0 −2.4 −2.3 −2.1

−104.7 −136.2 −70.1 −62.6 −52.1 −53.0 −86.2 −75.4

− 1.00 (1.00) 0.76 (0.75) 0.44 (−0.27) 0.56 (0.78) 0.71 (0.85) 0.25 (0.35) 0.20 (0.16)

− NH4 NH3 N NH3 NH3 NH N

a

The short-hand notation for the cations represents guanidinium (Gdm+), ammonium (Am+), methylammonium (MeAm+), tetramethylammonium (Me4N+), tert-butyl ammonium (BuAm+), anilinium (Ani+), pyridinium (Py+), and methyl pyridinium (MePy+). Rmin denotes the interionic distance at contact geometry with the reaction coordinate chosen as R = R(C−C) for the Gdm+ pair and R = R(N−N) otherwise. Energies are in kcal/mol and distances are in Å. bSeparation between Cζ atoms. cDistance at van der Waals contact is used because the PMF profile is purely repulsive. d Difference of ΔG calculated at the CCSD(T)/3D-RISM and HF/3D-RISM levels. eDifference of cavitation energy between Rmin and R∞ = 20 Å obtained from PCM calculation. fDifference of average solute−solvent interaction energy between Rmin and R∞ = 20 Å obtained from 3D-RISM calculation. gPotential derived partial charge of a given group. Numbers in parentheses represent Mulliken charge. hGroup for which the partial charge is calculated.

Figure 4. Spatial distribution functions (SDFs) of water oxygen atoms for solvated Gdm+ ions obtained from the 3D-RISM calculation: (a) isolated Gdm+ ion and (b) contact Gdm+ ion pair. The isosurface is plotted for ρO = 3.5.

(−2.3 kcal/mol). Comparison of the PMFs obtained from the HF and CCSD(T) methods also suggests that the dispersion effect makes a significant contribution to the stabilization of the contact pair, as in the case of the Gdm+ pair. 3.2. Alkyl-Substituted Ammonium Ion Pairs. In the above section we have seen that the Am+ ions do not form a stable contact pair because the electrostatic repulsion between the ions predominates over the effective attraction induced by ion−water interactions. A natural question then arises to what extent the PMF can be made more attractive by making an alkyl substitution for the Am+ ion and thus attaching some degree of amphiphilic character. Here we address this question by considering three kinds of alkyl-substituted ions, namely, methylammonium (MeAm+), tetramethylammonium (Me4N+ or TMA), and tert-butyl ammonium (BuAm+). These ions have been selected because of their comparable size to Gdm+ and also because the TMA is a component of ionic surfactants such as alkyl trimethylammonium ions and some biologically important molecules such as phospholipids. Our focus here is to examine to what extent the amphiphilic character of these ions can stabilize the contact ion pair in water. Figure 6 depicts the optimized geometries of the substituted Am+ pairs obtained from the MP2/3D-RISM calculation. It is seen that the cation pair takes a face-to-face configuration with respect to the alkyl groups while orienting the charged Am+ groups in the opposite direction. This is clearly advantageous in lowering the electrostatic repulsion between the charged Am+ moieties while maximizing the dispersion interaction of the alkyl groups. A reaction path was then obtained by separating the cation pair as shown in Figure 6. The PMF calculated along the reaction path is depicted in Figure 7(a), which shows that

Figure 5. PMF profiles for the guanidinium (Gdm+) and ammonium (Am+) pairs in water obtained from (a) PCM and (b) mean-field QM/ MM calculations. In panel (a), “ES” stands for the PCM calculation with only the electrostatic term included, and “ES + CAV” stands for the PCM calculation with the electrostatic and cavitation terms. The SMD method further includes solute−solvent dispersion and solvent− structure terms. In all the calculations the solute molecule is treated at the CCSD(T) level with BSSE corrections included.

the PMF for the MeAm+ and BuAm+ pairs exhibits a local free energy minimum with an activation barrier of 0.56 and 2.22 kcal/mol, respectively. On the other hand, the PMF for the Me4N+ pair is found to be purely repulsive with no contact minimum. This is consistent with previous studies on aqueous solution of Me4N+Cl− at infinite dilution.28−30 Experimentally, it has been reported that tetraalkyl ammonium ions can exhibit an association tendency only when the alkyl group is more 5503

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

positive charge of these ions is essentially localized on the amine group, and therefore multipole interactions are of minor importance. Table 1 compares the association properties of these cations. The value of ΔEslt for the contact MeAm+, Me4N+, and BuAm+ pairs (29−38 kcal/mol) is much smaller than that for the contact Am+ pair (74 kcal/mol) by reflecting the greater N−N distance. On the other hand, the value of ΔΔμ for the substituted Am+ pairs is smaller in magnitude than that for the Am+ pair. This is because the increased N−N distance makes the electrostatic field of the ion pair weaker, thereby decreasing the strength of ion−solvent interactions. As noted above, the compensation of these two opposing factors determines the stability of the contact ion pair. Table 1 also reveals that the BuAm+ pair benefits substantially from the dispersion and hydrophobic effects of tert-butyl groups. Specifically, Edisp and ΔEcav for the BuAm+ pair amount to −3.7 and −3.0 kcal/mol, respectively, which are much greater than those for the MeAm+ and Me4N+ pairs (at most −1 kcal/mol). These two terms contribute significantly to the relatively large activation barrier of the BuAm+ pair. Nevertheless, we emphasize that the ΔG for the BuAm+ pair (0.1 kcal/mol) is positive, and hence the stability of the BuAm+ pair is smaller than that of the Gdm+ pair with ΔG = −2.3 kcal/mol. This indicates that a rather large alkyl group is necessary for the substituted Am+ ions to exhibit an association propensity comparable to the Gdm+ ions. Figure S2 (in the Supporting Information) displays the solvent SDFs obtained from the 3D-RISM calculation. As expected, the Am+ pair is most strongly hydrated with water oxygen atoms bridging the two ions. The strong hydration of the Am+ pair is reflected in the negatively large value of ΔΔμ (−68.3 kcal/mol). The hydration pattern of the MeAm+ and BuAm+ pairs is very characteristic of amphiphilic molecules; that is, the polar NH3 groups are strongly hydrated by water oxygen atoms, whereas only a small solvent density is found in the vicinity of nonpolar alkyl groups. The solvation pattern becomes qualitatively different for the Me4N+ pair. This is because the four methyl groups prevent a direct interaction of nearby waters with the nitrogen atom on which the positive charge is essentially localized. In addition to the face-to-face configurations shown in Figure 6, we have also calculated the PMF for configurations in which the charged amine groups face each other. The PMF thus obtained (not shown) was found to be purely repulsive, which is reasonable considering the increased electrostatic repulsion between the amine groups and the loss of favorable dispersion and solvent-exclusion effects associated with contact alkyl groups. 3.3. Aromatic Amine Cation Pairs. The stability of the contact Gdm+ pair stems partly from the favorable dispersion/ π-stacking interactions between the Gdm+ monomers. One can thus naturally expect that the introduction of aromatic groups would substantially increase the stability of a cation−cation pair as compared to nonaromatic cations studied above. In this section we thus study three kinds of aromatic amine cation pairs that are of comparable size to Gdm+, namely, anilinium (Ani+), pyridinium (Py+), and methylpyridinium (MePy+) ions. Figure 8 depicts the optimized geometry of these ion pairs in water obtained from the MP2/3D-RISM calculation. This figure shows that the Ani+ pair takes a stacked structure with the NH3 groups orienting in the opposite directions. This configuration is reasonable because it maximizes the dispersion/π-stacking interactions between the benzene rings

Figure 6. Optimized structures of alkyl-substituted ammonium ion pairs in water obtained from MP2/3D-RISM calculation: (a) methylammonium (MeAm + ) pair, (b) tetramethylammonium (Me4N+) pair, (c) tert-butyl ammonium (BuAm+) pair.

Figure 7. PMF profiles for the methylammonium (MeAm+), trimethylammonium (Me4N+ or TMA), and tert-butyl ammonium (BuAm+) pairs in water obtained from CCSD(T)/3D-RISM calculations: (a) ΔG(R), (b) ΔEslt(R) and ΔΔμ(R). In panel (b), circles and squares represent the Coulomb potential between point charges (1/R and −1/R, respectively). All the profiles are anchored at R = 20 Å.

bulky than the ethyl group and that the association of Me4N+ ions is observed only at very high concentration.28 This indicates that the hydrophobicity and dispersion interactions of four methyl groups of Me4N+ are not sufficient for overcoming the residual (or shielded) electrostatic repulsion between the ions. Figure 7(b) displays the energy decomposition of ΔG(R) into solute and solvent contributions. It is seen that the ΔEslt curves for these cations overlap well with the point-charge interaction (i.e., 1/R) except for the contact region (e.g., R < 7.5 Å for the BuAm+ pair). This is to be expected because the 5504

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

that the Ani+ pair has a global free energy minimum with ΔG = −2.8 kcal/mol, which is comparable to that of the Gdm+ pair (ΔG = −2.3 kcal/mol). Several factors contribute to the stability of the contact Ani+ pair. First, the dispersion interaction for the Ani+ pair (Edisp = −6.9 kcal/mol) is the greatest among all the cation pairs studied in this paper (Table 1). For comparison, Edisp for the Gdm+ and BuAm+ pairs is −4.1 and −3.7 kcal/mol, respectively. Second, the distance between the charged NH3 groups is very large (R(N−N) = 7.90 Å), which reduces the electrostatic repulsion significantly as manifested by the small value of ΔEslt (25 kcal/mol). In addition, the hydrophobic effect caused by the stacking of benzene rings further contributes to the pair formation, as seen from the relatively large value of ΔEcav (−2.4 kcal/mol). The combination of all these factors results in the stability of the contact Ani+ pair. Figure 10(a) displays the solvation pattern of

Figure 8. Optimized structures of aromatic amine cation pairs in water obtained from MP2/3D-RISM calculations: (a) anilinium (Ani+) pair, (b) pyridinium (Py+) pair, (c) methylpyridinium (MePy+) pair.

while minimizing the (residual) electrostatic repulsion between the charged amine groups (here it should be recalled that the Am+ ions in water repel each other, as shown in Figure 2). The reaction path was defined by separating the Ani+ pair in the direction perpendicular to the benzene rings. Figure 9 displays the PMF calculated at the CCSD(T)/3D-RISM level. It is seen

Figure 10. Spatial distribution functions (SDFs) of water oxygen atoms around aromatic amine cations obtained from 3D-RISM calculation: (a) anilinium (Ani+) pair, (b) pyridinium (Py+) pair, (c) methylpyridinium (MePy+) pair. The isosurface is plotted for ρO = 3.5. Comparison between panels (b) and (c) shows that N-methylation significantly reduces the solvent density around the cations, which leads to a decreased stability of the MePy+ pair in aqueous solution.

the Ani+ pair in water. As seen, the density of water oxygen is very high around the charged amine groups by reflecting the solute−solvent hydrogen bonds, whereas the solvent density is much lower above and below the hydrophobic faces of benzene rings, as observed for the Gdm+ pair. Because the Ani+ ion is somewhat greater in size than the Gdm+ ion, we next consider a more compact aromatic cation, namely, the pyridinium (Py+) ion (Figure 8). The major difference between the Ani+ and Py+ ions is that the Py+ ion has a significantly delocalized charge distribution over the entire heterocycle. A partial charge analysis shows that only a fraction of total positive charge (+0.25e) resides on the NH moiety, while the remaining charge is delocalized over the other atoms. In this sense, the Py+ ion is of intermediate character between the Ani+ and Gdm+ ions. Figure 8(b) displays the optimized geometry of the Py+ pair, which takes a stacked structure similar to the Ani+ pair. However, the free energy minimum for the Py+ pair is found to be very shallow (ΔG = −0.4 kcal/mol) as compared to the Ani+ pair. The decreased stability of the Py+

Figure 9. PMF profiles for aromatic amine cations in water obtained from CCSD(T)/3D-RISM calculations: (a) anilinium (Ani+) pair, (b) pyridinium (Py+) pair, (c) methylpyridinium (MePy+) pair. 5505

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

pair results from several different factors. First, the N−N distance for the Py+ pair (5.16 Å) is significantly shorter than that for the Ani+ pair (7.9 Å), thus increasing the electrostatic repulsion between the cations. The delocalization of the +1e charge over the pyridine ring also contributes to the increased electrostatic repulsion, as seen from the much greater value of ΔEslt for the Py+ pair (44 kcal) than that for the Ani+ pair (25 kcal). The dispersion and hydrophobic interactions (as measured by Edisp and ΔEcav) are very similar between the Py+ and Ani+ pairs as expected from the similar π-stacking structures. The SDF in Figure 10 shows that the hydration pattern becomes slightly more isotropic by reflecting the more delocalized charge distribution. The overall strength of solvation (as measured by ΔΔμ) is stronger for the Py+ pair because of the proximity of the two nitrogen atoms and the resulting stronger electric field produced by the ions. The sum of all these contributions makes the Py+ pair less stable by 2.4 kcal/mol than the Ani+ pair. In addition to the Py+ pair, we also examined the stability of the methylpyridinium (MePy+) pair in water (Figure 8). At first, we expected that N-methylation of the Py+ ion would make the MePy+ pair more hydrophobic and hence a deeper free energy minimum. In reality, this is not the case, and the contact MePy+ pair is found to be less stable than the Py+ pair. Figure 9(c) displays the PMF profile for the MePy+ pair, which shows that it is the least stable among all the aromatic cations studied in this section. The decreased stability of the MePy+ pair is attributed to a significantly altered hydration pattern illustrated in Figure 10(c). The relatively small density of water oxygen around the nitrogen atoms indicates that the solute− solvent hydrogen bonds are substantially diminished by the methyl groups. Because of this, the solvation of the MePy+ pair is much weaker than that of the Py+ pair, as seen from the comparison of ΔΔμ for these pairs (−38.1 and −44.3 kcal/mol, respectively). On the other hand, the dispersion and hydrophobic interactions as measured by the Edisp and ΔEcav in Table 1 are essentially the same between the Py+ and MePy+ pairs. The combination of all these factors (particularly the significant loss of solute−solvent hydrogen bonds) leads to the lesser stability of the MePy+ pair as compared to the Ani+ and Py+ pairs. Overall, the above results suggest that even a slight modification of molecular cations can substantially alter the propensity for like-charge association. 3.4. Arg-Ala-Arg and Lys-Ala-Lys Tripeptides. As mentioned in the Introduction, a rather large number of protein structures in the PDB database contain argininearginine (Arg-Arg) pairs at close contact (e.g., with the pair distance less than 5 Å). Moreover, most of the Arg-Arg pairs are known to be solvent exposed, suggesting that the aqueous solvation contributes significantly to the stability of the ion pairs. As such, we have also investigated the Arg-Arg pairing in solution by using the 3D-RISM-SCF method. Specifically, the PMF profile was calculated for the Arg-Ala-Arg and Lys-Ala-Lys tripeptides, with the reaction coordinate chosen as the Cζ−Cζ distance of arginines and the N−N distance of lysines, respectively (see Figure 11). Note that the head groups of Arg and Lys are essentially the Gdm+ and Am+ ions, respectively, which are connected by short aliphatic chains and peptide backbones. Those tripeptides were terminated with acetyl and N-methyl residues at the N-terminus and Cterminus, respectively. The reaction path was obtained by optimizing the internal coordinates of the tripeptides for each value of the reaction coordinate at the M06-2X/PCM level.

Figure 11. Geometries of Arg-Ala-Arg and Lys-Ala-Lys in water obtained from M06-2X/PCM calculation. The reaction coordinate is chosen as R = R(C−C) for Arg-Ala-Arg and R = R(N−N) for Lys-AlaLys. The reaction path was obtained by optimizing all the internal coordinates other than R.

The PMF was then calculated along the obtained reaction path by using the M06-2X/3D-RISM method. Figure 12 displays the PMF profiles obtained for these systems. The Arg-Ala-Arg system exhibits a free energy

Figure 12. PMF profiles for the Arg-Ala-Arg and Lys-Ala-Lys tripeptides in water obtained from 3D-RISM-SCF, PCM, and meanfield QM/MM methods. The M06-2X method is used for the solute molecule in all the calculations.

minimum at R(C−C) = 3.2 Å with a well depth of about 6 kcal/mol. This is in contrast to the Lys-Ala-Lys system, which gives a simply repulsive profile with no contact minimum. For more comparison, we have also calculated the PMF using the PCM and mean-field QM/MM methods (Figure 12). Unlike the case of Gdm+ and Am+ pairs, we see that the calculated PMFs are rather insensitive to the solvation models used. The fact that the PMF is attractive for Arg-Ala-Arg and repulsive for Lys-Ala-Lys is qualitatively consistent with several previous studies on relevant systems.17−19 Specifically, Vazdar et al.19 demonstrated from classical MD simulations that only the ArgAla-Arg system forms a stable contact pair at R(C−C) ≃ 4 Å. Note that this distance is somewhat greater than that of the free energy minimum shown in Figure 12 (i.e., 3.2 Å). There are two possible reasons for this difference. First, the present study describes the intramolecular interactions of the peptides using the M06-2X method, whereas the MD simulation in ref 19 utilizes the AMBER ff99SB force field. Second, the solute geometry is frozen in the present calculations, whereas it is statistically sampled in ref 19. This indicates that a more extended form of Arg-Ala-Arg is preferred in the latter study because of the solute entropic effect. Yuzlenko and Lazaridis18 also performed classical MD simulations on the pairing of charged amino acid residues in bulk water and at a lipid−water interface. The obtained PMF for the arginine pair in bulk water 5506

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

supports the idea that like-charge attraction of solvent-exposed arginine groups may be partly responsible for protein−protein interactions.

has a free energy minimum of 4.4 kcal/mol, which is rather similar to the present result (about 6 kcal/mol). It is interesting to note that the well depth for the contact arginine pair shown in Figure 12 is considerably greater than that for the Gdm+ pair shown in Figure 2(a) (2.3 kcal/mol). Apart from the difference in the electronic structure methods used (i.e., M06-2X versus CCSD(T)), it is likely that the intramolecular potential as well as hydrophobicity of the aliphatic chains contribute somewhat to the increased well depth of the arginine pair. We also emphasize that all the PMFs shown above do not include the solute entropic effect because the solute molecule is fixed in the 3D-RISM-SCF calculation. If the entropic effect is taken into account, the free energy of association for the Gdm+ pair may become nearly zero or even positive as a result of greater configurational entropy for the isolated ions. This observation may be related to the fact that the formation of Gdm+ pairs in water remains controversial in the previous experimental studies.23−27 On the other hand, the entropic effect is expected to be rather minor for the Arg-AlaArg system because the Gdm+ headgroups are sterically restricted by the aliphatic chains. Thus, given that the PMF for Arg-Ala-Arg remains attractive after solute entropic correction, the present result suggests that arginine pairs frequently found on protein surfaces are stabilized significantly by solvation effects.

■

ASSOCIATED CONTENT

S Supporting Information *

More computational details on the 3D-RISM method; PMF profiles and solvent distribution functions for the nitrate (NO3−) anion pair in water obtained from CCSD(T)/3DRISM calculation; solvent distribution functions for the alkylsubstituted ammonium pairs in water obtained from 3D-RISM calculation. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +81-75753-4066. Fax: +81-75-753-4000. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by Grant-in-aid for Scientific Research (21350010) and that on Innovative Areas (“Dynamical Ordering of Biomolecular Systems”, 25102002) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan.

4. CONCLUDING REMARKS In this paper we have investigated the like-charge pairing of molecular cations in aqueous solution by using the 3D-RISMSCF method. The free energy profiles for a series of relatively small cations were calculated including guanidinium, alkylsubstituted ammonium, and aromatic amine cations. Our motivation for the present work has been to examine to what extent the Gdm+ ion is so special about the propensity for likecharge association. At first, we expected that the pair of aromatic amine cations can be more stable than the Gdm+ pair because of the greater dispersion, π-stacking, and cavitation (or solvent-exclusion) energies for the former, as verified from Table 1. Nevertheless, it turns out that only the anilinium pair is as stable as the pair with ΔG ranging from −2 to −3 kcal/ mol. This is somewhat surprising because the electrostatic repulsion for the aromatic cation pairs is generally smaller than that for the Gdm+ pair by reflecting the greater separation of charge centers (see ΔEslt in Table 1). On the other hand, the longer distance between the charged amine groups makes the strength of ionic solvation weaker for aromatic amine cations, as seen from ΔΔμ. This is because the electrostatic field produced by the cation pair becomes weaker with decreasing charge density of the solute molecule. This leads to the fact that the effective interionic attraction induced by the ion−solvent interactions (as given by −∂ΔΔμ(R)/∂R) becomes weaker for the aromatic amine cations. The overall balance among all these contributions leads to the fact that only the anilinium pair is as stable as the Gdm+ pair in water. The above observation suggests that the stability of contact ion pairs is determined by the subtle balance between interionic interactions and ionic solvation effects and that the Gdm+ ions possess a rather favorable character for like-charge association in water. Furthermore, we have also seen that charged arginine groups can form a relatively stable pair in water with a binding free energy of about 6 kcal/mol. This result is qualitatively consistent with the experimental fact that Arg-Arg pairs are frequently found on protein surfaces and interfaces and

■

REFERENCES

(1) Zangi, R. Attraction between like-charged monovalent ions. J. Chem. Phys. 2012, 136, 184501−184511. (2) Rouzina, I.; Bloomfield, V. A. Macroion Attraction Due to Electrostatic Correlation between Screening Counterions. 1. Mobile Surface-Adsorbed Ions and Diffuse Ion Cloud. J. Phys. Chem. 1996, 100, 9977−9989. (3) Tanaka, M.; Yu Grosberg, A. Giant charge inversion of a macroion due to multivalent counterions and monovalent coions: Molecular dynamics study. J. Chem. Phys. 2001, 115, 567. (4) Besteman, K.; Zevenbergen, M. A. G.; Lemay, S. G. Charge inversion by multivalent ions: Dependence on dielectric constant and surface-charge density. Phys. Rev. E 2005, 72, 061501. (5) Boroudjerdi, H.; Kim, Y.-W.; Naji, A.; Netz, R.; Schlagberger, X.; Serr, A. Statics and dynamics of strongly charged soft matter. Phys. Rep. 2005, 416, 129−199. (6) Magalhaes, A.; Maigret, B.; Hoflack, J.; Gomes, J.; Scheraga, H. Contribution of unusual Arginine-Arginine short-range interactions to stabilization and recognition in proteins. J. Protein Chem. 1994, 13, 195−215. (7) Soetens, J.-C.; Millot, C.; Chipot, C.; Jansen, G.; Angyan, J. G.; Maigret, B. Effect of Polarizability on the Potential of Mean Force of Two Cations. The Guanidinium-Guanidinium Ion Pair in Water. J. Phys. Chem. B 1997, 101, 10910−10917. (8) Zhang, Z.; Xu, Z.; Yang, Z.; Liu, Y.; Wang, J.; Shao, Q.; Li, S.; Lu, Y.; Zhu, W. The Stabilization Effect of Dielectric Constant and Acidic Amino Acids on Arginine-Arginine (Arg-Arg) Pairings: Database Survey and Computational Studies. J. Phys. Chem. B 2013, 117, 4827− 4835. (9) Lee, D.; Lee, J.; Seok, C. What stabilizes close arginine pairing in proteins? Phys. Chem. Chem. Phys. 2013, 15, 5844−5853. (10) Marsili, S.; Chelli, R.; Schettine, V.; Procacci, P. Thermodynamics of stacking interactions in proteins. Phys. Chem. Chem. Phys. 2008, 10, 2673−2685. (11) Mukherjee, A.; Bagchi, B. Anomalous Orientation-Dependent Effective Pair Interaction among Histidine and Other Amino Acid Residues in Metalloproteins: Breakdown of the Hydropathy Scale Index. Biochemistry 2006, 45, 5129−5139.

5507

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508

The Journal of Physical Chemistry B

Article

(12) Neves, M. A. C.; Yeager, M.; Abagyan, R. Unusual Arginine Formations in Protein Function and Assembly: Rings, Strings, and Stacks. J. Phys. Chem. B 2012, 116, 7006−7013. (13) Boudon, S.; Wipff, G.; Maigret, B. Monte Carlo Simulations on the Like-Charged Guanidinium-Guanidinium Ion Pair in Water. J. Phys. Chem. 1990, 94, 6056−6061. (14) No, K. T.; Nam, K.-Y.; Scheraga, H. A. Stability of Like and Oppositely Charged Organic Ion Pairs in Aqueous Solution. J. Am. Chem. Soc. 1997, 119, 12917−12922. (15) Maksimiak, K.; Rodziewicz-Motowidlo, S.; Czaplewski, C.; Liwo, A.; Scheraga, H. A. Molecular Simulation Study of the Potentials of Mean Force for the Interactions between Models of Like-Charged and between Charged and Nonpolar Amino Acid Side Chains in Water. J. Phys. Chem. B 2003, 107, 13496−13504. (16) Makowski, M.; Liwo, A.; Sobolewski, E.; Scheraga, H. A. Simple Physics-Based Analytical Formulas for the Potentials of Mean Force of the Interaction of Amino-Acid Side Chains in Water. V. Like-Charged Side Chains. J. Phys. Chem. B 2011, 115, 6119−6129. (17) Masunov, A.; Lazaridis, T. Potentials of Mean Force between Ionizable Amino Acid Side Chains in Water. J. Am. Chem. Soc. 2003, 125, 1722−1730. (18) Yuzlenko, O.; Lazaridis, T. Interactions between Ionizable Amino Acid Side Chains at a Lipid Bilayer-Water Interface. J. Phys. Chem. B 2011, 115, 13674−13684. (19) Vazdar, M.; Vymetal, J.; Heyda, J.; Vondrasek, J.; Jungwirth, P. Like-Charge Guanidinium Pairing from Molecular Dynamics and Ab Initio Calculations. J. Phys. Chem. A 2011, 115, 11193−11201. (20) Vondrasek, J.; Mason, P. E.; Heyda, J.; Collins, K. D.; Jungwirth, P. The Molecular Origin of Like-Charge Arginine-Arginine Pairing in Water. J. Phys. Chem. B 2009, 113, 9041−9045. (21) Vazdar, M.; Uhlig, F.; Jungwirth, P. Like-Charge Ion Pairing in Water: An Ab Initio Molecular Dynamics Study of Aqueous Guanidinium Cations. J. Phys. Chem. Lett. 2012, 3, 2021−2024. (22) Villarreal, M.; Montich, G. Energetic and entropic contributions to the interactions between like-charged groups in cationic peptides: A molecular dynamics simulation study. Protein Sci. 2002, 11, 2001− 2009. (23) Mason, P. E.; Neilson, G. W.; Enderby, J. E.; Saboungi, M.-L.; Dempsey, C. E.; MacKerell, A. D.; Brady, J. W. The Structure of Aqueous Guanidinium Chloride Solutions. J. Am. Chem. Soc. 2004, 126, 11462−11470. (24) Kubickova, A.; Krizek, T.; Coufal, P.; Wernersson, E.; Heyda, J.; Jungwirth, P. Guanidinium Cations Pair with Positively Charged Arginine Side Chains in Water. J. Phys. Chem. Lett. 2011, 2, 1387− 1389. (25) Hunger, J.; Niedermayer, S.; Buchner, R.; Hefter, G. Are Nanoscale Ion Aggregates Present in Aqueous Solutions of Guanidinium Salts? J. Phys. Chem. B 2010, 114, 13617−13627. (26) Hunger, J.; Neueder, R.; Buchner, R.; Apelblat, A. A Conductance Study of Guanidinium Chloride, Thiocyanate, Sulfate, and Carbonate in Dilute Aqueous Solutions: Ion-Association and Carbonate Hydrolysis Effects. J. Phys. Chem. B 2013, 117, 615−622. (27) Shih, O.; England, A. H.; Dallinger, G. C.; Smith, J. W.; Duffey, K. C.; Cohen, R. C.; Prendergast, D.; Saykally, R. J. Cation-cation contact pairing in water: Guanidinium. J. Chem. Phys. 2013, 139, 035104−035111. (28) Buckner, J. K.; Jorgensen, W. L. Energetics and Hydration of the Constituent Ion Pairs of Tetramethylammonium Chloride. J. Am. Chem. Soc. 1989, 111, 2507−2516. (29) Shinto, H.; Morisada, S.; Higashitani, K. A Reexamination of Mean Force Potentials for the Constituent Ion Pairs of Tetramethylammonium Chloride in Water. J. Chem. Eng. Jpn. 2004, 37, 1345− 1356. (30) Krienke, H.; Vlachy, V.; Ahn-Ercan, G.; Bako, I. Modeling Tetraalkylammonium Halide Salts in Water: How Hydrophobic and Electrostatic Interactions Shape the Thermodynamic Properties. J. Phys. Chem. B 2009, 113, 4360−4371. (31) Kovalenko, A.; Ten-no, S.; Hirata, F. Solution of threedimensional reference interaction site model and hypernetted chain

equations for simple point charge water by modified method of direct inversion in iterative subspace. J. Comput. Chem. 1999, 20, 928−936. (32) Sato, H.; Kovalenko, A.; Hirata, F. Self-consistent field, ab initio molecular orbital and three-dimensional reference interaction site model study for solvation effect on carbon monoxide in aqueous solution. J. Chem. Phys. 2000, 112, 9463−9468. (33) Hirata, F., Ed. Molecular Theory of Solvation; Kluwer, Dordrecht, 2003. (34) Yoshida, N.; Imai, T.; Phongphanphanee, S.; Kovalenko, A.; Hirata, F. Molecular Recognition in Biomolecules Studied by Statistical-Mechanical Integral-Equation Theory of Liquids. J. Phys. Chem. B 2009, 113, 873−886. (35) Kovalenko, A.; Hirata, F. Self-consistent description of a metalwater interface by the Kohn-Sham density functional theory and the three-dimensional reference interaction site model. J. Chem. Phys. 1999, 110, 10095−10112. (36) Kovalenko, A.; Hirata, F. Potentials of mean force of simple ions in ambient aqueous solution. I. Three-dimensional reference interaction site model approach. J. Chem. Phys. 2000, 112, 10391− 10402. (37) Kovalenko, A.; Hirata, F. Potentials of mean force of simple ions in ambient aqueous solution. II. Solvation structure from the threedimensional reference interaction site model approach, and comparison with simulations. J. Chem. Phys. 2000, 112, 10403−10417. (38) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (39) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (40) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157−167. (41) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model. J. Comput. Chem. 2003, 24, 669−681. (42) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (43) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347−1363. (44) Yamamoto, T. Variational and perturbative formulations of quantum mechanical/molecular mechanical free energy with meanfield embedding and its analytical gradients. J. Chem. Phys. 2008, 129, 244104. (45) Nakano, H.; Yamamoto, T. Variational calculation of quantum mechanical/molecular mechanical free energy with electronic polarization of solvent. J. Chem. Phys. 2012, 136, 134107. (46) Nakano, H.; Yamamoto, T. Accurate and Efficient Treatment of Continuous Solute Charge Density in the Mean-Field QM/MM Free Energy Calculation. J. Chem. Theory Comput. 2013, 9, 188−203. (47) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (48) Boys, S.; Bernardi, F. The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. Mol. Phys. 1970, 19, 553−566.

5508

dx.doi.org/10.1021/jp501212y | J. Phys. Chem. B 2014, 118, 5499−5508