Computer Simulation of Water near the Surface of Oligo(ethylene

Langmuir 2000, 16, 8829-8841

8829

Computer Simulation of Water near the Surface of Oligo(ethylene glycol)-Terminated Alkanethiol Self-Assembled Monolayers† Alexander J. Pertsin and Michael Grunze* Angewandte Physikalische Chemie, Universita¨ t Heidelberg, Im Neuenheimer Feld 253, D-69120 Heidelberg, Germany Received March 7, 2000. In Final Form: June 8, 2000 A Grand canonical Monte Carlo technique is used to simulate the behavior of the TIP4P model of water near the surface of an oligo(ethylene glycol) (OEG)-terminated alkanethiol self-assembled monolayer (SAM). Two structural modifications of the SAM are studied in an effort to shed light on the difference in their resistance to adsorb proteins. One modification (hereafter h-SAM), which is stable on the Au substrate, is characterized by a helical conformation of the OEG tails and a lower areal density. The other, higher density modification (hereafter t-SAM), which is experimentally observed on the Ag substrate, has their OEG tails in the all-trans conformation. Simulations show that water molecules can penetrate into the near surface region of the h-SAM in appreciable amounts to bring about substantial conformational disordering of the SAM, in agreement with recent sum frequency generation experiments. About 75% of the topmost oxygen atoms and even 2% of the next to the topmost oxygen atoms of the OEG tails prove to be involved in hydrogen bonds with water molecules. By contrast, the t-SAM is much more resistant to the penetration of water molecules. On interaction with water, it retains a good conformational order and shows a noticeably lower surface density of hydrogen bonds with water molecules. Both the h- and t-SAM surfaces produce water density oscillations that extend at least 3 to 4 molecular diameters into the water bulk. The density oscillations are less pronounced near the h-SAM, whose surface is ill-defined because of the conformational disordering of the OEG tails. The effect of the SAM on the structural characteristics of water is short range and practically vanishes at separations larger than 1 to 2 molecular diameters away from the SAM surface.

1. Introduction Self-assembled monolayers (SAMs) formed by oligo(ethylene glycol) (OEG)-terminated alkanethiol molecules are known for their ability to resist protein adsorption.1 Unlike poly(ethylene glycol) (PEG) brushes, whose protein resistance is rationalized in terms of a high conformational freedom of the PEG chains,2 the OEG-terminated SAMs represent closely packed molecular assemblies with a restricted conformational freedom of their constituent molecules. A remarkable property of the SAMs is that their protein adsorption properties are strongly dependent on the conformation of the OEG tails.3 In monolayers prepared on an Au substrate, the OEG tails assume a conformation close to a 7/2 helix and the resulting monolayer is protein resistant. By contrast, monolayers self-assembled on Ag show a planar all-trans conformation of the OEG tails, and the monolayer does adsorb protein. (Hereafter, the two modifications of the SAM will be referred to as the helical and all-trans ones or, simply, the h-SAM and t-SAM, respectively). An attempt to explain the observed difference in protein adsorption properties between the h- and t-SAMs was undertaken by Wang et al.4 using ab initio Hartree-Fock † Part of the special Issue “Colloid Science Matured, Four Colloid Scientists Turn 60 at the Millennium”.

(1) Prime, K. L.; Whitesides, G. M. J. Am. Chem. Soc. 1993, 115, 10714. (2) (a) Jeon, S. I.; Andrade, J. D. J. Colloid Interface Sci. 1991, 142, 159. (b) Taunto, H. J.; Toprakcioglu, C.; Fetters, L. J.; Klein, J. Nature 1988, 332, 712. (c) Sleifzer, I. Curr. Opin. Solid State Mater. Sci. 1996, 2, 337. (d) Halperin, A. Langmuir 1999, 15, 2525. (3) Harder, P.; Grunze, M.; Dahint, R.; Whitesides, G. M.; Laibinis, P. E. J. Phys. Chem. B 1998, 102, 426. (4) Wang, R. L. C.; Kreuzer, H. J.; Grunze, M. J. Phys. Chem. 1997, 101, 9767.

calculations. It was found that a single helical OEG strand was capable of forming strong bridged hydrogen bonds with a water molecule, which involve two successive oxygen atoms of the strand. By contrast, the planar alltrans conformation only allowed the formation of single hydrogen bonds involving one OEG oxygen atom and one water hydrogen atom. The behavior of water near the SAM surface was modeled with small clusters comprising up to 20 water molecules and up to 12 rigid OEG strands packed on a hexagonal lattice as observed in real SAMs. It turned out that the rigid t-SAM model could not form hydrogen bonds with water at all, except for weak C-H‚‚‚O bonds with the terminal methyl hydrogens. The h-SAM was capable of forming single hydrogen bonds between the topmost oxygen atoms and water hydrogen atoms. The ability of the helical OEG strands to form bridged hydrogen bonds with water was lost because the densely packed rigid SAM structure did not allow water molecules to penetrate into the SAM to reach the next to the topmost oxygen atoms. It was speculated that the bridged hydrogen bonds can nevertheless be formed in less dense regions of the h-SAM, e.g., near defects. Recent sum frequency generation (SFG) experiments5 have provided clear indications that the conformations of the OEG tails in the OEG-terminated SAMs are drastically affected when the SAMs are exposed to liquid water. An analysis of the SFG spectra has shown that water molecules penetrate fairly deep into the SAM to cause conformational disordering of the OEG moieties. Clearly, the assumption of a rigid regular SAM structure in modeling the properties of the SAMs is oversimplified and more realistic models, which involve flexible OEG chains (5) Zolk, M.; Eisert, F.; Buck, M.; Eck, W.; Herrwerth, S.; Pipper, J.; Grunze, M. Submitted for publication in Langmuir.

10.1021/la000340y CCC: $19.00 © 2000 American Chemical Society Published on Web 08/24/2000

8830

Langmuir, Vol. 16, No. 23, 2000

and thermal motion, are called for. Further arguments against the use of rigid models for the OEG strands have been provided by recent ab initio calculations6 based on density functional theory (DFT) with nonlocal gradient corrections. It has been shown that the interaction with water molecules alters substantially the stability order for the OEG conformers, and so the OEG chains can in no case be treated as rigid. In the present paper, which is issued together with the above-cited SFG and quantum chemical studies,5,6 we construct a molecular-level description of the OEGterminated SAM/water interphase region based on the Monte Carlo (MC) computer simulation technique. The simulation model involves up to ∼1000 water molecules and 36 flexible OEG chains periodically replicated in two dimensions parallel to the SAM surface. The intra- and intermolecular interactions in the system are described using a classical atomistic force field based on the TIP4P model of water7 and ab initio MP2 level calculation results for short OEG homologues8 and an OEG/water complex.9 The MC technique is favored over the more popular method of molecular dynamics (MD) for two reasons. One is that the chain topology of OEG-terminated alkanethiol molecules allows an efficient sampling of the molecular configurational space using a rotational displacement algorithm suggested in our early work.10 Since the rotational displacements are “unphysical” (i.e., inconsistent with the natural time evolution of the system), they can only be implemented within the MC technique. The other reason to choose the MC technique is that it can be readily adapted to calculation of averages in the Grand canonical (GC) ensemble, where the chemical potential µ, volume V, and temperature T are fixed, while the number of molecules, N, fluctuates. The use of the GC ensemble is of exceptional importance in simulation of the SAM/ water interphase region because of large density and pressure gradients. An additional advantage of the GCMC sampling is that it will occasionally destroy water molecules strongly bound to the SAM through hydrogen bonds and then create unbound molecules far away from the interface. This will increase the sampling efficiency and, in addition, break down metastable hydrogen-bonded structures near the SAM surface. All the calculations are made with SAMs composed of methoxy tri(ethylene glycol) (hereafter, EG3-OMe)-terminated alkanethiol molecules, S(CH2)n-(OCH2CH2)3OCH3. For the sake of computational simplicity, the alkane chain is taken short, n ) 3. A few comparative simulations of the SAMs with n ) 3 and n ) 10 did not reveal a perceptible effect of the alkane chain length on the behavior of the EG3-OMe tails. 2. Methods 2.1. Sampling Procedure. In our simulations, an MC run consists of n passes, each composed of N + Ns moves, where N is the current number of water molecules and Ns the number of OEG-terminated alkanethiol molecules in the SAM. There are four types of moves: (1) a displacement of a water molecule, which combines a translational and rotational displacement; (2) an insertion of a water (6) Wang, R. L. C.; Kreuzer, H. J.; Grunze, M. Phys. Chem. Chem. Phys., in press. (7) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (8) Smith, G. D.; Jaffe, R. L.; Yoon, D. Y. J. Phys. Chem. 1993, 97, 12752. (9) Bedrov, D.; Pekny, M.; Smith, G. D. J. Phys. Chem. B 1998, 102, 996. (10) Pertsin, A. J.; Hahn, J.; Grossmann, H.-P. J. Comput. Chem. 1994, 15, 1121.

