6056
J . Phys. Chem. 1990, 94, 6056-6061
Monte Carlo Smu)atlons on the Lke-Charged Guanidinium-Guanidinlum Ion Pair in Water Stiphane Boudon,* George Wipff, and Bernard Maigret Institut le Bel., Llniuersite Louis Pasteur de Strasbourg, 4, rue Blaise Pascal, Strasbourg 67000. France (Received: October 23, 1989: I n Final Form: January 29, 1990)
Monte Carlo simulations using the statistical perturbation theory were carried out for the guanidinium-guanidinium ion pair in water at 25 OC and 1 atm using the isothermal-isobaric (NPT) ensemble to obtain interionic potentials of mean force and to evaluate "effective" interactions in water at short distances. The system is found to be attractive with a deep and large minimum at a C-.C distance of approximately 3.30 A. Our results address fundamental questions about interionic interactions in recognition processes involving like-charged binding sites. Contrary to the simple electrostatic picture, they may attract locally due to related induced solvent effects.
Electrostatic interactions are important determinants for protein function and structure1 and for a wide variety of cellular and biochemical processes.* The long-range character of electrostatic interactions makes them unique in their involvement in major biological phenomena. Similar electrostatic interactions also drive the complexation and recognition of ionic substrates by natural and synthetic macrocyclic receptor^.^ Stability constants for a series of complexes demonstrate the primary role of electrostatics in addition to size and shape ~omplementarity.~Converging like-charged binding sites of molecular receptors, although destabilizing the receptor itself, lead to enhanced stabilities for the complexes. In these compounds, however, most of the energy cost for such an organization has been already paid in the course of the synthesis. Accumulation of like charges at the surface or in active sites of proteins may be observed. For instance, in the active site of ribonuclease there is a large number of positively charged groups which may be involved in binding the ~ u b s t r a t e . ~ In the BPTI/BPT6 complex like-charged groups belonging to both partners are found in the binding site region in close p r ~ x i m i t y . ~ Such proximities are surprising given the extremely high association constant of BPTI with BPT (>1013 M-I at neutral pH). This suggests that a simple electrostatic model in vacuo does not account for structure or stability and that environment and solvation effects play a major role. In fact, structural* and theoreticalgJOstudies demonstrate that water networks stabilize clusters of like charges. Tabushi et al. calculated that four water molecules located in the vicinity of Arg+ and Lys+ significantly reduced the electrostatic repulsion^.^ Molecular dynamics simulations on the active site of ribonuclease A without anions show stabilization of positively charged groups (Arg-39, Lys-41, Lys-7, His-1 19) by bridging water molecules not observed in crystallographic s t ~ d i e s . Similar ~ simulations on the active site of lysozyme beautifully showed the formation of a like-charged ion pair between Arg-61 and Arg-73 mediated by bridges of water molec u l e ~ .In~ ~these complex systems it is not clear whether such proximities of like-charged groups are due to structural constraints, ( I ) Perutz, M. F. Science 1978, 201, 1187-1 191. (2) Warshel, A . Acc. Chem. Res. 1981, 14, 284-290. (3) Lehn, J. M. Angew. Chem., Int. Ed. Engl. 1988, 27, 89-112, and references therein. (4) (a) Schneider, H. J.; Theis, 1. Angew. Chem., Int. Ed. Engl. 1989, 28, 753-754. (b) Izatt, R. M.; Bradshaw, J. S.; Nielsen, S. A.; Lamb, J. D.; Christensen, J. J. Chem. Reo. 1985, 85, 271-339. ( 5 ) Matthew, J. B.; Richards, F. M. Biochemistry 1982,21,4989-4993. (6) Bovine pancreatic trypsin inhibitor/bovine pancreatic trypsin. (7) Tabushi, 1.; Kiyosuke, Y.; Yamamura, K. J. Am. Chem. SOC.1981, 103, 5255-5257, and references therein. (8) Baker, E. N.; Hubbard R. E. Prog. Biophys. Mol. Biof. 1984, 44, 97-1 79. (9) Briinger, A. T.; Brooks 111, C. L.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 1985. 82. 8458-8462. (10) Br&ks-iIl, C. L.; Karplus, M. J . Mol. Biol. 1989, 208, 159-181
0022-3654/90/2094-6056$02.50/0
TABLE I: Full Geometry Optimization at the Hartree-Foek Level of the Guanidinium Ion Using Various Basis Sets"
basis set
energy
C-N
N-H
H-N-H
STO-3G 3-21G 6-31G(d)
-201.960712948 -203.40I 115429 -204.521 536401
1.3551 1.3245 1.3216
1.0201 0.9995 0.9965
120.90 121.77 121.61
"Energy in hartrees, distances in an stroms, angle in degrees. Experimental C-N distance 1.323-1.335 .!
result from local field effects due to the solutes, or result from solvent effects. In order to gain insight into the effective interaction between like charges at short distance in water, we decided to simulate the guanidinium-guanidinium (Gu+/Gu+) ion pair in water and calculate the free energy profile as a function of the separation between the ions (potential of mean force). Statistical mechanical theories as well as molecular dynamics calculations have been able to describe ion pairing in water involving CI-/CI-. The computation of the potential of mean force (pmf) suggests that the Cl-/Cl- pair could be weakly bound via the solvent and that the stability of this 'solvent-bridged" ion pair could be related to water molecules that simultaneously hydrogen-bond to both anions."J* In contrast, Monte Carlo simulations showed that the NMe4+/NMe4+ pair is repulsive12since the pmf shows no clear minimum. Another study using the reference hypernetted chain (RHNC) approximation applied to aqueous ionic solution^'^ showed attractive wells in the pmf for negative as well as positive like-charged ion pairs. Compared to CI- and NMe4+ the Gu+ ion is aspherical, and we felt it interesting to characterize its hydration pattern. In addition, Gu+ is of biological interest as a model for arginine side chains. Computational Procedure The procedure is comparable to that employed by Jorgensen to study the hydration of aqueous tert-butyl ch10ride.I~ A series of ab initio molecular orbital calculations on Gu+/H20 and Gu+/Gu+ interactions were performed to calibrate the force field used in a Monte Carlo (MC) simulation. To obtain the potentials of mean force, statistical perturbation theory was used. The geometrical parameters of the guanidinium ion have been determined from ab initio SCF-LCAO-MO calculations using the ~~
(1 I ) (a) Dang, L. X.;Pettitt, B. M. J. Chem. Phys. 1987,86 6560-6561. (b) Dang, L. X.;Pettitt, B. M. J. Am. Chem. SOC.1987, 109, 5531-5534. (12) Buckner, K.; Jorgensen, W. L. J . Am. Chem. SOC.1989, 1 1 1 , 2507-25 16. (13) (a) Patey, G. N.; Carnie, S. L. J. Chem. Phys. 1983,78, 5183-5190. (b) Pettitt, B. M.; Karplus, M. J. Chem. Phys. 1985,83, 781-789. (c) Pettitt, B. M.; Rmky, P. J. J. Chem. P h p . 1986,84, 5836-5844. (d) Pettitt, B. M.; Rossky, P. J. Isr. J. Chem. 1986, 27, 156-162. (14) Jorgensen, W. L.; Buckner, J. K.; Huston, S. E.; Rossky, P. J. J . Am. Chem. SOC.1987, 109, 1891-1899.
0 1990 American Chemical Society
Guanidinium-Guanidinium Ion Pair in Water
The Journal of Physical Chemistry, Vol. 94, No. 15, 1990 6051
TABLE 11: Comparison of the ab Initio and Empirically Computed Guanidinium-Water Interaction Energies (See Figure 1)“ structure a b initio empirical structure a b initio emoirical I 14.1 12.9 IV 10.2 10.0 I1 18.1 15.6 V 7.7 6.3 111 10.3 9.2 ~~
a
N k
1 !
NY
NY
0
H’
dco=3.6A
‘H
L
NHZ
0
H’
dco=3 4h
H‘
Energy in kcal/mol.
TABLE 111: Geometric and Energetic Parameters model water
TIP4P
site 0 H
M guanidinium
C N
H OAii
geometry roH = 0.9572
A
aHOH 104.52’ roM = 0.150 A rCN = 1.321 A aNCN = 120.0’ rNH = 1.000 A
4i +O.OOO +0.520 -1.040 +0.640 -0.800 +0.460
Aii
Gii
600003 0 0
610
3367 944300 0
NY
0 0 26 800
NY
H r
0
H
and Cii are expressed in kcal Ai2/mol and kcal A6/mol.
extended 6-31G(d) basis set which includes a set of d-type polarization functions on all atoms other than hydrogen.15 Neither the correlation energy nor a basis set superposition error estimate was calculated. The optimized parameters are reported in Table I an can be compared t o values obtained with smaller basis setsI6 or experimental values.” The guanidinium ion was kept rigid in all the simulations, since this ion belongs to a class of compounds often referred as “Y aromatics” in which large rotational barriers (IO kcal/mol) have to be overcome to interconvert two stable D3,, forms via a CZUform. The ion-water and water-water interactions are described by Coulomb and Lennard-Jones interactions between sites centered on nuclei. The intermolecular interaction energy is given by the sum of interactions between sites i and j of the two molecules a and b (eq 1) with combining rules for the Lennard-Jones interactions such that A , = (AijAjj)Il2and C, = (CiiC,j)l/2.
The energy parameters for the guanidinium ion were arrived at through a b initio SCF-LCAO-MO calculations with the 63 IG(d) basis set for five low-energy forms of complexes between water and ions (Figure I ) . The charge and Lennard-Jones parameters were chosen to reproduce the resultant optimal ab initio geometries of various G u + / H 2 0 dimers. The calculated interaction energies with water (Table 11) are deliberately designed to be slightly less than the interaction energies determined by 6-31G(d) quantum mechanical calculations since the latter are usually somewhat greater (more exothermic) than the experimental data for closely related systems.’* The large polarity of the C-N and N-H bonds in this model is noticeable (Table 111), and 36% of the positive charge is spread over the -NH2 groups. The TIP4P model of water was used in conjunction with the above-described parameters (Table 111). The intermolecular potential function is pairwise additive and includes neither higher order effects due to three-body interactions nor polarization effects. Nevertheless, pair potentials can yield thermodynamical results (hydration enthalpies, relative hydration free energies, partial molar volumes, etc.) comparable with experiment when dealing with various ions in ~ a t e r . ~ ~ , ~ ~ (15) Hehre, W. J.; Radom, L.; Schleyer, P.v.R.; Pople, J. A. Ab-Initio Molecular Orbital Theory; Wiley: 1986. (16) (a) Kollman, P.; McKelvey, J.; Grund, P. J . Am. Chem. SOC.1975, 97, 1640-1645. (b) Capitani, J. F.; Pedersen, L. Chem. Phys. Lett. 1978,54, 547-550. (c) Sapse, A. M.; Massa, L. J. J . Org. Chem. 1980, 45, 719-721. (d) Sapse, A. M.; Snyder, G.;Santoro, A. V. J . Phys. Chem. 1981, 85, 662-665. (e) William, M. L.; Gready, J. E. J . Compuf. Chem. 1989, 10,
35-54. (17) (a) Adams, J. M.; Small, R. W. H. Acra Crysfallogr.1974, B30, 2191-2193. (b) Curtis, R. M.; Pasternak, R. A. Acra Crystallogr. 1955,8, 675-681. ( 18) For example, for ammonium-water or methylammonium-water complexes 6-31G(d) calculations overestimate the interaction energies by about 1-2.5 kcal/mol or 10% of the total. See ref 19.
H
J-“. NY
0/
NY
dc0=3 O A
H ‘
(U)
Figure 1. Five geometries optimized via a b initio calculations. Structures I and 111 are in-plane and structures 11, IV, and V are perpendicular. The N-Ow,,, and C a w a mdistances were optimized at the 6-3 1(d) level. The previously determined geometries of the ion and water were kept fixed.
The system consisted of 212 water molecules with a single guanidinium ion or of 310 water molecules plus the ion pair each in a box with eriodic boundary conditions of 19 X 19 X 19 and respectively. The Monte Carlo Metropolis (MC) 18 X 18 X 28 procedure was used2’ for these isobaric-isothermal (NPT) ensembles22modified by preferential ampl ling.^^,^^ The solutesolvent statistics were enhanced by trying to move the solute every 60 configurations. The volume changes were attempted on every 2500 configurations and involved scaling all the intermolecular separations. Configurations were generated by moving a water molecule randomly in all three Cartesian directions and by rotating it about a randomly chosen axis25or by moving one of the two ions while keeping the C-C distance fixed. The perturbation is done along the C-C distance of the ions running from 2.20 to 6.10 8, by steps of 0.15 A. Like in a comparable pmf calculation26relative orientations of the solutes have been allowed to vary at each separation. Uniform sampling of all rotational degrees of freedom for each separation could not be achieved within reasonable computer time. We however preferred a “relaxed” to a fully ‘constrained” pmf calculation in
i3,
(19) Jorgensen, W. L.; Gao, J. J . Phys. Chem. 1986, 90, 2174-2182. (20) (a) Singh, U.C.; Brown, F. K.; Bash, P. A,; Kollman, P.A. J . Am. Chem. Soc. 1987,109,1607-1614. (b) Rao, B. G.; Singh, U.C. J . Am. Chem. SOC.1989, 111, 3125-3133.
(21) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M.N.; Teller, A. H. J . Chem. Phys. 1953, 21, 1087-1092. (22) McDonald, I. R. Mol. Phys. 1972, 23, 41-58. (23) (a) Owicki, J. C.; Scheraga, H. A. Chem. Phys. Lerr. 1979, 47, 600-601. (b) Bigot, B.; Jorgensen, W. L. J . Chem. Phys. 1981, 74, 1944-1 952. (24) A I/($ + c) wei hting factor was chosen where c = 180 A2for the small system and c = 5 0 0 5 1 ~for the large system in order to have an adequate
sampling. (25) The ranges for these motions were chosen to provide acceptance probabilities of ca. 40% for new configurations and were 0.15 A and 15’ for water monomers and 0.06 A and 6O for the solutes. (26) Jorgensen, W.L. J. Am. Chem. SOC.1989, I l l , 3770-3771.
6058 The Journal o j Physical Chemistry, Vol. 94, No. 15, 1990 TABLE IV: Tbermodynnmic Properties ET -22 14.2 E& -2144.8 E,, -1 32.71 AE%i -69.3 E, -208 1.5 V 6401.7
VO AVwI
Boudon et al. Energy Pair Distribution
6360.5 41.2
5-
4 -
'Energies in kcal/mol and volumes in A'. The reference energy Esse and volume VO have been calculated on an equilibrated water box by averaging the thermodynamic parameters over 2 X IO6 configurations.
order to avoid being trapped artificially in local minima and allow for smoother perturbations. Indeed, we experienced on an in-plane constrained (C, symmetry) acetate-guanidinium pmf calculation27 that perturbations of the solvent-separated pair caused artifically large fluctuations in calculated free energy. This was related to the prevented solventsolute concomitant reorganization. Ideally, one should compute the potential of mean force over all degrees of freedom of the system. Double wide sampling was used so that 0.3 8, could be covered in one MC simulation, and a total of 13 simulations with nonoverlapping windows were realized. At each C-C separation the MC runs involved equilibration over IO6 configurations and averaging over 2.4 X IO6 configurations obtained from three subaverages over 8 X IO5 configuration batches. At the final step the solvation free energy is still decreasing (Table V); nevertheless, this point will be taken as reference and set to 0 kcal/mol. The atom-atom interactions of the system were calculated with a 8.5-A cutoff as used previously by Jorgensen for studies on CI-/CI- ion pairing.I4 It has been pointed out that the necessary truncation of the potential energy function could introduce artifacts into the simulation,28because the solvent might be disrupted in the cutoff region. However, neglecting long-range electrostatic interacitons has less influence on relative free energies when the charge of the solute is kept constant rather than mutated.29 The results of perturbation studies on relative free energy of hydration for solutes of constant charges with a similar c u t o f p are in overall agreement with experimental values. Here the charge of the system remains constant, and we expect no significant change in the shape of the calculated pmf using a longer cutoff. Additional testing of the effect of this cutoff on the pmf would require intensive calculations which are beyond our actual computational resources. The free energy profile for the Gu+/Gu+ ion pair is computed in the context of statistical perturbation t h e ~ r y . ~This ' approach is based on eq 2 in which the free energy difference between two G I - GO= -kJ In(exp(-(H, - H o ) / k s T ) ) o (2) states 0 and I , 0 being the reference state, is given as a function of a cumulant expansion in the energy difference between the two systems. At each step the change in free energy of hydration was evaluated as the difference in ion-water interactions. The total energy change is computed by adding the cumulative difference in free energy of hydration to the mean gas-phase interaction energy between the two ions. Results and Discussion
Hydration of a Single Guanidinium ion. Since the experimental thermodynamic data relative to the hydration of the guanidinium ion do not appear to be available, the solvation will here be characterized only by using a MC simulation of the ion in water. The total potential energy of the solution can be written as the sum of the solute-solvent energy and the solvent-solvent energy: ET = CEs,x + E X E S , or ET = E5x + E , (3) i
i pi
The index i runs over all water molecules of the box. The energy (27) Boudon, S.;Wipff, G. Unpublished results. (28) Madura, J. D.; Pettitt, B. M. Chem. Phys. Lerr. 1988,150, 105-108. (29) Pearlman, D. A.; Kollman, P. A. In Computer Simulation of Biomolecular Sysrems; ESCOM: Leiden, 1989. (30) Beveridge, D. L.; DiCapua, F. M. Annu. Reu. Biophys. Chem. 1989, 18. 431-492, and references therein. (31) Zwanzig, R. W. J . Chem. Phys. 1954, 22, 1420-1426.
3-
-10
-20
.30
0
10
Interaction Energy , kcal/mol Figure 2. Energy pair distribution for the Gu+ ion in water. Units for the ordinate are number of solvent molecules per kcal/mol.
+ I ; ; ; ; =
;
;
;
;
5
6
7
. l
0
1
2
3
4
Distance,
Figure 3. N-O,,,,,
A
radial distribution function of Gu+.
change on transferring the solute from the gas phase to the solution (partial molar internal energy) can be expressed by
+
AESl = Esx E5s - E5S0 (4) where ESSOis the energy for the pure solvent (212 water molecules) and E% - ESSothe solvent reorganization energy. The solvation enthalpy may then be written AHSO, = AE%l+ PAVSi - R T (5) where AVSl is the partial molar volume and R T the gas-phase contribution for the solute in the ideal-gas approximation. The calculated partial molar volume has to be considered cautiously since large standard deviations (about 50 A3) are observed on computed volumes similar to those already reported. Nevertheless, the magnitude of the computed values (Table IV) seems to be acceptable, and the solvent density remained at the experimental value of 1.0 g/cm3. The results compare well to similar s t ~ d i e son l ~ammonium ~~~ ions. On a crude hydrophobicity scale the Gu+ ion is in between CH3NH4+and (CH3)4N+. The solvation enthalpy AHsl for the guanidinium ion (-66.8 kcal/mol) is lower than for CH3NH3+33*34 (-75, -79 kcal/mol) but higher than for (CH3)4N+33*35 (-53, -41 kcal/mol). The distribution of individual solute-solvent interaction energies is shown in Figure 2. The energy band that stretches from the minima near -10 kcal/mol to around -17 kcal/mol results from the interactions of Gu+ with the nearest water molecules, which give also the highest interactions and form mostly in-plane arrangements (of type I t o IV in Figure 2). The energy of fa(32) Alagona, G.; Ghio, C.; Kollman, P. J . Am. Chem. Soc. 1986, 108, 185-1 91. (33) Aue, D. H.; Webb, H. M.; Bowers, M. T. J . Am. Chem. Soc. 1976, 98, 318-329. (34) Amett, E. M.; Jones, F. M.; Taagpera, M.; Henderson, W. G.; Beauchamp, J. L.; Holtz, D.; Taft, R. W. J . Am. Chem. Soc. 1972, 94, 4724-4726. (35) Boyd, R . H . J . Chem. Phys. 1969,51, 1470-1474.
The Journal of Physical Chemistry, Vol. 94, No. 15, 1990 6059
Guanidinium-Guanidinium Ion Pair in Water
3-
."
-." E
Ol
l
$ 1
.= . .-........................ .........
. . . .... .'.
.--0 U
2-.
. .
.........
-."m
1
-.
0
-.
U
............
E
.
z
-1:
I
1
2
3
I
1
I
I
4
5
6
7
I
I
A Figure 4. H-Owe,,, radial distribution function of Gu'.
,
0
8
1
i
i
m
'p
E 0-.
Figure 6. N-H,,,,,
.. .
......
..................
;
i
A radial distribution function of Gu'.
.. -."
;
i Distance,
Distance,
P e,
...... ..............
..........
z
0
T
L
I
.
.'.=.......... ...... ....
I
9.
9
?
i, -
I
0
1
;
;
:
4
4
;
- 1
8
A radial distribution function of Gu'. Distance,
Figure 5. C-O,,,,
cial-coordinated water (form V in Figure 1) being significantly smaller might contribute to the population of the main band. Thc solute-solvent structure can be further characterized by radial distribution functions (rdf s) and running coordination numbers for the three different atoms of the guanidinium ion. The rdf is defined as the probability of finding an atom of type y at a distance r from an atom of type x normalized for the bulk density p of atoms y and for the volume element 47rr2 dr. This can be written as g J r ) = dN/(47rp$ dr) = (N(r,r+Ar)) /(47rpr2Ar) (6)
In-plane hydrogen bonding is revealed by the sharp peaks of the N-Ow,,,, or H-Owate,centered at 2.75 and 1.75 A, respectively (Figures 3 and 4). Integration up to the first minima at 3.15 and 2.45 A gives 2.1 water molecules for the first hydration shell around each of the -NH2 groups. The C-Owate,rdf (Figure 5) shows a broad peak centered at 3.75 A, a result consistent with the energy distribution for in-plane arrangements (Figure 2). Much less structure is observed in the second shell of water around the carbon atom than around -NH2. The N-Owate,rdf clearly shows a second and third peak centered at 3.85 and 4.85 A. Integration yields 8.0 and 13.1 water molecules in these successive water shells. Comparison of the first peaks in the N-O,,,,, and N-H,,,,, rdfs (Figure 6) respectively centered at 2.75 and 3.35 A reveals the orientation of water molecules. They are coordinated to the -NH2 group via -NH2--OH2 in-plane interactions rather than perpendicular -H2N-H-OH interactions. The relatively broad and short peaks show that the related water molecules have a large orientational freedom. Similarly, the C-OWate,rdf and the C-H,,, rdf first peaks respectively centered at 3.75 and 4.1 5 A (Figures 5 and 7) reveal that the hydrogen atoms are pointing away from the cationic center. The mean water positions in crystallized proteins36and our analysis of the solvation show that the binding of water to guanidinium is preferred within the plane of the group. No water molecules were found in close contact (36) Thanki, N.; Thornton, J. M.; Goodfellow, J. M. J . Mol. Biol. 1988, 202,637-657.
Figure 7. C-H,,,,,
J
Distance, A radial distribution function of Gu'.
I
A
12 8-
P
4-
WC 0-
4-
8-
C...C Distance, A
Figure 8. Calculated potentials of mean force for two Gu' ions in water along the C-.C distance coordinate (kcal/mol). Insets show the structure of the ion pair without the solvent at 3.45 and 5.95 A.
with the central carbon atom in the crystalline state. Pmf for the Gu+/Gu+ Ion Pair in Water. The calculated pmfs are displayed in Figure 8, and the related numerical values are given in Table V. The fluctuations (a) in the free energy changes averaged 0.29 kcal/mol per step. The end point of the pmf curve at a 6.10-A distance is taken as reference of the perturbation energy. Although the uncertainities in the results are difficult to establish unequivocally, they show a broad minimum for the contact ion pair occurring at a carbon-carbon distance of 3.30 A. Around 5.50 A we observe a change in the slope of the pmf curve, whether this is an artifact due to approaching the cutoff region or reflecting the real dissociation cannot be assessed without additional calculations. The difference in energy between the contact ion pair and the system at a distance of 5.50 A is significant and calculated to be about 9.5 kcal/mol. Given that the system should be re-
Boudon et al.
6060 The Journal of Physical Chemistry, Vol. 94, No. 15, 1990
------
TABLE V: Computed Variations of the Free Energy at Each Stepa rcc rcc AG fU 2.35 2.20 -2.34 0.06 2.35 2.50 +2.46 0.07 2.65 2.50 -2.44 0.5 1 2.65 2.80 +2.3 1 0.38 2.95 2.80 -2.45 0.20 2.95 3.10 +2.46 0.37 3.25 3.10 -2.09 0.09 3.25 3.40 +2.04 0.07 3.55 3.40 -1.88 0.14 3.55 3.70 +2.18 0.28 3.85 3.70 -2.79 0.34 3.85 4.00 +1.68 0.02 4.15 4.00 +2.41 0.52 4.15 4.30 +2.26 0.33 4.45 4.30 -1.99 0.10 4.45 4.60 +2.28 0.09 4.75 4.60 -1.75 0.40 4.75 4.90 +1.75 0.67 5.05 4.90 -1.53 0.50 5.05 5.20 +1.86 0.30 5.35 5.20 -1.83 0.25 5.35 5.50 +2.08 0.07 5.65 5.50 -0.68 0.02 5.65 5.80 +0.72 0.06 5.95 5.80 -1.53 0.20 5.95 6.10 -1.58 0.25 ~~~
~
----------
OThe distances are in
A
n
v Figure 10. Plots of a Gu+/Gu+ ion pair configuration at a C-C distance of 5.95 A with surrounding water molecules. A
A and energies in kcal/mol.
Figure 9. Plots of a Gu+/Gu+ ion pair configuration at a C-C distance of 3.45 A with surrounding water molecules.
pulsive at large distance according to the Coulombic picture, this would correspond to the dissociation barrier. In contrast to the NMe4+/NMe4+pmf calculated by a similar technique,I4 we find an attractive behavior at short distances. It is noticeable that the well depth calculated for the Gu+/Gu+ pair is significantly lower than the one for the Cl-/Cl- pairI3*l4and large compared to statistical uncertainties. These results support the formation of an intimate ion pair when the system is constrained at less than 5.0 A apart. This intimate ion pair is formed by stacking of the two ions similar to stacking observed for phenyl rings in water.37 By analogy, the Gu+/Gu+ interaction may display hydrophobic character. In addition, the staggeredlike configuration which minimizes the electrostatic and steric repulsion of the ions at this distance is preferred (Figure 9). Without water and by use of the empirical potential function, we calculate at 3.35 and 5.50 A respectively average cation-cation repulsion energies of 68.4 and 50.4 kcal/mol. Thus, the solvation effect overcomes the 18.0 kcal/mol difference in repulsion going from 5.50 to 3.35 A. The in-plane hydration of the -NH2 groups of Gu+ more favorable than the facial hydration is preserved within that intimate ion pair. Interestingly, various X-ray structures of proteins or organic molecules involving guanidinium derivatives in close proximity (37) Ravishanker, G.; Beveridge, D. L. J. Am. Chem. Soc. 1985, 107, 2565-2566.
Figure 11. Plot of the side chains of Arg-33 and Arg-123 belonging to two different protein coils extracted from the rhinovirus crystal structure. The stereoplot is perpendicular to the above plot.
display similar structural features. This situation is quite well illustrated for the human r h i n o ~ i r u s ~(HRV * * ~ ~ 14) where two proteins of the capsid envelope are held together with two close staggered and stacked guanidinium residues (Figure 1l), the C--C distance being 4.1 A. In the complex between immunoglobulin IgA FAB and lysozyme,39@guanidinium residues of Arg-45 and Arg-68 of lysozyme face carboxylate residues of Glu-35 and Glu-50 from FAB. They are in parallel planes with a C-C separation of 4.5 A. In methylguanidinium dihydrogen orthophosphate4' cation-cation stacking is observed with a distance between the mean ion planes of 3.58 A. The guanidinium hydrogen oxalate monohydrate crystals4*the ions are also arranged in a coplanar fashion at 3.20 A apart. The driving forces leading -
(38) Arnold, E.; Vriend, G.; Luo, M.; Griffith, J. P.; Kamer, G.; Erickson, J. W.;Johnson, J. E.; Rossmann, M. G. Acta Crystallogr., Sect. A 1987,43, 346-36 1. (39) Protein Data Bank, Department of Chemistry, Brookhaven National Laboratory, Upton, NJ 11973. (40) Sheriff, S.; Silverton, E. W.; Padlan, E. A.; Cohen, G. H.; Smith-Gill, S. J.; Finzel, B. C.; Davies, D. R. Proc. Natl. Acad. Sei. U.S.A. 1987,84, 8075-8079. (41) Cotton, F. A.; Day, V. W.;Hazen, E. E.; Larsen, S. J. Am. Chem. SOC.1973,95,4834-4840. (42) Adams, J. M. Acta Crystallogr. 1978,B34, 1218-1220.
J . Phys. Chem. 1990, 94, 6061-6069 to such arrangements in the above crystals are different from the forces leading to an intimate ion pair in water. Nevertheless, observation of close proximities between these groups suggests that electrostatic repulsions can be overcome by environmental effects in both cases. The striking similarity between those unrelated structures is noteworthy. At a C-C distance of 5.95 A the relative orientation of the Gu+ ions differs from that in the intimate ion pair. In typical structures perpendicular arrangements (Figure 10) which facilitate electrostatic quadrupolequadrupole interactions between the two groups are found.
Conclusions We have calculated the free energy profile for two Gu+ ions in water and found a clear energy minimum at short distances corresponding to an intimate stacked ion pair. Because of computer limitations, larger separations were not explored and it cannot be concluded whether this intimate pair would be more
6061
stable than a solvent or “infinitely” separated pair. Contrary to the gas-phase electrostatic model, approaching these two cations below approximately 5 A in water led to stabilization. It is stressed that such arrangements will be very sensitive to the molecular environment (e.g., counterions, field of the protein, etc.) and that fluctuations of this environment may either disrupt or enhance this association. Molecular recognition involving reversible binding processes might take advantage of such “counterelectrostatic” attraction. Acknowledgment. Gratitude is expressed to the “groupement scientifique” IBM-CNRS, to the C N R S Computing Center of Strasbourg for support of this research, and to Prof. J. M. Lehn and Dr. Y. Hofflack for helpful discussions. S. Boudon greatly acknowledges fellowship support kindly provided by Rh6nePoulenc and Prof. W. L. Jorgensen for giving the last version of the BOSS computer program. Registry No. Guanidinium, 25215-10-5.
Adsorption in Carbon Micropores at Supercritlcal Temperatures Ziming Tan and Keith E. Gubbins* School of Chemical Engineering, Cornell University, Ithaca, New York 14853 (Received: February 12, 1990)
Nonlocal density functional theory and grand canonical Monte Carlo simulations are used to investigate the adsorption behavior of gases of simple spherical molecules in model carbon micropores at temperatures above the critical value for the gas. In most of the calculations the parameters are chosen to model methane as the adsorbed gas, but some calculations are reported for a model of ethylene. The excess adsorption isotherms (which measure the increased density in the pore, relative to that of the bulk fluid) show a maximum at a particular value of the bulk gas density (pressure). Near the capillary critical temperature these maxima have a cusplike nature, similar to that observed experimentally. We investigate the effect of temperature and pore size on the excess adsorption isotherms and on the maximum excess adsorption. Heats of adsorption for an ideal gas phase are also calculated.
1. Introduction Over the past few years powerful theoretical and computer simulation techniques have been brought to bear on the study of adsorption of fluids on single solid surfaces and in micropores of various simple geometries (slits, cylinders, and spherical cavities). Most of these studies have explored low temperatures, below the critical point of the adsorbate, where a wide variety of novel phase transitions (two-dimensional melting and vaporization, orientational transitions, layering transitions, wetting and prewetting transitions, capillary condensation and melting, and so on) and hysteresis effects occur. Less attention has been paid to adsorption effects at higher temperatures. However, many important applications of adsorption involve adsorbate gases near or above their critical temperature, and there is a need to develop a theoretical understanding of the effects of pore size, geometry, nature of surface, and state conditions on the adsorption, density profiles, diffusion rates, and so on. One example of such an application is the storage of methane (the primary constituent of natural gas) at high densities at near-ambient temperatures, well above the critical temperature of the gas (191 K). Such storage methods are needed for both the transportation of natural gas and for its use as a fuel for vehicles; the very high pressures needed to compress the bulk fluid to high densities at ambient temperature are not as economical as storage at low pressure. The use of microporous materials as a storage medium offers the possibility of achieving high methane densities at low pressures of the bulk gas. Questions relevant to such an application are, what are the optimum temperature and pressure of the bulk methane gas in contact with the micropores, and what is the optimum pore size for such an application? In this paper we examine such questions for a simple model of the methane-carbon micropore system, in
which the carbon pores are of parallel wall geometry and have surfaces that approximate the basal plane of graphite. Other relevant questions, not addressed here, include the effects of modifying the shape and composition of the solid surface. For our calculations we use density functional theory, together with grand canonical Monte Carlo simulations. The densityfunctional theory treats the fluidsolid system as an inhomogeneous fluid in a static external field exerted by the solid. It involves two major approximations, the smoothed density approximation (SDA),’q2which is used to calculate the contribution to the free energy from the short-range repulsive forces, and the mean field approximation (MFA),’ which is used to estimate the contribution from the longer-range attractive forces. In the SDA the free energy density for the fluid with repulsive forces is given by that for a uniform hard-sphere fluid at a nonlocal smoothed density. We refer to this theory as the nonlocal mean-field density functional theory (NLMFT). The NLMFT has been employed successfully in studies of the structure of confined fluids$J capillary condensation,6-’ layering transitions,6 and the structure and se(1) Tarazona, P. Phys. Rev. A 1985,31,2672; 1985,32,3148. Tarazona, P.; Marini Bettolo Maconi, U.; Evans, R. Mol. Phys. 1987, 60, 573. (2) Tan, Z.; Marini Bettolo Marconi, U.; van Swol, F.; Gubbins, K. E. J . Chem. Phys. 1989, 90, 3704. (3) Sullivan, D. E.; Telo da Gama, M. M. Wetting Transitions and Multilayer Adsorption at Fluid Interfaces. In Fluid Interfacial Phenomena; Croxton, C. A., Ed.; Wiley: New York, 1986, and references therein. (4) Heffelfinger, G.S.;Tan, Z.; Marini Bettolo Marconi, U.; van Swol, F.; Gubbins, K. E. Mol. Simul. 1988, 2, 393. ( 5 ) Ball, P. C.; Evans, R. J . Chem. Phys. 1988, 89, 4412. (6) Tarazona, P.; Marini Bettolo Marconi, U.; Evans, R. Mol. Phys. 1987, 60, 573. (7) Peterson, E. K.; Gubbins, K. E.; Heffelfinger, G.S.;Marini Bettolo Marconi, U.; van Swol, F.J. Chem. Phys. 1988,88,6487.
0022-3654/90/2094-6061%02.50/00 1990 American Chemical Society
