J. Phys. Chem. B 2008, 112, 14221–14229
14221
Molecular Dynamics Description of Grafted Monolayers: Effect of the Surface Coverage F. Goujon,† C. Bonal,*,† B. Limoges,‡ and P. Malfreyt† Laboratoire de Thermodynamique et Interactions Mole´culaires de l’UniVersite´ Blaise Pascal (Clermont-Ferrand II), FRE 3099, 24 AVenue des Landais, 63177 Aubie`re Cedex, and Laboratoire d’Electrochimie Mole´culaire de l’UniVersite´ Denis Diderot (Paris 7), UMR CNRS 7591, 2 Place Jussieu, 75251 Paris Cedex 05, France ReceiVed: April 3, 2008; ReVised Manuscript ReceiVed: August 28, 2008
Molecular dynamics simulations of monolayers of metal-chelating ligands grafted onto a graphite surface in water are carried out to calculate structural (density profiles, radius of gyration, and asphericity coefficients), dynamical (diffusion coefficients), and energetical properties as a function of the surface coverage. The purpose is to provide a better understanding of the dependence of various properties of these monolayers on the surface coverage. A critical value of the surface coverage from which all structural properties derive a limiting value has been established. It also appears that the chains rather adopt an elongated conformation along the direction normal to the surface from this critical surface coverage. The hydrogen-bonding structure and dynamics of water molecules are reported. An ordered structure of water in the region close to the terminal groups of the grafted molecules is shown at a relatively high surface coverage. This ordering is similar to that observed in the case of water in interaction with a solid surface. 1. Introduction Controlling the molecular structure and organization of macrobiomolecules such as proteins immobilized on an electrode surface without distressing their native functionality is a key issue in the development of applications such as electrochemical biosensors,1 biomolecular electronic devices,2 or biofuel cells.3 An essential step for protein immobilization is to construct a well-defined surface exhibiting selective affinity for the desired protein. One attractive method is based on the formation of ternary metal-chelate complexes between a chelating agent anchored on a conductive surface and a histidine-tagged recombinant protein. Such a strategy was proved useful for controlling protein orientation at interfaces and thus improving the recognition between a receptor immobilized on a surface and a ligand in solution.4 It also allows modulation of the electrical communication between an electrode and a redox protein site specifically oriented on its surface.2,5 Self-assembled monolayers (SAMs) of nitrilotriacetic (NTA)-terminated alkanethiol on a gold surface, loaded with divalent transition metal ions such as Cu(II) or Ni(II), have been by far the most investigated platforms in the applications mentioned above.4-6 However, the use of gold electrodes has a number of disadvantages such as a narrow potential window and a slow degradation of the chemisorbed monolayer upon storage.7 In addition, the details of the interaction between the thiol headgroup and the gold surface remain uncertain to some extent, and the bond is not very strong.8 Carbon electrodes are an attractive alternative to gold since they offer a wider potential window and their surface can be functionalized by a true covalent linkage of the metal-chelating species. Such a possibility was recently illustrated by Blankespoor et al.,9 who developed a convenient method to electrochemically address the functionalization of a carbon surface with a ligand. The * To whom correspondence should be addressed. E-mail: Christine.Bonal@ univ-bpclermont.fr. † Universite ´ Blaise Pascal (Clermont-Ferrand II). ‡ Universite ´ Denis Diderot (Paris 7).
method takes advantage of the electrochemical reduction of an aryl diazonium salt substituted with a NTA or an imminodiacetic (IDA) ligand at the para position. The transient aryl radical was generated during the reduction; it thus forms a covalent carbon-carbon bond with the electrode surface, producing a densely packed monolayer of the ligand. The NTA- and IDAmodified electrodes were then shown to efficiently chelate Cu(II) and Ni(II) ions and, once loaded with a metal ion, were able to specifically bind histidine-tagged proteins (see Figure 1). Optimization of the interactions between immobilized molecules and their corresponding target molecules in solution requires a better understanding of the properties governing these monolayers. Earlier theoretical and experimental studies10,11 have already suggested that several factors such as the packing density, the chain lengths of the monolayers, or the chemical composition of the monolayers can modify the interfacial structure and dynamics. For instance, with a standard surface of oligo(ethylene oxide) SAMs that resists protein adsorption, experimental measurements10 have shown that the protein resistance depends on the packing density. Before a more sophisticated investigation of the interaction between the immobilized molecules and their corresponding target molecules in solution, it is thus of primary importance to focus on the impact of the surface coverage. In fact, it appears that the thermodynamic of the system is a key issue that has not been addressed in past studies. The structural organization of molecules assembling on a surface is strongly affected by the various energy contributions within the system.12 In addition, in spite of a large panel of surface analysis techniques, many questions still remain unanswered. For instance, a quantitative understanding of the details of molecularscale surface architecture is lacking. Computational chemistry, which simulates chemical structures and reactions numerically, is now a vital adjunct to experimental studies.13 Computer simulation also provides a natural complementary method to answer some of the questions about the general properties governing these systems and has already made
10.1021/jp8028825 CCC: $40.75 2008 American Chemical Society Published on Web 10/18/2008
14222 J. Phys. Chem. B, Vol. 112, No. 45, 2008
Goujon et al. chemical studies of monolayers on others surfaces than gold have been reported,30-34 and none of them has focused on negatively charged monolayers of chelating metal compounds immobilized onto the graphite surface via the tetradentate NTA acid (Figure 1). The aim of this work is to investigate in detail these grafted systems by molecular simulation. Recently,35 a methodology for the molecular dynamics simulations of systems with a finite length in the third dimension has been established. Our main objective is to use this methodology to perform a series of comparisons on the structure, dynamics, and energetic properties as a function of surface coverage. Structural properties such as density profiles, radius of gyration, and asphericity coefficients are investigated. The dynamic properties of the water molecules through transport properties such as diffusion coefficients and the hydrogen-bonding structure are also reported. To conclude, the energetics of these systems is studied as a function of surface coverage. The outline of the paper is as follows. In the section concerning the experimental methods, we give the details of the computational models, while in the next section, we present the results of the simulation and compare the data at low and high surface coverage. Finally, in the last section, we draw the main conclusions from this work. 2. Computational Procedures The phenyl-NTA molecule was modeled using the all-atom (AA) version of the Cornell force field AMBER.36 The general potential function is of the form
U)
∑ kb(r - r0)2 + ∑ kθ(θ - θ0)2+ ∑
bonds
angles
N-1
cos(lφ + δ)] +
kφ[1 +
dihedrals
N
∑∑
i)1 j)i+1
{ [( ) ( ) ] 4εij
σij rij
12
∞
σij rij
6
qq
+
j ∑ ’ |rij +i nL| l
Figure 1. (a) NTA-Cu(II) complexes grafted onto the solid surface. (b) Snapshot of the MD configuration of NTA-Cu(II) complexes grafted onto the HOPG surface without the water molecules at 55 pmol cm-2 (oxygen, red; carbon, black; nitrogen, blue; copper, yellow; and hydrogen, white). (c) Snapshot of the MD configuration of NTA-Cu(II) complexes grafted onto the HOPG surface without the water molecules at 297 pmol cm-2 (oxygen, red; carbon, black; nitrogen, blue; copper, yellow; and hydrogen, white).
an important contribution to the understanding of the processes involved. Following the pioneering work of Klein’s group,14-16 a number of simulations employed force field to explore the structure and dynamics of SAMs. However, the majority of the work was limited to the study of alkanethiols of different tail groups, self-assembled on the Au(III) surface in the presence of water.14-24 Several molecular dynamics and Monte Carlo computer simulation studies of carboxylic acid-functionalized SAMs25-27 and of oligo(ethylene oxide) or poly(ethylene oxide) SAMs have recently been reported.28,29 However, few theoretical
}
(1)
where kb, kθ, and kΦ are the force constants for the deformation of bonds, angles, and dihedrals, respectively. The equilibrium values of bond distances and valence angles correspond to r0 and θ0, respectively. In the dihedral angle term, l is the periodicity, and δ is the phase factor. The C-H and O-H covalent bonds were kept at fixed lengths with the use of the SHAKE algorithm,37 and the aromatic rings were kept planar by using six improper torsional potentials. The intermolecular and intramolecular interactions consisted of a van der Waals repulsion-dispersion term calculated from the Lennard-Jones (6-12) potential, represented by the penultimate term in eq 1. In the AMBER force field, the nonbonded interactions between atoms separated by exactly three bonds (van der Waals 1-4 interactions) were reduced by a factor of 0.5.36 The Lennard-Jones potential parameters for the interactions between unlike atoms were calculated by using the Lorentz-Berthelot mixing rules quadratic and arithmetic rules for εij and σij parameters, respectively. The water molecules were represented by the TIP4P/2005 model.38 This representation consists of four representation sites: Three of them are placed at the oxygen and hydrogen atom positions, respectively. The fourth site is placed at the bisector of the H-O-H angle. The metal-ligand bond was treated as a covalent one, and it was modeled using the following Morse potential:39
Molecular Dynamics Description of Grafted Monolayers
J. Phys. Chem. B, Vol. 112, No. 45, 2008 14223
Umetal · · · ligand ) E0[{1 - exp[-k(rij - r0)]}2 - 1] where E0 is equal to 359 and 918 kJ mol-1 and k is 1.85 and 1.2 Å-1 for Cu-O and Cu-N bonds, respectively. The last term in eq 1 corresponds to the electrostatic energy (UELEC) of the system. Consider N molecules in a volume V ) LxLyLz with a center of mass ri; each molecule contains ni charges qia at position ria. The electrostatic interactions, handled with the Ewald sum method40 in a box with an orthogonal axis that has an overall neutral charge, are given by the following contributions:
UELEC )
∑
1 Q(h)S(h)S(-h) + 2ε0V k*0 1 8πε0
R 3⁄ 2
4π ε0
∑ ∑ ∑ qia∑ qjberfc(Rriajb ⁄ riajb) i
a
j*i
b
qiaqib erf(Rriaib) + r b*a iaib
1 ∑ ∑ qia2 - 8πε ∑∑∑ 0 i
a
i
a
1 M 2 (2) 2ε0V z where erfc(x) is the complementary error function and erf(x) is the error function. R was chosen so that only pair interactions in the central cell need to be considered in evaluating the second term in eq 2. In the last term of eq 2, Mz is the net dipole moment of the simulation cell given by ∑Ni)1qiri. This is the contribution term from Yeh and Berkowitz,41 which results by the planewise summation method proposed by Smith.40 Adding this term to the total energy amounts to using a z-component force for each atom, which is given by
Fi,z ) -
qi M ε0V z
(3)
The EW3DC method only differs from the standard EW3D method by the presence of this dipole correction. The functions S(h) and Q(h) are defined respectively by:
S(h) )
∑ ∑ qiaexp(ih · ria) i
Q(h) )
a
( )
1 h2 exp - 2 2 h 4R
(4) (5)
The reciprocal lattice vector h is defined as h ) 2π(l/Lx, m/Ly, n/Lz), where l, m, and n can take values of 0, (1, (2,..., (∞. The atomic charges of NTA-Cu(II) complexes grafted onto solid surface were calculated from the density functional theory (DFT)42,43 (B3LYP)44-46 with effective core potential (SDDALL) Gaussian basis sets using the Gaussian 03 package47 and the CHELPG48 procedure as a grid-based method. The equations of motions were integrated using the Verlet Leapfrog algorithm scheme at T ) 298 K with a time step of 2 fs. The simulations were performed in the NVT ensemble using the Hoover algorithm49 with a coupling constant of 0.5 ps. The SHAKE algorithm was used to preserve the rigidity of the bonds involving hydrogen atoms so that a higher time step value can be used. The cutoff for the Lennard-Jones contributions was fixed to 12 Å, whereas it was equal to 18 Å for the real space of the electrostatic interactions. We also used a reciprocal space cutoff radius of 1.14 Å-1. The configurations were generated by the parallel version of the modified DL_POLY_MD package50 by using up to eight processors at a time.
TABLE 1: Number of Water Molecules as a Function of the Grafting Density no. of molecules grafted
no. of water molecules
surface coverage calculated (pmol cm-2)
1 5 10 14 18 23 27 36
3835 3775 3715 3659 3603 3533 3477 3351
011 055 110 154 198 253 297 396
The surface consisted of a single layer of aromatic carbon atoms modeling a surface of highly oriented pyrolitic graphite (HOPG). The molecules were chemically grafted by a covalent bond between the aryl group and one carbon of the graphite surface. Depending on the surface coverage, these grafting points were chosen at random using an additional distance criterion to avoid overlaps between neighbor molecules. The size of the surface was 16 × 18 unit cells, which corresponded to 576 carbon atoms, leading to a box dimension of 39.4 × 38.3 Å2 for the xy plane. The dimension along the z-axis was set to 80.0 Å so that the water was large enough to exhibit a bulk behavior.35 However, as we considered this system as 2Dperiodic, the system was surrounded by a vacuum along z-axis, increasing the box dimension to 320 Å. The inclusion of empty space into the central box was carried out to dampen out interslab interactions and to avoid consequently an artificial influence from the periodic images in the z direction. Considering this extended system, the number of k-vectors used for solving the electrostatic interactions with the EW3DC method was 7, 7, and 58 for the x, y, and z directions, respectively. The parameters of the long-range electrostatic interactions evaluated by the Ewald summation technique were calculated to satisfy a relative error of 10-6. This relative error in the reciprocal term 2 2 was approximately expressed as exp(-hmax /4R2)/hmax , where hmax is calculated from the largest h-vector in reciprocal space. This gave a value of R ) 0.17329 Å-1 for the Ewald convergence parameter. We checked that the Ewald sum was correctly converged in the reciprocal space by comparing the total electrostatic energy and the electrostatic atomic virial. The deviation between these two energy contributions did not exceed 0.003% for our system sizes, as expected for a good convergence of the Ewald sum. Because the simulations were carried out at constant volume, the number of water molecules was adjusted to get a normal bulk density in the upper part of the box. For all surface coverage studied, the number of molecules grafted on the surface and the number of water molecules are listed in Table 1. One simulation consisted of an equilibration period of 500 ps during which the velocities were rescaled every 100 steps. Data were then collected during the acquisition period of 1 ns, and one configuration was stored every 0.1 ps. This led to 5000 saved configurations from which the energetic and structural properties were averaged. 3. Results and Discussion For grafted systems, the density profiles along the direction normal to the surface are a basic characterization of the local structure. Figure 2a shows the molecular density profiles of the grafted NTA-modified complexes as a function of the surface coverage. The first five well-defined peaks are easily attributed to the carbon atoms of the phenyl group covalently attached to the graphite surface. The next two peaks correspond to the
14224 J. Phys. Chem. B, Vol. 112, No. 45, 2008
Figure 2. (a) Atomic density profiles of the grafted molecules along the z direction for different values of surface coverage. For clarity, density profiles at various grafting densities are offset by 1 unit. (b) Water molecular density profiles along the z direction for different surface coverages. (c) Average number of water molecules along the z direction for different surface coverages.
carbon atoms close to the amide function (Figure 1). At higher surface coverages, there is a significant increase in the magnitude of these peaks due to the increase of the number of grafted molecules. However, the position of the peaks is not affected. At lower surface coverage, these well-marked peaks are followed by weak oscillations in the range of 8 Å < z < 18 Å. They correspond to the terminal group of the grafted molecules and indicate no privileged orientation. From a surface coverage of 154 pmol cm-2, these weak oscillations are replaced by peaks. This means that the orientation of the terminal group is less free when the number of the grafted molecules increases. Nevertheless, this profile highlights that the first part of the grafted molecule has a preferential orientation along the direction normal to the surface due to the covalent bond between the aryl group and the graphite surface. This contrasts with the behavior of the terminal group, which assumes a more random orientation. Figure 2b depicts the density profiles of the water molecules along the z direction. We first observe two distinct water layers propagating from the solid surfaces. The resulting peaks have already been observed in water molecules in interaction with metal surfaces and hydrophobic walls.51 In agreement with experimental results,52 computer simulations have shown that the interaction of water with a solid surface induces layering in the confined liquid and consequently changes its structural and dynamical properties.53,54 In addition, strong oscillations of the water density at the surface, lateral positional, and orientational ordering of the absorbed water molecules have been reported. A zone of interfacial water molecules characterized by reduced mobility as compared to those within the liquid bulk has also
Goujon et al. been observed. After these two main peaks, oscillations are observed in the range of 8 Å < z < 14 Å, showing less distinct layers and a value of the water density lower than that of the bulk density. At approximately 14 Å to the surface, the water density increases again until it reaches its bulk density at z ) 25 Å, as it was shown in our previous paper.35 This behavior is clearly related to the structure of the grafted NTA-Cu (II) complexes. In fact, it appears that the grafted molecules have a hydrophobic headgroup (corresponding to the region where the water density decreases with the surface coverage) and a large hydrophilic charged tail group (corresponding to the region where the water density increases with the surface coverage). As the surface coverage increases, we observe in the range of 8 Å < z < 14 Å that the water density quickly drops. This involves less and less penetration of the water into the region surrounding the hydrophobic headgroup of the grafted molecule. As expected in the case of the dense packing of the chains, the free volume is weak and hinders the diffusion of the water molecule in the monolayers. This is confirmed in Figure 2c, where the average number of water molecules is calculated as a function of z for different surface coverages. Indeed, at a distance of 14 Å from the surface, the number of water molecules is equal to 600 for a surface coverage of 11 pmol cm-2 and decreases to 200 for a surface coverage of 396 pmol cm-2. The emergence of two peaks located at 18 and 22 Å from a surface coverage of 110 pmol cm2 is also observed in Figure 2b. The height of these peaks increases with the surface coverage. This leads to two well-defined peaks at the highest coverage. These latters correspond to two distinct water layers into the region of the large hydrophilic charged tail group. This result indicates that the dense packing of the negatively charged monolayers creates a wall that induces layering in the water molecules. We then observe an ordering of water identical to the one close to the solid surfaces. Figure 3 represents the two-dimensional density profiles of water contained in the box slab from z ) 15 Å to z ) 18 Å at low (55.0 pmol cm-2) and high (297 pmol cm-2) surface coverages, respectively. As the yellow color corresponds to the average bulk density, we can observe (Figure 3a) that for a low surface coverage, the water molecules are bulklike, showing a homogeneous structure. The low-density zones represent the regions occupied by the grafted molecules. Figure 3b shows the equivalent density profile at a high surface coverage. Because of the large number of grafted molecules, less water molecules are present in the grafting regions. The grafted molecules cannot move a lot around the grafting point because of the steric repulsion between molecules. Besides, we can observe small sharp peaks for water density with an average value greater than the bulk density, especially in small water zones joining fixed grafted molecules. Such results lead to two hypotheses: The first one is that the grafted molecules may have a fixed geometry at high surface coverage; this will be checked in the next section by calculating structural quantities. The second one concerns the water molecules: Hydrogen bonding with grafted molecules may constrain water molecules at some places, leading to a decrease of the water diffusion in the grafted region. To obtain a quantitative estimate of the structural properties such as the size of the molecule, we computed the mean square radius of gyration given by:
〈RG2 〉 )
1 MT
〈
N
∑ mi|ri - rcm|2 i)1
〉
(6)
where MT and N are the total mass and total number of atoms, respectively, mi and ri are the mass and position of the ith atom,
Molecular Dynamics Description of Grafted Monolayers
J. Phys. Chem. B, Vol. 112, No. 45, 2008 14225 A useful descriptor of anisometry is the asphericity55 order parameter, A, defined as
1 A) 2
〈∑ (λ - λ ) 〉 i
j
〈(∑ ) 〉
i>j
3
2
2
(8)
λi
i)1
A measures the deviation from a spherical form. If A tends to zero, the molecule adopts a spherical conformation, while A tends to 1 as it becomes a thin rod. The values of the mean square radius of gyration as a function of surface coverage are plotted in Figure 4a with the components of RG2 along the x, y, and z directions. At the lowest surface coverage (11 pmol cm-2), the value of RG2 (17.5 ( 2.7 Å2) is almost equal to the one calculated in the homogeneous system constituted by only the NTA-modified molecule chelated by Cu(II) in water (17.6 ( 3.4 Å2). At low surface coverage, the three components of RG2 are similar, whereas the x and y components are smaller than the z component at higher surface coverages. This indicates that the conformation of the molecule changes with the surface coverage; an elongation of the chain along the direction normal to the surface is observed. Moreover, these results not only show this expansion of RG2 in the z direction but also that the main square radius of gyration rapidly reaches a limiting value of 28 Å2. This value is attained from a surface coverage of 154 pmol cm-2.
Figure 3. Two-dimensional water density profiles calculated in the zone 15 Å < z < 18 Å for two values of surface coverages: (a) 55 and (b) 297 pmol cm-2.
and rcm refers to the position of the center of mass. In fact, the radius of gyration of a group of atoms is defined as the rootmean-square distance from each atom of the molecule to their centroid. The shape of the system can also be estimated from the inertia tensor S. Diagonalization of S results in three eigenvalues λ, the sum of which is the mean-squared radius of gyration and the largest of which corresponds to an eigenvalue vector representing the long axis of the molecule. The eigenvalues are related to the radius of gyration from the following equation: 2 〈RG2 〉 ) R2Gx + R2Gy + RGz ) λ1 + λ2 + λ3
(7)
The three eigenvalues of the inertia tensor define the lengths of the principal axes of the molecule and provide information about the absolute shape of the molecule. These eigenvalues cannot be related to the components of the square radius of gyration that indicate a relative shape of the molecule within the Cartesian coordinate frame.
Figure 4. (a) Mean squared radius of gyration. (b) Eigenvalues of the inertia tensor. (c) Asphericity coefficients as a function of the surface coverage.
14226 J. Phys. Chem. B, Vol. 112, No. 45, 2008
Figure 5. Hydrogen bond densities along the z direction for a surface coverage of 297 pmol cm-2.
The plot of the three eigenvalues of the inertia tensor as a function of grafting density (Figure 4b) presents a profile similar to the one obtained for the main square radius of gyration. It is also possible to check the changes in the molecular geometry of the grafted molecules from the estimation of the tilt angle, defined as the angle between λ1 and the z direction. From a 2 surface coverage of 154 pmol cm2, λ1 is almost equal to RGz , indicating that the tilt angle is very small. This feature is in agreement with the fact that the monolayers adopt an elongated conformation normal to the surface. Whereas Figure 4a shows the equivalence between the x and y components of the square radius of gyration from 110 pmol cm-2, Figure 4b establishes two distinct values for the second and third eigenvalues λ2 and λ3. This comparison indicates that the absolute shape of the grafted molecules presents three different lengths for the principal axes, but the statistical distribution of the shape is identical along the x and y directions as expected for an equilibrium conformation with no external constraint in these two directions. The calculation of the asphericity coefficient gives a value of 0.24 ( 0.07 at the lowest surface coverage (Figure 4c). The value of the asphericity coefficient at the lowest surface coverage is found to be equal to that calculated in the homogeneous system (0.24 ( 0.07). From 154 pmol cm-2, the asphericity coefficient reaches a limiting value of 0.495. It is interesting to associate the description of the conformation of the grafted monolayers with those of grafted polymer chains by using the “mushroom” and “brush” regimes.56 We are aware that the different lengths and time scales spanned by the polymer chains are totally different from those of the grafted molecules studied here. However, it appears that at the lowest surface coverage (11 pmol cm-2), the conformation of our monolayer can be compared with the “mushroom” regime where the chains are isolated from each other and form “mushrooms” on the surface. For surface coverage varying from 55 to 110 pmol cm-2, the concentration of grafted molecules (brush concentration) allows us to observe the scaling of some structural properties as a function of the surface coverage. From 154 pmol cm-2, the structure still remains a brush type with chains that stretch away from the interface to avoid overlapping, but no scaling law can be observed. Let us specify that it is not possible to reach this surface coverage in the case of grafting polymer chains. The hydrogen-bonding structure and dynamics for water molecules in the system are also studied. Hydrogen bonds can be defined on the basis of either energetic or geometric criteria. Hydrogen bonds occur if the donor-acceptor distance is lower than 2.5 Å and the hydrogen donor-acceptor angle is smaller than 180° ( 30°. Figure 5 shows the density profiles of
Goujon et al. hydrogen bonds along the direction normal to the surface at the highest surface coverage (396 pmol cm-2). A detailed analysis of the hydrogen-bonding structure between atoms of the grafted molecule and water molecules reveals two pronounced peaks at 7 and 17 Å. The first peak corresponds to the strong hydrogen bonds between the atoms of the amide group and the water molecules, and the second peak corresponds to the hydrogen bonds between the oxygen atoms of COOheadgroup and the water molecules. As expected from the density profile of the water molecules (Figure 2b), the hydrogen bond density in the range of 8 Å < z < 14 Å is weak because this region corresponds to the hydrophobic group of the grafted molecules. Moreover, at higher surface coverages, the water penetration into the layers is hindered. It should also be noted that in the vicinity of the COO- headgroup (17 Å), the total number of hydrogen bonds is mainly due to the interaction between the water and the grafted molecules, whereas close to the solid surface, the hydrogen-bonding structure is only due to hydrogen bonds between water molecules. To examine in detail the dynamic behavior of hydration water, the self-diffusion coefficients were also calculated. Self-diffusion coefficients provide insight into the water dynamics and mobility near the surfaces. The translational self-diffusion coefficients were determined by two methods: from mean squared displacement of the molecular centers of mass using the Einstein relation
Dmsd ) lim tf∞
1 6NN0t
〈
N0
N
∑ ∑ [rcom,i(t + t0) - rcom,i(t0)]2
t0)1 i)1
〉
(9)
or from the integration of the velocity autocorrelation function57
Dvacf )
1 3NN0
〈
N0
N
∫0 ∑ ∑ [Vcom,i(t + t0) · Vcom,i(t0)] ∞
t0)1 i)1
〉
dt
(10) where rcom,i and Vcom,i are the position and the velocity of the center of mass of water molecule, respectively. The water diffusion coefficient was calculated in the grafted zone (0 Å < z < 20 Å) and in the water zone (35 Å < z < 80 Å). We took into account only the water molecules that never left the region of interest during the calculation. The mean square displacement and the velocity autocorrelation functions were calculated over a trajectory of 50 ps using N0 ) 60 time origins. The diffusion coefficient calculated from eq 9 is reported in Figure 6a. The values of the diffusion coefficient calculated from eq 10 are in line with those calculated from the mean square displacement within the statistical fluctuations. The statistical fluctuation were calculated using the block average method. The maximum standard deviation is about 0.1 × 10-9 m2 s-1. The error bars are not shown in Figure 6a for clarity. As expected, the parallel and perpendicular components of the diffusion coefficient in the water zone are not dependent on the surface coverage. Small differences between the parallel and the perpendicular components are observed. The value of the lateral component of the diffusion coefficient is approximately equal to that of the isotropic diffusion coefficient for the TIP4P/2005 water model at 298 K (2.4 × 10-9 m2 s-1), whereas the value of the perpendicular component is slightly weaker. The fact that the bulk diffusive behavior of water molecules is not completely recovered and retarded in the direction perpendicular to the surface indicates that the region between 35 and 80 Å is weakly affected by the presence of the graphite surface. The values of the components of the diffusion coefficient of water molecules embedded in the modified surface as a function
Molecular Dynamics Description of Grafted Monolayers
Figure 6. (a) Diffusion coefficients of the water in the grafted zone and in the bulk zone as a function of the grafting density. (b) Velocity autocorrelation functions as a function of time. The velocity autocorrelation function of bulk water molecules is given for comparison. (c) Mean square displacement as a function of time.
of the surface coverage were calculated. First, it appears that the parallel and perpendicular components follow the same trend with respect to the grafting density, and they become equal in magnitude from 154 pmol cm2. At high grafting densities, the diffusion coefficient is similar in magnitude to that of water molecules in monohydrated montmorillonites.58 When the coverage increases, the interfacial region becomes more compact, and the water molecules are trapped between the grafted molecules. The grafted molecules behave as if they are a part of the graphite surface, significantly reducing the diffusion of trapped water molecules. This observation is in line with the ordering of water in the region close to the terminal groups of the grafted molecules. Part b of Figure 6 shows the center of mass velocity autocorrelation functions of water molecules of the grafted region against of the surface coverage. For comparison, the velocity autocorrelation function of bulk water molecules at 298 K is reported. Interestingly, significant changes between the correlations functions of water molecules in the grafted zone and in the isotropic bulk water are observed. As the surface
J. Phys. Chem. B, Vol. 112, No. 45, 2008 14227
Figure 7. (a) Dispersive and repulsive parts of the Lennard-Jones potential and electrostatic part for the molecule-molecule interaction energy as a function of the surface coverage. (b) Dispersive and repulsive parts of the Lennard-Jones potential and electrostatic part for the molecule-water interaction energy as a function of the surface coverage. For comparison, the corresponding contributions of the homogeneous system constituted by a molecule of NTA-modified complex in water are reported in dotted lines. (c) Number of interacting water molecules within two cutoff radii for different values of surface coverage. The dotted lines also give the number of interacting water molecules calculated in isotropic bulk water for these two cutoff radii.
coverage increases, the first peak of the velocity autocorrelation function becomes more negative. The negative value of this peak is attributed to the reversal of the velocities due to a cage effect. The increasing magnitude of the negative peak with the surface coverage confirms that the water molecules rattle in a decreasing accessible volume. Part c of Figure 6 presents the mean square displacement of water molecule in the grafted zone for three surface coverages and in an isotropic bulk water. The intermolecular molecule-molecule energy was decomposed into the Coulombic, dispersive, and repulsive energy contributions. We report these contributions as a function of the surface coverage in Figure 7a. The Lennard-Jones potential is composed of two terms: One term depending on the separation by 1/r12 represents the very short-range repulsion, and a second term, depending on -1/r6, is used to model the van der Waals interactions based on dipole/dipole, dipole/induced dipole, and induced dipole/induced dipole interactions (Keesom, Debye, and London dispersion energies). It appears in Figure 7a that the electrostatic and the Lennard-Jones repulsive energy contributions increase with the surface coverage, whereas the LennardJones dispersive contribution decreases. Results show that the
14228 J. Phys. Chem. B, Vol. 112, No. 45, 2008 unfavorable repulsive and electrostatic interactions due to a relatively local high density are compensated by favorable van der Waals interactions between the groups of the chains. This compensation acts at all coverages. The same energy contributions were calculated for the molecule-water pair as a function of surface coverage (Figure 7b). Significant differences appear between the contributions calculated in heterogeneous and homogeneous (dotted lines in Figure 7b) systems. As expected, the magnitude of these differences increases with the surface coverage. The electrostatic and the Lennard-Jones dispersive energies increase with the surface coverage, whereas the Lennard-Jones repulsive contribution decreases. These observations are in line with the fact that there are significantly less neighboring water molecules in interaction with the grafted molecules at high surface coverage. This is established in Figure 7c. A quasi compensation between the two Lennard-Jones contributions is obtained. The major contribution to the total molecule-water energy comes from the electrostatic energy, this energy being in relation with the polar nature of the water molecule. However, it is very difficult to extract from these simulations a quantitative energy balance between the enthalpy and the entropic contributions. Two different entropic contributions are expected as follows: a favorable term due to the consequent release of water molecules when the surface coverage increases and an unfavorable term due to the loss of conformational flexibility when the number of chains grafted increase. The entropy change resulting from these two different contributions is then difficult to evaluate from standard molecular simulations. 4. Conclusions Negatively charged monolayers immobilized onto graphite surface have been investigated using molecular simulations. The density profiles along the direction normal to the surface have been reported. The first part of the grafted molecule has a preferential orientation along the direction normal to the surface due to the covalent bond between the aryl group and the graphite surface. This contrasts with that of the terminal group, which is further random. The density profile of the water molecules shows that the dense packing of the negatively charged monolayers forms a “wall” at high surface coverages and induces a layering in the water molecules. The two-dimensional density profiles of water molecules have been reported in the lateral plane with respect to the surface coverage. At low surface coverages, the water molecules are bulklike, showing a homogeneous density distribution, whereas at higher surface coverages, distinct layers corresponding to highly ordered water are observed. A series of comparisons on the structure properties such as radius of gyration and asphericity coefficients have been carried out as a function of the grafting density. The chains adopt rather an elongated conformation along the direction normal to the surface from 154 pmol cm-2. From this value, all structural parameters tend to a limiting value independent of the surface coverage. A high degree of ordered and dense nature of SAMs is then obtained from this critical value. In fact, in an ideal closepacked monolayer, the grafted molecules are assumed to have a homogeneous density and to be tethered perpendicularly to the surface. Our results interestingly show that these factors are highly dependent on surface coverage. It should also be emphasized that many applications actually depend crucially on structural details. Within this context, fundamental questions are addressed concerning the influence of the surface coverage and the impact of the various degrees of freedom of the grafted molecules.
Goujon et al. The hydrogen-bonding structure and dynamics for water molecules in the system have also been studied. We show that the interfacial region becomes more compact and that the water molecules trapped between the grafted molecules are less mobile as the surface coverage is increased. At high surface coverages, the grafted molecules behave as if they are a part of the graphite surface, reducing significantly the diffusion of trapped water molecules. This observation is in line with the ordering of water in the region close to the terminal groups. The molecule-molecule and molecule-water contributions to the potential energy have been calculated and decomposed into electrostatic, repulsive, and dispersive energies. Our results interestingly show that the molecule-molecule total energy is not dependent on the surface coverage because there is compensation between the favorable van der Waals interactions and the unfavorable repulsive and electrostatic energy parts. However, the total effect of the surface coverage on the energetic properties is difficult to predict because the entropic contributions cannot be estimated from our simulations. References and Notes (1) Gooding, J. J.; King, G. C. J. Mater. Chem. 2005, 15, 4876. (2) Benson, D. E.; Conrad, D. W.; De Lorimier, R. M.; Trammell, S. A.; Hellinga, H. W. Science 2001, 293, 1641. (3) Willner, I.; Katz, E. Angew. Chem., Int. Ed. 2000, 39, 1181. (4) (a) Sigal, G. B.; Bamdad, C.; Barberis, A.; Strominger, J.; Whitesides, G. M. Anal. Chem. 1996, 68, 490. (b) Klenkar, G.; Valiokas, R.; Lundstroem, I.; Tinazli, A.; Tampe, R.; Piehler, J.; Liedberg, B. Anal. Chem. 2006, 78, 3643. (c) Lata, S.; Piehler, J. Anal. Chem. 2005, 77, 1096. (5) (a) Madoz-Gurpide, J.; Abad, J. M.; Fernandez-Recio, J.; Velez, M.; Vazquez, L.; Gomez-Moreno, L.; Fernandez, V. M. J. Am. Chem. Soc. 2000, 122, 9808. (b) Johnson, D. L.; Martin, L. L. J. Am. Chem. Soc. 2005, 127, 2018. (6) (a) Lee, J. K.; Kim, Y.-G.; Chi, Y. S.; Yun, W. S.; Choi, I. S. J. Phys. Chem. B 2004, 108, 7665. (b) Luk, Y.-Y.; Tingey, M. L.; Hall, D. J.; Israel, B. A.; Murphy, C. J.; Bertics, P. J.; Abott, N. L. Langmuir 2003, 19, 1671. (c) Kro¨ger, D.; Liley, M.; Schiweck, W.; Skerra, A.; Vogel, H. Biosens. Bioelectron. 1999, 14, 155. (d) Keller, T. A.; Duschl, C.; Kro¨ger, D.; Se´vin-Landais, A.-F.; Vogel, H.; Cervigni, S. E.; Dumy, P. Supramol. Sci. 1995, 2, 155. (7) Beulen, M. W.; Kastenberg, M. I.; van Veggel, C. J. M.; Reinhoudt, D. N. Langmuir 1998, 14, 7643. (8) Bernard, M.-C.; Chausse´, A.; Cabet-Deliry, E.; Chehimi, M. M.; Pinson, J.; Podvorica, F.; Vautrin-Ui, C. Chem. Mater. 2003, 15, 3450. (9) Blankespoor, R.; Limoges, B.; Scho¨llhorn, B.; Syssa-Magale´, J.L.; Yazidi, D. Langmuir 2005, 21, 3362. (10) Ismail, E. A.; Grest, G. S.; Stevens, M. J. Langmuir 2007, 23, 8508. (11) Herrwerth, S.; Eck, W.; Reinhardt, S.; Grunze, M. J. Am. Chem. Soc. 2003, 125, 9359. (12) Anderson, M. R.; Evaniak, M. N.; Zhang, M. Langmuir 1996, 12, 2327. (13) Leach, A. R. Molecular Modelling Principles and Applications; Pearson Education EMA: Harlow and New York, 2001. (14) Hautman, J.; Klein, M. L. J. Chem. Phys. 1989, 91, 4994. (15) Hautman, J.; Klein, M. L. J. Chem. Phys. 1990, 93, 7483. (16) Mar, W.; Klein, M. L. Langmuir 1994, 10, 188–196. (17) Zhang, J.; Bilic, A.; Reimers, R. J.; Hush, S. N.; Ulstrup, J. J. Phys. Chem. B 2005, 109, 15355. (18) Srivastava, P.; Chapman, W. G.; Laibinis, E. P. Langmuir 2005, 12171. (19) Kim, Y.-H.; Jang, S. S.; Goddard, W. A. J. Chem. Phys. 2005, 122, 244703. (20) Rousseau, R.; Mazarello, R.; Scandolo, S. ChemPhysChem 2005, 6, 1756. (21) Zheng, J.; He, Y.; Chen, S.; Li, L.; Bernards, M. T.; Jiang, S. J. Chem. Phys. 2006, 125, 174714. (22) Wang, Y.-T.; Cheng, C.-L.; Shih, Y.-C.; Kan, H.-C.; Chen, C.-H.; Hu, J.-J.; Su, Z.-Y. Can. J. Chem. 2007, 25, 1090–1093. (23) Rai, B.; Sathish, P.; Malhotra, C. P.; Pradip, C. P.; Ayappa, K. G. Langmuir 2004, 20, 3138. (24) Mazarello, R.; Cossaro, A.; Verdini, A.; Rousseau, R.; Casalis, L.; Danisman, M. F.; Floreano, L.; Scandolo, S.; Morgante, A.; Scoles, G. Phys. ReV. Lett. 2007, 98, 016102. (25) Pei, Y.; Ma, J. J. Am. Chem. Soc. 2005, 127, 6802. (26) Lahann, J.; Mitragotri, S.; Tran, T.-N.; Kaido, H.; Sundaram, J.; Choi, I. S.; Hoffer, S.; Somorjai, G. A.; Langer, R. Science 2003, 299, 371.
Molecular Dynamics Description of Grafted Monolayers (27) Duffy, D. M.; Harding, H. J. Langmuir 2005, 21, 3850. (28) Pertsin, J. A.; Grunze, M. Langmuir 2000, 16, 8829. (29) Zheng, J.; Li, L.; Tsao, H.-K.; Sheng, Y.-J.; Chen, S.; Jiang, S. Biophys. J. 2005, 89, 158. (30) Kong, D.-S.; Yuan, S.-L.; Sun, Y.-X.; Yu, Z.-Y. Surf. Sci. 2004, 573, 272. (31) Mikulski, P. T.; Herman, A. L.; Harrison, A. J. Langmuir 2005, 12197. (32) Srinivas, G.; Nielsen, S. O.; Moore, P. B.; Klein, M. L. J. Am. Chem. Soc. 2006, 128, 848. (33) Winter, N.; Vieceli, J.; Benjamin, I. J. Phys. Chem. B 2008, 112, 227. (34) Alexiadis, O.; Harmandaris, V. A.; Mavrantzas, V. G.; Site, L. D. J. Phys. Chem. C 2007, 111, 6380. (35) Goujon, F.; Bonal, C.; Limoges, B.; Malfreyt, P. Mol. Phys. 2008, 106, 1397. (36) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. M.; Fox, T.; Caldwell, J. W.; Kollman, P. J. Am. Chem. Soc. 1995, 117, 5179. (37) Ryckaert, J. P.; Cicotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (38) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (39) Sabalovic, J.; Liedl, K. R. Inorg. Chem. 1999, 38, 2764. (40) Smith, E. R. Proc. R. Soc. London, Ser. A 1981, 375, 475. (41) Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155. (42) Honenberg, P.; Kohn, W. Phys. ReV. A 1964, 136, 864. (43) Kohn, W.; Sham, L. J. Phys. ReV. A 1965, 140, 1133. (44) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (45) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (46) Perdew, J. P.; Wang, Y. Phys. ReV. E 1991, 45, 13244. (47) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.;
J. Phys. Chem. B, Vol. 112, No. 45, 2008 14229 Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian; Gaussian, Inc.: Wallingford, CT, 2004. (48) Breneman, C. P.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (49) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (50) DL_POLY is a parallel molecular dynamics simulation package developed at the Daresbury Laboratory Project for Computer Simulation under the auspices of the EPSRC for the Collaborative Computational Project for Computer Simulation of Condensed Phases (CCP5) and the Advanced Research Computing Group (ARCG) at the Daresbury Laboratory. (51) Marti, J.; Nagy, G.; Gordillo, M. C.; Guardia, E. J. Chem. Phys. 2006, 124, 094703. (52) Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press: London, 1985. (53) Spohr, E. J. Phys. Chem. 1989, 93, 6171. (54) Liu, X. Y.; Boek, E. S.; Briels, W. J.; Bennema, P. Nature 1995, 374, 342. (55) Rudnick, J.; Gaspari, G. Science 1987, 237, 384. (56) Milner, S. T. Science 1991, 251, 905. (57) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (58) Marry, V.; Turq, P.; Cartailler, T.; Levesque, D. J. Chem. Phys. 2002, 117, 3454.
JP8028825