Pertsin and Grunze

molecule; (3) a deletion of a water molecule; (4) a move of a SAM molecule, which changes both the position and all the conformational parameters of the molecule using a rotational displacement procedure10 described in the next subsection. The number of molecules in the SAM is fixed, Ns ) const, so that no attempt to delete or create a molecule comprising the SAM is undertaken. The four allowed types of moves are attempted with probabilities Pw, Pi, Pd, and Ps, respectively, such that

Pw + Pi + Pd + Ps ) 1

(1)

The probabilities, P, are calculated from the following input parameters:

Rw/s ) (Pw + Pi + Pd)/Ps

(2)

Ri ) Pi/(Pw + Pi + Pd)

(3)

Ri/d ) Pi/Pd

(4)

The meaning of these parameters is apparent: Rw/s specifies the distribution of attempted moves between the water and SAM molecules, Ri represents the fraction of insertions among all moves of the water molecules, and Ri/d specifies the ratio between water insertion and deletion attempts. The numerical values used for these parameters are Rw/s ) 10, Ri ) 0.7-0.8, and Ri/d ) 15. Displacements. A water molecule, considered as a rigid body, is displaced in a standard way, by translating its center of mass by a random vector δr and then rotating it by a random angle δγ about one of the three space-fixed axes chosen at random.11 The OEG-terminated alkanethiol chains are treated as being flexible but subject to bondlength constraints. In addition, the position of the two hydrogen atoms attached to a carbon atom, C*, is constrained by maintaining the H-H vector perpendicular to the C-C*-X plane and the midpoint of this vector on the bisector of the C-C*-X angle (X ) S, C, or O). The bond-length constraints are implemented using the rotational displacement procedure.10 In this procedure, an elementary displacement represents a rotation of a chain backbone atom, C*, around the virtual bond linking the two neighboring backbone atoms, C and X. Because of the constraints imposed upon H-H vectors, the hydrogens attached to C* are rotated together with C* as a rigid body. For the same reason, the hydrogens at C and X (when X ) C) prove to be also involved in the rotation. The headgroup sulfur atom is rotated around the closest C-C bond, while the terminal carbon atom is rotated around the C-O bond linking the -OMe group to the rest of the chain. One move of a SAM molecule is defined to include successive rotational displacements of all backbone atoms, starting with the headgroup sulfur and ending with the terminal carbon atom. Unlike rotations around single bonds, which may result in large displacements of the molecular tails and in intermolecular overlaps, rotations around the above-defined virtual bonds affect the structure of a long-chain molecule only locally. The rotational displacement procedure is therefore particularly suitable for dense systems such as SAMs, where a representative sampling involves serious problems because of “pseudononergodicity”12 of their configurational space.13 (11) Barker, J. A.; Watts, R. O. Chem. Phys. Lett. 1969, 3, 144. (12) Wood, W. W. In Physics of Simple Liquids; Temperley, H. N. V., Rowlinson, J. S., Rushbrooke, G. S., Eds.; North-Holland: Amsterdam, 1968. (13) Pertsin, A. J.; Grunze, M. Langmuir 1994, 10, 3668.

Water Behavior near SAMs

Langmuir, Vol. 16, No. 23, 2000 8831

All the displacement moves are handled using the normal Metropolis algorithm; i.e., a displacement is accepted with a probability equal to

min(1, exp(-∆UN/kT))

(5)

where ∆UN is the change in the potential energy of the system due to the given displacement. Insertions. In the conventional grand canonical Monte Carlo (GCMC) technique,14 as applied to a system of N rigid molecules, an insertion move is accepted with a probability given by

(

min 1,

)

N exp[(µ′ + UN - UN+1)/kT] N+1

(6)

In this equation, µ′ is the excess (nonideal) part of the chemical potential, µ′ ) µ - µ0, where µ0 is the contribution to µ due to the translational and rotational degrees of freedom of the water molecules (with the intermolecular interactions switched off). The explicit expression for µ0 in terms of the mass, m, and the moments of inertia of the water molecule, Ii (i ) 1, ..., 3), is

µ0 ) kT ln

[ (

)]

)(

N h 8π2V (2πmkT)1/2

3

h (2πkT(I1I2I3)1/3)1/2

3

(7)

At high densities, the usual GCMC technique becomes useless because the probability of making a successful insertion attempt becomes very small. To improve the efficiency of insertions, we employ an excluded volume mapping method.15 In this method, the whole volume accessible to the system is divided into a large number of small cubes. A cube is considered to be forbidden for insertion of a new water molecule if there is an oxygen or a carbon atom closer than a specified distance to all points within the cube. (The particular distances used are 2.5 Å for O and 2.8 Å for C.) As the simulation proceeds, the state of the cubes is checked and updated. The probability of accepting an insertion attempt is modified to

(

min 1,

)

NfN exp[(µ′ + UN - UN+1)/kT] (N + 1)Ri/d

(8)

where fN is the fraction of cubes allowed for insertion and Ri/d is the factor that corrects for the bias when insertions and deletions are attempted with different frequencies. Deletions. A deletion attempt is accepted with a probability given by

(

min 1,

)

Ri/d exp[(-µ′ + UN - UN-1)/kT] fN-1

(9)

Compared to the case of insertion, a deletion of a water molecule involves an additional problem because of the necessity of estimating fN-1 prior to the deletion attempt. This problem arises when no state with (N - 1) water molecules has been yet sampled and hence no estimate is available for fN-1. In such cases, we evaluate fN-1 by linear extrapolation from the equation

fN-1 ) fN + 〈fN〉 - 〈fN+1〉

(10)

(14) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1993. (15) Stapleton, M. R.; Panagiotopoulos, A. J. Chem. Phys. 1990, 92, 1285.

Figure 1. The cell used in simulations of the SAM/water interphase region. The hydrogen atoms in the OEG-terminated alkanethiol molecules are not shown for clarity.

where the angular brackets denote ensemble averaging. (The event when both fN-1 and 〈fN+1〉 are not available may happen only once in the beginning of the GCMC run. In that event, fN-1 is simply taken to be equal to fN.) Boundary Conditions. The system comprising N water molecules and Ns EG3-OMe terminated alkanethiol molecules is confined within a rectangular cell with dimensions Lx, Ly, and Lz, as shown in Figure 1. The SAM substrate coincides with the bottom face of the cell. The side lengths in the x and y dimensions are multiples of the SAM periods a and b ) a31/2. In the x and y directions, standard periodic boundary conditions14 are used. As to the z direction, the usual practice in simulation of liquids at planar interfaces is to confine liquid between two parallel identical interfaces perpendicular to z. In our case, this would require simulation of two SAMs, one at the bottom and the other at the top of the simulation cell, which is computationally expensive. Another possible solution is to confine water at the top with a hard wall. In this case, however, the behavior of water in the SAM/water interphase region may be distorted by the structuring effect of the wall.16,17 To avoid these problems, we here suggest “gliding plane” boundary conditions, which combine a mirror plane at z ) Lz with a translation by half period Lx along x. That is, the space above the top face of the simulation cell is filled with a mirror image of the system shifted along x by Lx/2. (The shift is necessary to avoid straightforward correlations of a particle and its image across the mirror plane.) The corresponding transformation of coordinates is (16) Valleau, J. P.; Gardner, A. A. J. Chem. Phys. 1987, 86, 4162. (17) Grigera, J. R.; Kalko, S. G.; Fischbarg, J. Langmuir 1996, 12, 154.

8832


x′ ) x + kLx/2

Pertsin and Grunze

(11)

(k ) 1 for x < Lx/2 and k ) -1 otherwise).

y′ ) y

(12)

z′ ) 2Lz - z

(13)

As a molecule leaves the simulation cell through the top face, its image, whose coordinates are given by eqs 1113, enters the cell through the other half of the top face. Compared to the case of two parallel interfaces, the gliding plane boundary conditions save nearly half the CPU time, though they may, in principle, induce artificial correlations between the positions of particles near the top of the simulation cell. 2.2. Force Field. The choice of the force field is of critical importance in computer simulations using classical analytical potentials. This particularly applies to PEG and its low-molecular homologues because the preferred conformations of their flexible chains are governed by very delicate effects. Experimentally, there is a preference for a gauche conformer about the C-C bond (the so-called gauche effect) and for a trans-gauche-trans (tgt) conformer of the C-O, C-C, and C-O bond sequence.18 Recent ab initio quantum chemical studies6,19,20 of short OEG homologues have shown that the tgt conformer has practically the same or a slightly lower energy compared to the ttt one. In our simulations, the intra- and intermolecular interactions within the OEG-terminated SAM are described based on a force field developed by Smith, Jaffe, and Yoon8 (SJY) by fitting to ab ibitio MP2 level calculation results for 1,2-dimethoxyethane (DME) and diethyl ether (DEE).20 The SJY force field comprises bond stretching, bond angle bending, torsional, and nonbonded contributions. The nonbonded terms, in their turn, include electrostatic and van der Waals atom-atom interactions, as described by Coulomb and Buckingham (6-exp) potentials, respectively. The partial atomic charges in the SJY force field were first calculated from Mulliken populations and then proportionally scaled to fit the dipole moment of the tt conformer of DEE. The other parameters of the force field were adjusted so as to attain the best agreement between the force field and ab initio results for the energies and geometries of DME and DEE. The resulting force field was finally tested in a stochastic dynamics simulation of gas-phase DME by comparing the calculated conformer populations and pair distribution functions with those deduced from an electron diffraction experiment. We subjected the SJY force field to an independent test by calculating the energy of various conformers of dyglime and found an excellent agreement with ab initio MP2 level results.19 For the C-C-C-C and C-C-C-O torsions and for the C-C-C bend, which are absent in the SJY force field but are required to treat the alkane chain and its junction with the EG3-OMe tail, we use the analytical potentials from the alkyl ether force field,21 which was employed as the starting point in deriving the SJY potentials. To describe the interaction of the monolayer with the metal (18) (a) Yoshida, H.; Kaneko, I.; Matsuura, H.; Ogawa, Y.; Tasumi, M. Chem. Phys. Lett. 1992, 196, 601. (b) Inomate, K.; Abe, A. J. Phys. Chem. 1992, 96, 7934. (c) Takasi, K.; Abe, A. Polym. J. 1985, 17, 641. (19) Gejji, S. P.; Tegenfeldt, J.; Lindgren, J. Chem. Phys. Lett. 1994, 226, 427. (20) Jaffe, R. L.; Smith, G. D.; Yoon, D. Y. J. Phys. Chem. 1993, 97, 12745. (21) Sorensen, R. A.; Liau, W. B.; Kesner, L.; Boyd, R. H. Macromolecules 1988, 21, 200.

substrate, the force field is complemented by analytical potentials from our previous study of alkanethiol/Au(111) SAMs13 (force field I). The headgroup-substrate interactions are described by a periodic “surface corrugation” potential, while the methylene groups of the alkanethiol chain interact with the substrate through (3-12) potential functions. The interactions of the EG3-OMe tails with the metal substrate are neglected. The goodness of the whole force field describing the interactions within the SAM was tested in our recent calculations of the equilibrium structure of the h- and t-SAMs.22 The calculations correctly reproduced the transition from the helical conformation of the EG3-OMe chains on Au to the all-trans conformation on Ag, as well as the relevant changes in the vibrational spectra. For water-water interactions, we employ the TIP4P model potential,7 which is perhaps the most widely used and carefully tested potential in simulations of liquid water. The TIP4P model for water comprises four interaction sites roughly corresponding to the oxygen, two hydrogen atoms, and a quasi-atom located 0.15 Å away from the oxygen atom along the bisector of the H-O-H bond angle. The two hydrogen atoms carry a charge of 0.52e each, while the quasi-atom carries a balancing charge of -1.04e. The oxygen atoms interact through a Lennard-Jones potential. For cross interactions involving water and EG3-OMe alkanethiol molecules, use is made of a force field recently derived by Bedrov et al.9 by fitting their ab initio MP2 level results for the DME/water complex. We tested this force field independently by comparing the force field results for single and bridge-mode hydrogen bonds between a EG3-OMe tail and a water molecule with the relevant ab initio MP2 level data.4 The discrepancies proved to be within 15%. More recent ab initio calculations6 using a larger basis set have however shown that the binding energies found in the earlier work4,9 for the single hydrogen bond are 5-8 kJ/mol too large (see ref 6 for more detail). Fortunately, the parametrization of the Bedrov et al. force field9 was made based on a bridged hydrogen-bonding configuration, where the new6 and old4,9 binding energies agree to within 10%. Nevertheless, we are inclined to regard the force field used in our work in describing the OEG/water interactions only as a first approximation. One of the reasons is that this force field was parametrized not specifically for EG3-OMe but for shorter OEG homologues. In our subsequent publications, our GCMC simulations will be extended using a refined force field based on accurate ab initio DFT results for EG3-OMe/n‚H2O (n ) 0, 1, 2) complexes.6 Since the force field chosen for our simulations contains Coulombic terms, which are characterized by a slow distance dependence (∼r-1), the problem of efficient summation of the electrostatic interactions takes on great importance. In systems with three-dimensional periodic boundary conditions, this problem can be successfully solved using various “accelerated convergence” techniques based on the ideas of Ewald and Bertaut.23 The convergence problem is more difficult to cope with when the system is periodic in two dimensions, while being of a finite thickness in the third dimension. The pertinent formulas24 are rather complicated and provide a satisfactory convergence only when the system size in the “nonperiodic” dimension is much less than the periods in (22) Pertsin, A. J.; Grunze, M.; Garbuzova, I. A. J. Phys. Chem. B 1998, 102, 4843. (23) Bertaut, F. J. Phys. Radium 1952, 13, 499. (24) Hautman, J.; Klein, M. L. Mol. Phys. 1992, 75, 379.


Langmuir, Vol. 16, No. 23, 2000 8833

the other two dimensions. In our case, this condition obviously fails. The convergence problem is particularly important in evaluation of the electrostatic energy of the SAM because of the presence of the long-range translational order and a nonzero dipole moment of the unit cell. In our simulations, the electrostatic interactions within the SAM are handled using a modified direct summation scheme suggested by Heyes and van Swol.25 In this scheme, the electrostatic potential at a point r ) {X, Y, Z} from a system of M charges qj replicated in two dimensions with periods Lx and Ly is represented by a truncated direct sum M

qj ∑ j)1

∑

|r - rj - i1Lx - i2Ly|-1

i1,i2(i12+i22eic2)

and a correction term, C(ic), which provides an approximate estimate for the discarded long-range contributions (i12 + i22 > ic2). The approximation involves, in particular, the replacement of the discrete charges in each jth sublattice by a uniform two-dimensional distribution of density qj/|Lx × Ly|. The final equation for the correction term is of the form M 1 C(ic) ) - Ly-3ic-1 qj[Cx(x - xj)2 + Cy(y - yj)2 + 2 j)1 Cz(z - zj)2] (14)

∑

where the coefficients Cx, Cy, and Cz are expressed in terms of complete elliptical integrals of the first and second kind, K(2) and E(2), with ) (Lx2 + Ly2)1/2/Lx (Lx > Ly). As more and more periodic cells are included in the direct sum, the magnitude of the correction term decreases, which is reflected in the inverse ic dependence in eq 14. This feature is only intrinsic in the special case of two-dimensional periodicity and a finite thickness in the third dimension. For a system periodic in three dimensions, the corresponding long-range correction is independent of ic and, in addition, depends on the shape of the unit cell, as defined by the periods Lx and Ly. The efficiency of the Heyes-van Swol scheme, when applied to EG3-OMe terminated SAMs, has been demonstrated in our recent work.22 For the interactions involving water molecules, we use a spherical cutoff (SC) scheme, which simply discards the interactions of molecules lying farther than a specified distance from each other. To avoid discontinuity in the interaction potential, the electrostatic interactions are multiplied by a switching function of the form

{

r < rl, 1, 2 2 2 2 S(r) ) (r - rl )/(rh - rl ), rh > r > rl, r > rh, 0,

2.3. Ensemble Averages and Structural Quantities. In the course of the GCMC run we calculate the average potential energy of the system, 〈U〉, the average number of particles, 〈N〉, and also the z component of the pressure, Pz, at z ) Lz. The structure of the SAM-water interphase region is monitored using the following distribution functions and parameters: ‚The z profile of the density of water, d(z), defined in terms of the density of the water O sites as d(z) ) do(z)18/16. ‚The xy maps of the distribution of water O atoms in slices perpendicular to z, d(x,y,z). ‚The pair distribution functions of the water O atoms in slices perpendicular to z, d(x,y,z), and their averages, Fz(r), obtained by averaging F(r,z) over z in some selected intervals. ‚The radial distribution functions of the water O and H atoms around the different O atoms of the EG3-OMe tails, Fi(r) (i ) 1, ..., 4), where i ) 1 refers to the deepest oxygen atom (at the junction of the alkanethiol chain with the EG3-OMe tail) and i ) 4 to the topmost oxygen atom. ‚The z dependence of the distribution of the angle θµ formed by the molecular dipole moment and the z axis

p(θµ,z) ) 〈δ(θµ - θµ′(z)〉/sin θµ

(16)

‚The z dependence of the orientational order parameters

S01(z) ) 〈cos θµ(z)〉

(17)

S02(z) ) 〈3 cos2 θµ(z) - 1〉/2

(18)

which represent, in essence, the first and second moments of the distribution of cos θµ.27 ‚The distribution of the dihedral angles, τk, in the EG3OMe tails

pk(τ)

(k ) 1, ..., 10)

(19)

where k ) 1 refers to the C-C-O-C dihedral angle at the junction of the alkanethiol chain and the EG3-OMe tail, k ) 2 refers to the next dihedral, C-O-C-C, in the direction of the tail end, and so on, up to k ) 10 corresponding to the C-C-O-Me dihedral in the terminal methoxy group. ‚The z dependence of the fraction of water molecules involved in n hydrogen bonds with its neighbors, hn(z) (n ) 0, ..., 6), and the average number of hydrogen bonds per water molecule, 6

〈n(z)〉 )

∑ nhn(z)

(20)

n)0

(15)

where rl ) 8.0 Å and rh ) 8.5 Å. The applicability of the SC scheme to simulation of water between planar walls has been demonstrated by Shelley and Patey.26 It was shown that this simple scheme provides nearly the same results as the Ewald method and allows one to avoid a pathological anisotropy in the orientational distribution of water molecules, which is observed with the minimal image and cylindrical cutoff schemes.16 (25) Heyes, D. M.; van Swol, F. J. Chem. Phys. 1981, 75, 5051. (26) Shelley, J. C.; Patey, G. N. Mol. Phys. 1996, 88, 385.

two water molecules are considered to form a hydrogen bond if their interaction energy is lower than -9.41 kJ/ mol.19 In the calculation of the z dependences, the whole height of the simulation cell, Lz, is divided into 80 slices parallel to the xy plane. During the GCMC run, the above-listed quantities are averaged within each individual slice and then referred to the z coordinate of its center. (27) The utility of S02 can be readily understood by considering (i) a distribution with two preferred orientations strictly along z and -z, (ii) a uniform distribution with no preferred orientation, and (iii) a distribution where all the molecular dipole moments are perpendicular to z. In all the three cases, S01 ) 0 whereas S02 ) 1 in case (i), 0 in case (ii), and -0.5 in case (iii).

8834


Figure 2. Hydrogen bond distribution per water molecule for bulk water. The crosshatched bars and full squares refer to this work and Jorgensen’s results,7 respectively.

3. Results and Discussion 3.1. Reference Systems. Prior to simulations of the SAM/water interphase region, we studied three systems which can be regarded as references. These are an isolated h-SAM, bulk water, and water at a hard wall. Isolated SAM. The results for an isolated h-SAM have been described in detail elsewhere,28 and here we only note that at room temperature the h-SAM retained a good translational, orientational, and conformational order, in agreement with FTIRAS experiments in a waterless atmosphere.3 The nine dihedrals in the OEG tails each had a Gaussian-like distribution with a width ranging from 20° to 30°. The mean values of the dihedrals formed a tgt tgt tgt sequence, although all three O-C-C-O angles were systematically larger (by 20°-25°) than the ideal gauche angle (60°). Bulk Water. The simulations of bulk water were carried out using periodic boundary conditions in all three dimensions with periods Lx ) 30, Ly ) 26, and Lz ) 32.7 Å. The simulations were performed at some selected values of µ′ in the range from -23.0 to -25.5 kJ/mol. A value of Pz close to the atmospheric pressure was reproduced at µ′ ) -25.1 kJ/mol, in accordance with the available literature data for the TIP4P model of water.26 The average number of water molecules in the simulation cell proved to be 853, which corresponds to a density of 0.996 g/cm3, in agreement with experiment. The range of sampled values of N in the production part of the GCMC run was from 836 to 896. The calculated average potential energy was -42.3 kcal/mol, also in agreement with previous simulations.26 No anisotropy in the angular distribution of water molecules was observed: Both S01 and S02 showed only slight oscillations around zero, while p(θµ,z) oscillated around 0.5. The hydrogen bond distributions hn(z) were independent of z and showed an excellent agreement with the constant-NPT MC results reported by Jorgensen et al.7 for the TIP4P model of water at room temperature and atmospheric pressure (Figure 2). The average number of hydrogen bonds per water molecule, 〈n〉, calculated from eq 20 was 3.5. Water at a Hard Wall. The simulation of water near a hard wall was carried out using the same dimensions of the simulation cell as specified above for bulk water. The hard wall was placed at the bottom face of the cell. The molecules near the top face of the cell were treated using the gliding plane boundary conditions given by eqs 11-13. A nearly atmospheric pressure at Lz was observed (28) Pertsin, A. J.; Grunze, M.; Kreuzer, H. J.; Wang, R. L. C. Phys. Chem. Chem. Phys. 2000, 2, 1721.

Pertsin and Grunze

at µ′ ) -24.7 kJ/mol. This value is close to µ′ for bulk water, which suggests that the system size is adequately large to describe the transition from the hard wall/water interface to bulk water. The average number of water molecules at µ′ ) -24.7 kJ/mol proved to be close to 800. The calculated water density profile and the distribution functions were similar to those found in simulations of water confined between two planar hydrophobic walls.16,17,26 In particular, near the hard wall, S01 was slightly positive, which showed a tendency for the molecular dipole moment µ to point away from the surface, while S02 was negative, which pointed to a trend to have a significant component parallel to the surface (Figure 3a). The same trends were apparent from the behavior of p(θµ,z) at small z (Figure 3b). In slices close to the wall, a clear maximum near θµ ) 90° was observed, as well as a more frequent occurrence of angles close to 0° compared to angles close to 180°. Far from the wall, the surface-induced orientational ordering vanished. This supported the conclusion26 that the pathological orientational anisotropy observed in the early simulations16 was due to the minimal image convention used in calculating the water-water interactions. Close to the wall, hn(z), corresponding to three, four, and five hydrogen bonds per water molecule, rapidly reduced and 〈n〉 tended to 2. An important result was that the pair distribution function Fz(r) averaged over slices adjacent to the top face of the simulation cell was practically identical with that averaged over the bulk slices and also with that found in the simulations of bulk water. This suggests that the correlations introduced by the gliding plane boundary conditions were not significant for the system size used. 3.2. SAM/Water Interphase Region. In the simulations of the SAM/water interphase region, the starting configuration of the system was constructed by combining the lowest energy configuration of the SAM, as found in our static lattice energy calculations,22 with a typical configuration of water molecules found in our simulations of bulk water. The simulation cell included 36 OEGterminated alkanethiol molecules, while the average number of water molecules was dependent on the lattice period of the SAM (〈N〉 ) 950 for the h-SAM and 〈N〉 ) 715 for the t-SAM). The starting configuration of water molecules was initially attached to the starting configuration of the SAM so that the closest distance between a water molecule and a non-hydrogen atom of the SAM was equal to the respective van der Waals separation. In the course of the GCMC run, the water molecules were allowed to penetrate into the SAM, and the number of water molecules in the system was allowed to fluctuate in accordance with the probabilities of insertion and deletion attempts, as given by eqs 8 and 9. The length of the GCMC run ranged from 106 to 107 passes depending on the convergence of the ensemble averages and structural quantities. A value of Pz close to the atmospheric pressure was obtained at µ′ ) -24.3 kJ/mol. Helical Modification. The basic results for the h-SAM are summarized in Figures 4-11. As seen from the distribution of the dihedral angles in the EG3-OMe tails (Figure 4) and from a typical snapshot of the interphase region (Figure 5), the contact with water brings about a substantial conformational disordering of the h-SAM. Although the gauche O-C-C-O and trans C-C-O-C angles remain most populated, the contribution of “wrong” conformers is quite appreciable. The disordering propagates deep into the h-SAM and affects not only the terminal fragments of the OEG tails, which are in direct contact with water, but also fragments close to the alkane chain. The density profile of water strongly overlaps with that


Langmuir, Vol. 16, No. 23, 2000 8835

Figure 3. Water at a hard wall: (a) the orientational order parameters; (b) the angular distributions p(θµ,z) for five ∼0.4 Å thick slices close to the hard wall.

Figure 4. Distribution of the dihedral angles in the EG3-OMe chains for the h-SAM in contact with water. The angles are numbered from 1 to 10 in succession, starting with the C-C-O-C dihedral at the junction of the alkanethiol chain and the OEG tail and ending by the C-C-O-CH3 dihedral in the terminal methoxy group.

of the non-hydrogen atoms of the h-SAM (Figure 6, solid lines), which indicates that water penetrates into the top layers of the h-SAM in appreciable amounts. The average density of water in the penetration region, calculated by the integration of d(z) over z in the respective interval, proves to be as large as 23% of the density of bulk water. Note that for a 25 Å thick water layer, including the penetration range, the average density is 91% bulk density,

which compares well with the result of recent neutron reflectivity measurements (86%).29 The fact of conformational disordering (amorphization) of the SAM when exposed to water is supported by the SFG measurements.5 The most striking result of these measurements is a strong reduction of the asymmetric (29) Steitz, R.; Pipper, J.; Schwendel, D.; Grunze, M. In preparation.

8836


Pertsin and Grunze

Figure 5. A snapshot of a fragment of the h-SAM/water interphase region close to the interface.

Figure 6. Density distribution of water and non-hydrogen atoms of the SAM for the h- and t-SAM/water interfaces.

in-plane methyl mode intensity after immersing the SAM in water. This is a clear indication that the methyl end group, which has a distinct preferential orientation in air, assumes a basically random orientation when the SAM is in contact with water. Perceptible changes are also observed in the intensity of the methylene modes,5 which suggests that water not only affects the terminal methyl groups but also penetrates deep into the SAM. Due to the conformational disordering of the OEG tails, the surface of the h-SAM is ill-defined. As seen from the snapshot in Figure 5, the OEG tails differ greatly in length, so that some individual tails protrude into the water. Such a “diffuseness” of the h-SAM surface can also be appreciated from a fairly gentle slope of the edge of the SAM

density distribution in Figure 6: The range of occurrence of the terminal methyl carbon atom extends from 17.8 to 20.5 Å. The water molecules easily reach the level of the topmost oxygen atoms (O4) in the h-SAM. The separation of the lowest water molecules from the level of the next to the topmost atoms (O3) is ∼3 Å, which is close enough to form a hydrogen bond. The formation of hydrogen bonds between the water molecules and the oxygen atoms of the EG3-OMe chains is evident from the behavior of F3(r) and F4(r), which show the radial distribution of the water O atoms around the two topmost oxygen atoms of the EG3OMe chains (Figure 7, solid lines). One can see that not only the O4 but also O3 atoms of the EG3-OMe chains form


Langmuir, Vol. 16, No. 23, 2000 8837

Figure 7. Radial distribution of water oxygen atoms around (a) O3 and (b) O4 atoms of the EG3-OMe chains for the h- and t-SAM/water interfaces.

hydrogen bonds with water molecules. The location of minima in the distribution curves provides a convenient criterion to distinguish water molecules involved in hydrogen bonds. Assuming that each oxygen atom of the EG3-OMe chain can form no more than one hydrogen bond with water, the percentage of O3 and O4 atoms involved in hydrogen bonding with water can be calculated by integrating properly weighted Fi(r). The integration results are 75% for O4 and 2% for O3. The fact of formation of hydrogen bonds between the water molecules and the topmost oxygen atoms of the OEG chains is also in agreement with the SFG measurements,5 which show a pronounced blue shift of the methyl end group modes after immersion of the SAM in water. A similar shift has been observed in methoxyoctadecanethiol30 and associated with hydrogen bonding of the water molecules with the oxygen atoms of the methoxy group. In contrast to water at a hard wall, where the molecular dipole moments have a slightly preferred orientation away from the surface (Figure 3), the water molecules in direct contact with the h-SAM have their dipole moments directed preferentially into the SAM (Figure 8). By “direct contact”, we mean the overlapping part of the density profiles in Figure 6, which corresponds to molecules penetrated into the topmost surface layers of the SAM. Most of these molecules have their hydrogen atoms involved in hydrogen bonds with the underlying O4 and O3 atoms of the EG3-OMe chains, and as a consequence, the molecular dipole moments have a preferentially downward orientation (S01 < 0 throughout the penetration region). As can be seen from Figure 6, the profiles of S01 and S02 in the penetration region have a fine structure. In the outer part of the region (18.5-20.5 Å), where the water molecules dominate and the SAM is represented by individual protruded OEG tails, S02 is slightly negative, (30) Ong, T. H.; Davies, P. D.; Bain, C. D. Langmuir 1993, 9, 1836.

and hence the molecular dipole moments have a noticeable component perpendicular to the z axis. This component seems to be mostly due to the molecules that are not involved in hydrogen bonding with the underlying OEG chains and are more or less free in the choice of their orientation. In the inner part of the penetration region, where most (if not all) of the water molecules are hydrogen bonded to the OEG tails, S01 is large negative and S02 is large positive; i.e., the molecular dipole moments strongly prefer the downward orientation. Outside the penetration region (z > 20.5 Å), the orientational effect of the SAM is very small and does not extend farther than one-two molecular diameters away from the SAM surface. In the penetration region, the number of water molecules involved in three hydrogen bonds, h3(z), increases with decreasing z (Figure 9) and does not drop, as occurs near a hard wall. The value of 〈n〉 in the slices within the penetration region proves to be about 3, fairly close to 〈n〉 for bulk water (3.5). That is, water molecules tend to maintain a hydrogen bonding network comparable to that in the bulk. Like the orientational effect of the SAM, its effect on the “outer” part of hn(z) is restricted by one to two molecular diameters from the SAM surface. The effect of the SAM on the pair distribution function of water can be appreciated from Figure 10, which compares four distribution curves Fz(r) obtained by averaging F(r,z) over several 0.5 Å thick slices selected in four different parts of the system: (1) in four slices within the penetration range; (2) in seven slices adjacent to the SAM and corresponding to the major maximum of the water density profile d(z) in Figure 6; (3) in the next five slices; (4) in the rest of water, where the density profile d(z) shows more or less uniform oscillations. Compared to the bulk region, the main maximum of Fz(r) in the penetration range is nearly twice as high, whereas the level of Fz(r) at larger separations is noticeably lowered. Such a distortion of the pair distribution function suggests that the water molecules in the top layers of the h-SAM

8838


Pertsin and Grunze

Figure 8. Orientational order parameters of water in the h-SAM/water interphase region. The position of the SAM surface (the dash-and-dot line) and is taken to be at z where the SAM density profile vanishes.

Figure 9. Distribution of the number of hydrogen bonds per water molecule in the h-SAM/water interphase region, as a function of the distance from the substrate.

tend to form clusters, such that their constituent molecules are rich in the nearest and short in more distant neighbors. The observed tendency for clustering seems to originate from the fact that the penetration of a water molecule into the (initially ordered) SAM leads to substantial local disordering in the SAM structure. As a result, it is easier for a next water molecule to penetrate into the SAM in the neighborhood of the already penetrated molecule. The pair distribution of water molecules averaged over seven 0.5 Å thick slices adjacent to the SAM differs from that for the bulk water only slightly, while the average over the next five slices, up to two molecular diameters away from the SAM surface, is practically coincident with

the bulk distribution. That is, in our simulation model the effect of the h-SAM on the structure of water turns out to be short range. Further evidence is provided by Figure 11, which shows the lateral distribution of water density in a slice corresponding to the main peak of the water density profile. No signs of a SAM-induced lateral order can be found in this density distribution. This result is opposite to what was found in the rigid SAM model,4 where the SAM induced a clear preference of water-water separations close to the SAM lattice period (5 Å). All-Trans Modification. Compared to the helical modification of the SAM, the all-trans modification is much more resistant to disordering with water molecules. This


Langmuir, Vol. 16, No. 23, 2000 8839

Figure 10. Pair distribution functions of water averaged over z in four different parts of the h-SAM/water interphase region.

Figure 11. Lateral distribution of water density in a slice corresponding to the main peak of the water density profile d(z) in Figure 6. The minimum and maximum densities are shown in dark blue and red, respectively.

can well be seen from Figures 12 and 13, which present, respectively, the distribution of the dihedral angles, pk(τ), and a snapshot of a part of the system near the interface. The higher resistivity of the t-SAM to disordering is not unexpected in view of its noticeably higher lateral density compared to the h-SAM (the difference in lateral density being ∼15%).3 Of the 10 dihedral angles specifying the conformation of the EG3-OMe chains, only the outermost O-C-C-O and C-C-O-C dihedrals (angles 9 and 10 in Figure 12) show perceptible populations of “wrong” conformers (20 and 8%, respectively).

The density profiles of water and non-hydrogen atoms of the t-SAM (Figure 6, dashed lines) overlap to a much lesser degree, as compared to the h-SAM, thus demonstrating a less deep and intense penetration of water into the monolayer. The average density of water in the penetration range is only 9% bulk density, to be compared with 23% bulk density for the h-SAM. Unlike the h-SAM, the surface of the t-SAM is well-defined. This can be appreciated both from the snapshot in Figure 13 and from a much sharper edge of the SAM density profile in Figure 6 (dashed lines). A consequence of this is that the t-SAM

8840


Pertsin and Grunze

Figure 12. Distribution of the dihedral angles in the EG3-OMe chains for the t-SAM in contact with water. (The numbering of the angles is the same as that described in the caption to Figure 4.)

Figure 13. A snapshot of a fragment of the t-SAM/water interphase region close to the interface.

has a stronger layering effect on the density of the nearby water region: The amplitudes of the three oscillations of d(z) next to the t-SAM surface are noticeably larger than those for the h-SAM (Figure 6). The water molecules penetrated into the near-surface region of the t-SAM interact with the topmost oxygen atoms (see Figure 7, dashed lines), though the percentage of such atoms (O4) which are involved in hydrogen bonds with water is only 31% (versus 75% in the h-SAM). This is at variance with the rigid model results,4 which predict that the t-SAM cannot form O-H‚‚‚O hydrogen bonds with water at all.

Unlike the h-SAM, no hydrogen bonds are formed between the water molecules and the next to the topmost oxygen atoms of the t-SAM. Despite these differences, the average number of hydrogen bonds per water molecule was found to be again close to 3 and the molecular dipole moments to have a preferentially downward orientation, as observed in the penetration region of the h-SAM/water interface. Also similar was the effect of the t-SAM on the pair distribution functions of water molecules, except that no tendency for clustering of water molecules in the penetration region was observed.


4. Conclusions The simulation of the behavior of water near the helical and all-trans modifications of the EG3-OMe terminated alkanethiol SAM has revealed both similarities and differences. In both cases, the effect of the SAM on the structure of the adjacent water layers turns out to be short range. This refers to the pair distribution function, the orientational distribution of the molecular dipole moments, and the distribution of the number of hydrogen bonds per water molecule, which all attain their bulk behavior one to two molecular diameters away from the surface. A somewhat longer range effect is observed in the water density profiles, which exhibit well-defined SAM-induced oscillations extending three to four molecular diameters into the water bulk. Four principal differences between the h- and t-SAM/ water interfaces have been found: First, due to its lower areal density, the h-SAM is much more penetrable to water molecules. Compared to the t-SAM, the water molecules penetrate into the h-SAM deeper and in noticeably larger amounts to form as if a microscopically thin overlayer of aqueous OEG solution at the SAM surface. Second, a distinguishing feature of the h-SAM is that its interaction with water leads to a substantial conformational disordering (amorphization) of the SAM, in agreement with SFG measurements.5 The disordering affects even inner dihedrals of the OEG tails, close to their junction to the alkanethiol chains. By contrast, the

Langmuir, Vol. 16, No. 23, 2000 8841

conformation of the OEG tails in the t-SAM is affected only marginally, mostly within the terminal methoxy groups. Third, due to the conformational disordering of the OEG tails, the surface of the h-SAM is ill-defined. As a consequence, the h-SAM has a noticeably weaker effect on the z distribution of water density near the SAM surface, compared to the t-SAM. Fourth, an important property that distinguishes the helical modification of the SAM from the all-trans one is a much higher surface density of hydrogen bonds with water molecules and, hence, much higher hydrophilicity. Of the four above differences between the h- and t-SAMs, only the last one has a direct bearing to protein resistance: Since the adsorption of a protein molecule on the SAM surface involves dehydration of the OEG tails, the higher the surface density of hydrogen bonds, the more resistant to protein adsorption should be the surface. The role of the three other distinctive features of the h-SAM in its protein resistance is not so evident and still remains to be understood. Acknowledgment. We thank the Deutsche Forschungsgemeinschaft and the Fond der Chemischen Industrie for the financial support of this work. Stimulating discussions with H. J. Kreuzer and M. Buck are gratefully acknowledged. LA000340Y

Computer Simulation of Water near the Surface of Oligo(ethylene

Recommend Documents