J. Phys. Chem. C 2007, 111, 2631-2642
2631
Interaction of Hydrated Amino Acids with Metal Surfaces: A Multiscale Modeling Description Pim Schravendijk, Luca M. Ghiringhelli, Luigi Delle Site, and Nico F.A. van der Vegt* Max-Planck-Institute for Polymer Research, Ackermannweg 10, D 55128 Mainz, Germany ReceiVed: August 28, 2006; In Final Form: NoVember 7, 2006
We present a multiscale modeling procedure that offers the opportunity to study hydrated biomolecules at metal surfaces. First principle DFT calculations and classical atomistic simulations are used interactively in order to account for both quantum and statistical aspects of molecular conformations at the surface. We present models for water, benzene, phenol, alanine, and phenylalanine at the (111) surface of nickel. These models are subsequently used in classical atomistic simulations to study physical-chemical aspects of amino acids at a Ni(111)/water interface. Application of this method to a larger set of “molecular building blocks” opens a computational route for molecular engineering of bio/inorganic hybrid systems.
I. Introduction A chemical realistic modeling procedure for describing peptide interactions with metal surfaces is of considerable scientific interest because of a complex interplay of quantum (electronic) and classical (solute and solvent) degrees of freedom that, so far, is only little understood. The ability to model these interactions in chemically realistic environments can support rational procedures for the on-demand design of peptides that specifically bind to inorganic surfaces (aptamers). This will help the design of bio/inorganic composite materials, for example, combining optimal solute binding specificity in proteins (e.g., antibodies) with the optimal signaling properties of inorganic materials.1 Surface binding up- and down-modulating amino acids in synthetic polypeptides have been identified recently by systematic experimental studies.2,3 These studies have only recently become possible because of technical developments such as combinatorial peptide engineering4,5 and have opened the way to predict, on an empirical basis, overall protein adhesion by considering adhesive properties of the constituting amino acids. In a previous paper, we outlined a procedure6 suitable to describe interactions of solvated molecules (i.e., benzene) with transition- and noble-metal surfaces. In the current paper, we extend upon that work by reporting the interaction of solvated amino acids with a Ni(111) surface. For this purpose, we use an iterative procedure that converges once the surface interactions described by the resulting classical model are consistent with those obtained from first principle quantum calculations, in terms of both conformational and energetic aspects. We apply this procedure to benzene, phenol, alanine, and phenylalanine. We, moreover, describe the modeling of water at the metal surface and combine that with the modeling of organic solutes to study amino acid surface interactions under aqueous conditions at 300 K. In view of the chemical complexity of biological molecules, as well as the system sizes that can be treated by first principle electronic structure calculations, it would be desirable to develop models describing surface interactions of isolated amino acid residues on a case-to-case basis. This means that we aim to * Corresponding author. E-mail:
[email protected].
develop molecular building blocks whose surface interactions are carefully parameterized and that can be combined to form any peptide sequence or larger biomolecule. This approach is similar in spirit to the development of biomolecular force fields in which reproducing condensed-phase properties of small molecules (e.g., resembling the amino acid residues) plays a key role. In contrast to biomolecular force fields, however, no general atomistic parameterization can be performed for the metal surface atoms because the delocalized electrons in the metal surface will be perturbed differently by different solutes, which can only be described with quantum mechanical calculations. The resulting parameter sets describing the interaction of specific solutes with the metal surface are therefore transferable only for chemical groups with the same chemical composition. They can be used to describe any peptide sequence, and interactions with the solvent can be incorporated using existing force fields. Obviously, this choice necessarily introduces approximations, which should be examined carefully whenever possible. We assume that the water-surface as well as all building block-surface interactions add in a pairwise additive fashion. This assumption should be reasonable as long as the chemically active groups on the biomolecule have relatively isolated electronic properties (e.g., a polar group in an alkane chain). Then, one may combine any number of those groups in a macromolecule, as long as the surface-adsorbed groups will be separated far enough to not disturb the electron distribution of the neighboring solute-metal interactions. Although the benzene and phenol models discussed in this paper can be used as building blocks in peptide sequences, the neutral amino acid models discussed in this work do not describe all features present in peptide bonds. In forthcoming work, we will therefore report the surface interaction of the peptide bond (based on a study of n-methyl acetamide), which may then be combined with further models for the amino acid side chain analogs to construct oligopeptides. This paper is organized as follows. After a description of computational details in Section II, we introduce in Section III the modeling procedure and the newly developed models and parameter sets for phenol, alanine, phenylalanine, and tyrosine interacting with the (111) surface of nickel. The modeling of water-surface interactions, described previously
10.1021/jp065568u CCC: $37.00 © 2007 American Chemical Society Published on Web 01/23/2007
2632 J. Phys. Chem. C, Vol. 111, No. 6, 2007 in ref 6, will be discussed more extensively in Section III as well. In Section IV, the newly developed models are used to study surface interactions under aqueous conditions at 300 K. II. Computational Details We extend the “molecular building block” approach that has already been applied successfully to organic-inorganic interfaces.6,7 In this approach, macromolecule-surface interactions are described at two levels. First, building-block molecules are chosen that describe the recurring parts in the macromolecule. The small sizes of these building blocks allow for quantum mechanical calculations of the building block-surface interactions, and the resulting data are used for parameterizing atomistic solute-surface interaction potentials. On the other, classical atomistic, level, this quantum-based parameterization describes building block-surface interactions, not only for the single building blocks but also for the structures such as (oligo)peptides or proteins that can be constructed from combinations of these building blocks. In the current paper, we will model phenol and the neutral forms of alanine, phenylalanine, and tyrosine and study their interaction with a Ni(111) surface from the bulkhydrated state. The Ni(111) metal surface has been chosen because its properties for similar systems have been well-tested in previous works.6-11 Extension to other surfaces (e.g., Pt, Pd, Au) is straightforward as will be discussed later on. As an ideal model surface, any oxidation effects at the surface will be ignored. Quantum calculations of the amino acids alanine and phenylalanine have been performed using the DFT based finiteelectronic temperature method of Alavi et al.12 (FEMD), as implemented in the CPMD code.13 Note that amino acids under physiological conditions exist mainly (>99%) in the zwitterionic state, whereas in the case of peptides the end groups are condensed to form peptide bonds and are therefore of different chemical nature. The neutral amino acids used here are therefore not the exact representations of any of those two states but are still of interest because the neutral state (as well as the zwitterionic state) has been shown to bind to a Ni(111) surface.14 The zwitterionic state of amino acids as well as the peptide form are currently under investigation14 and will be a matter of following publications. We used the Gromacs 3.3 code,15 adapted to allow for any combination of atom-wall potentials divided over the molecule, thus enabling a convenient parameterization of both configurational and energetic data from quantum calculations. Surfaces are introduced by atom-wall potentials acting in the z direction at both z borders of the box, and an empty space larger than any cutoff length is added to remove periodicity in the z direction. With surfaces at two sides of the box, the system actually represents a slit; however, the distance between the surfaces was chosen large enough (>7 nm) to enable a layer of bulk-like water in the center, several nanometers thick. The box size in the x and y direction extended 3.4 nm. Pressure coupling and particle mesh Ewald (PME) treatment of electrostatic interactions is performed by (semi) two-dimensional schemes that are part of the Gromacs 3.3 package. Pressure coupling was chosen to be coupled only to the box size in the x and y directions, parallel to the surface, and the box size in the z direction was set constant. The PME Fourier spacing was set at 0.12 nm, the real space and neighbor search cutoff was set at 0.9 nm, and van der Waals interactions were cut off at 1.4 nm. Simulations were carried out in a system similar to the one described in the previous paper,6 consisting of 3000 water molecules at constant pressure and temperature (NPT) (1 atm
Schravendijk et al. with a coupling constant of 1.0 ps) (300 K with a coupling constant of 0.5 ps), using Nose-Hoover temperature coupling16 and Parrinello-Rahman pressure coupling,17 but additionally introducing the molecules described above and using the LINCS bond constraint algorithm.18 All simulations were performed with a molecular dynamics (MD) integration time step of 2 fs. The solute-surface and solvent-surface potentials were treated independent from the solute-solute, solute-solvent, and solventsolvent intermolecular potentials, which were described by the GROMOS 43a1 force field.19 As a solvent, the SPC/E water model was used.20 The natural question arising at this point is if the use of a specific force field for water causes model dependencies. We address this question below and show that this is not the case. In addition to MD simulations, several Langevin dynamics (LD) simulations (300 K) were performed in which the solute interacts with the metal surface in a vacuum environment. The LD algorithm, available in Gromacs,15 was used with a friction coefficient of 1 ps-1. Solute-surface potentials of mean force (PMFs) in water were obtained via constraint-biased simulation with force averaging21 as provided by the Gromacs package. The constraint direction was chosen perpendicular to the surface (i.e., in the z direction). For each PMF, 142 defined constraint distances from the surface were chosen, in steps of 0.01 nm from 0.19 to 1.00 nm, in 0.02 nm steps from 1.00 to 2.00 nm and in 0.10 nm steps for distances up to 3.00 nm. In all cases, this proved to be far enough to reach a bulk-hydrated state of the solute molecule as indicated by a constant value for the PMF at these distances. Several constraint sites on the solute molecule were tested, each separately, to see if any site dependence for the PMF was present. The following interaction sites were chosen as constraint sites: for benzene and phenol the geometric center of the ring; for alanine the carbonylic oxygen and the amine nitrogen; for phenylalanine the ring center, carbonylic oxygen, and amine nitrogen. It is important that only the z component of the constraint site was kept fixed. The constraint site could move in the xy directions, and the remaining parts of the molecule could move in any direction as long as this movement would not displace the z distance of the constraint site. The starting conformations at the 142 distances mentioned above were generated by first solvating the molecule in the middle of the slit, following by pulling it to an adsorbed state (at z ) 0.19 nm). From there, the z distance was increased in a sequence of short pull runs. Once all starting points were generated, 3 ns production runs were performed for each z distance, enabling the calculation of distance-dependent mean forces and mean surface interaction energies. III. Modeling Surface Interactions The quantum-atomistic multiscale modeling of amino acidmetal surface interactions is more complex than our previous multiscale modeling of benzene/polycarbonate adsorption.6,7 This is largely due to the low symmetry and relatively large number of interaction sites in amino acids and the fact that these interaction sites are located too close to each other to be treated separately. The low symmetry greatly affects the number of quantum calculations necessary: even though alanine has fewer atoms than benzene, there are many more nonequivalent geometries in which it can be oriented relative to the surface. A complete quantum analysis of all of these configurations at all possible metal lattice sites would require a very expensive computational effort. An extension to our multiscale modeling procedure6 is introduced here to cope with the before-mentioned problems in a way that actually enhances the procedure’s general validity and applicability.
Hydrated Amino Acids at Metal Surfaces We start by emphasizing that the attention lies on finding global minimum energy conformations because these will be the dominant contribution to the sampling statistics during long runs at the atomistic level. Calculations on intermediate, higherenergy states would provide only superfluous information at the current modeling precision. Therefore, the starting point of the modeling will be a selection of initial configurations that aims to give us the extreme binding conformations. The electronic properties of the metal binding of these configurations were calculated by a series of quantum density functional calculations (geometry optimization),14 considering per configuration the various possible positions on the metal lattice. The data retrieved from these calculations give us the parameters (interaction sites, interaction strength, optimal distances), based on which an atomistic model is constructed. However, because we are studying a molecule with a large number of orientational degrees of freedom, the molecule is likely to get trapped in local energy minima during quantum calculations. We resolve this by performing LD simulations of the molecule in vacuum, applying the modeled molecule-surface interaction parameters, and sample its configuration space at a range of distances perpendicular to the surface, using the constraint method and sites described in the Methodology section. This set of LD runs can sample a larger amount of the phase space than could be done by quantum optimization runs, in a fraction of the time. If the LD simulation generates (only) the surface adsorbing conformations that were already found in our initial quantum calculations, then our classical model is considered complete. If, however, the LD simulation generates adsorbing (low-energy) conformations that were not previously found in quantum calculations, then these conformations are used as starting configurations in a new series of quantum-based structure optimizations. If this second set of quantum calculations results in interaction properties that differ from the LD simulations, then the atomistic modeling is adjusted. Then, the procedure is repeated (LD run with new modeling, quantum calculation of LD output), until consistency between the quantum and the atomistic level is reached and all LD lowest energy conformations correspond to stable conformations from quantum calculations with respect to the energy and the configuration. The newly modeled molecules are introduced below, starting with phenol which forms a logical extension to the benzene modeling we performed previously.6 First, however, we describe the modeling of water at the metal surface. a. Water. A large amount of research of either experimental or theoretical nature exists in the field of water-surface interaction.22-27,29 Recent development in this field focuses on the aqueous-solid interface.23,25,26 The basics of the atomistic modeling of the interaction of water with metal surfaces applied by us has been introduced in a previous publication.6 However, in that paper we considered only one classical water model. Here we briefly repeat how the water-surface interaction is modeled and then proceed further by showing that the relevant interfacial hydration properties are independent of the classical water model chosen to describe water-water interactions. To model the water-surface interaction, we make use of quantum density functional calculations of small clusters of water molecules at the metal surface. The water cluster configurations are chosen by considering all relevant ways for liquid water to interact with the surface. Examples are shown schematically in Figure 1. We base our modeling idea on a well accepted and by now proven statistical property of liquid water; that is, it is locally and instantaneously tetrahedral. This allows us to imagine that
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2633 locally at the surface water conformations must consist of full or half tetrahedra (see Figure 1) because of the confinement of the surface. Among possible arrangements like the ones shown in Figure 1, some are not allowed according to data available from quantum calculations, that is, those conformation displaying hydrogen down-like structures. Thus, we have a first screening of possible conformations that must be reproduced by our modeling. The second important point is that, because we want to explore many relevant configurations, we should carry out many calculations. Because first principles quantum calculations are computationally very demanding, a solution to this problem comes from the observation that the tetrahedral structure of Figure 1 can be described well by different combinations of substructures consisting of monomers, dimers, and trimers in different conformations. As shown in Figure 1, we then make what we would call “the first layer approximation’’; that is, only the molecule close to the surface and those directly bonded to that participate to the adsorption strength. The adsorption energy per water molecule in these allowed configurations is calculated using a quantum-based approach10 (see the caption of Figure 1 for the calculation used), and the results are used to parameterize a water-surface potential for classical simulations. We have chosen an attractive 10-4 potential to describe the interaction between water oxygen and the surface (see eq 1). No interaction was applied between water hydrogen and the surface. See Table 1 for exact values of the parameters used.
{
2π10-4
UAttr.10-4 )
[52(σz ) - (σz ) ] z e z 10
4
cutoff
z > zcutoff
0
(1)
In the classical simulation of the metal-water interface, waterwater interactions are accounted for by a classical force field, and therefore the energy of water-water hydrogen bonding should not be included in the water-surface interaction energy. b. Phenol. Phenol on Ni(111) has a maximum surface interaction energy that is nearly 15% lower than the maximum surface interaction energy of benzene (see Table 2 and ref 14). Its conformation at the surface has the hydroxyl oxygen lifted by 0.5 Å, which has some effect on the adjacent carbon atom as well. The hydroxyl hydrogen is pointing down, ending up at the same height as the phenylic carbon atoms (see Figure 2). As a basis of the atomistic modeling of phenol, the benzene model described previously6 is chosen. To describe the interaction with the surface, Morse potentials are used on the aromatic ring carbon atoms (Cr) (eq 2), and repulsive 10-4 potentials (eq 3) are used on the ring hydrogens (Hr). The model parameters are given in Table 1, and atom types are defined in Figure 3. For all attractive atom types, a zcutoff of 1.4 nm was used.
UAttr.Morse )
URep.10-4 )
{
{
M(1 - e-a(z-σ))2 - M z e zcutoff z > zcutoff 0
2π10-4 0
[52(σz ) - (σz ) + 53] z e σ 10
(2)
4
z>σ
(3)
Besides the substitution of one of the hydrogens of benzene by a hydroxyl group, only a reduction of interaction energy was needed, which was done by scaling down all carbon-wall interactions. The carbon next to the hydroxyl group (Cp) was scaled down to half the interaction energy of the other carbon atoms (Cr). A weak 10-4 repulsion was put on the hydroxyl
2634 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Schravendijk et al.
TABLE 1: Overview of Molecule-Ni(111) Surface Force Field Parametersa interaction
potential type
(kJ/mol)
σ (nm)
Water
Ow-Ni Hw-Ni
attractive 10-4 no interaction
Cr-Ni Hr-Ni
6.40
0.24
Benzene Morse (a ) 35 nm-1) repulsive 10-4
17.5 4.27
0.20 0.20
Cr-Ni Hr-Ni Cp-Ni Op-Ni Hp-Ni
Phenol Morse (a ) 35 nm-1) repulsive 10-4 Morse (a ) 35 nm-1) repulsive 10-4 attractive 10-4
15.8 4.27 7.96 1.00 0.70
0.20 0.20 0.20 0.25 0.22
Oc-Ni N-Ni CRβ-Ni Cv-Ni R-Ni
Alanine attractive 10-4 attractive 10-4 rep. Morse (a ) 6.0 nm-1) repulsive 10-4 repulsive 10-4
8.90 15.0 4.0 10.0 4.27
0.23 0.22 0.58 0.38 0.20
Oc-Ni N-Ni Hv-Ni Cv-Ni Cr-Ni Hr-Ni R-Ni
Phenylalanine attractive 10-4 attractive 10-4 repulsive 10-4 repulsive 10-4 Morse (a ) 35 nm-1) repulsive 10-4 repulsive 10-4
8.90 15.0 4.27 10.0 17.5 4.27 4.27
0.23 0.22 0.17 0.38 0.20 0.20 0.20
Oc-Ni N-Ni Hv-Ni Cv-Ni Cr-Ni Hr-Ni Cp-Ni Op-Ni Hp-Ni R-Ni
Tyrosine attractive 10-4 attractive 10-4 repulsive 10-4 repulsive 10-4 Morse (a ) 35 nm-1) repulsive 10-4 Morse (a ) 35 nm-1) repulsive 10-4 attractive 10-4 repulsive 10-4
8.90 15.0 4.27 10.0 15.8 4.27 7.96 1.00 0.70 4.27
0.23 0.22 0.17 0.38 0.20 0.20 0.20 0.25 0.22 0.20
a
See Figure 3 for the corresponding structures and eqs 1-4 for the potential functions used. The “R” covers all interactions with no specific interaction with the surface (unlabeled atoms in Figure 3).
TABLE 2: Ab Initio Interaction Energies and Optimal Distances Used to Model Interaction Potentials for Classical MDa interaction 10
water-Ni benzene-Ni 44 phenol-Ni 44 Ala-Ni Phe-Ni 9
Eads (eV)
dopt (nm)
-0.25 -1.05 -0.9 -0.57 -1.1
0.24 0.20 0.20 (ring) 0.22 (N) 0.20 (ring)
a For neutral phenylalanine (Phe-Ni), only the optimal configuration was evaluated. For neutral alanine (Ala-Ni), the energy of the configuration with highest surface interaction is given. A complete overview is given in Table 3.
oxygen (Op), and an attractive 10-4 potential (see eq 1) was put on the hydroxyl hydrogen (Hp). Configurations of phenol with the OH group pointing away from the surface did not need to be taken into account in QM calculations, because the perturbing effect of the OH group will be lower the further it is away from the surface, and any interaction would be similar to the inclination dependence of benzene modeled previously.6 The configuration where phenol is pointing with the OH group toward the Ni(111) surface is currently being investigated by quantum calculations.30 Preliminary calculations indicated a weak interaction, not signifi-
TABLE 3: Properties of Four Neutral Alanine Conformations That Were Evaluated in Quantum Mechanical DFT Calculations14 a conformation
Eads (eV)
N (nm)
Oc (nm)
Cβ (nm)
NdownCHup NdownCHdown OdownCHup OdownCHdown
-0.57 -0.32 -0.37 nonbonding
0.22 0.24 0.46 0.39
0.32 0.33 0.25 0.22
0.46 0.36 0.57 0.27
a
Shown are the adsorption energy on the Ni(111) surface, and the atom-to-top Ni layer distance of three of the atoms whose surface distance most strongly affects the surface interaction.
cantly higher than the error margin of the quantum calculation. Currently, no special modeling considerations were made to account for this orientation in the classical atomistic simulation. c. Neutral Alanine. Four distinct conformations were used in the initial quantum calculations of neutral alanine, as described in ref 14. As an example, the strongest binding conformation (0.57 eV) is shown in Figure 4; the energies and optimal surface interaction site distances of all configurations are shown in Table 3. The amine nitrogen (N) and the carbonyl oxygen atom (Oc) interact attractively with the surface. Surface attraction at these sites is of the order of one to two hydrogen bond strengths and is modeled with attractive 10-4 potentials. An interesting feature is the influence of the methyl side-group. Even when not directly interacting with the surface, the methylsurface distance is found to influence the total surface interaction energy. The origin of this influence might be an effect of the position of the methyl group on the surface orientation of the interacting amine and carboxyl groups. A simple scheme to model this effect is chosen that includes a weak repulsive Morse potential (eq 4) on the center of mass of the CR and Cβ atoms: this additional interaction site is referred to as CRβ (see Figure 3).
URep.Morse )
{
M(1 - e-a(z-σ))2 z e σ 0 z>σ
(4)
The optimal distances (σ values) for these potentials are retrieved from the binding conformations in our quantum calculations. Binding energies for the potentials ( values) are chosen such that the sum of all site-surface potentials (eqs 1-4) acting on the molecule will reproduce the total surface interaction energy as found in the quantum calculation, for every conformation tested. Atoms with no specific surface interaction were given a simple repulsive 10-4 potential (with parameters based on the repulsive hydrogens in the benzene ring, see group R in Table 1) to prevent the occurrence of unphysical conformations corresponding to a penetration of the surface. The first iteration of this modeling resulted in minimal energy LD conformations that had both the amino nitrogen and the carbonyl oxygen interaction sites at optimum distance. The total surface interaction energy corresponded to the sum of both interactions. Consecutive quantum calculations of the LD output structures showed that this was actually a nonbonding conformation. Consequently, a modification of our initial modeling was needed to ensure that only one interaction site at a time will be able to bind with optimal interaction energy. This was reached by introducing a repulsive site (Cv) at a position between the amino nitrogen and the carbonylic oxygen (at 40% of the bond length between carbonylic carbon and CR, see Figure 3), leading to a seesaw-like mechanism that allows for the binding of either the amino or the carbonyl group. LD runs with this modeling lead to optimal configurations close to the previous quantum calculations, see Figure 5, and this modeling was therefore
Hydrated Amino Acids at Metal Surfaces
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2635
Figure 1. Graphical aid showing the rationale behind the ab initio modeling. The water configurations studied by quantum mechanic DFT calculations were chosen by considering all relevant ways for liquid water to interact with the surface. Some examples are schematically shown here. In white circles are tetrahedral water substructures at the water-metal surface. I: A possible tetrahedral water structure with two waters interacting with the metal surface. II: A trimer, the upper part of a tetrahedral water structure, with one water interacting with the metal surface and two additional hydrogen-bonded water molecules. III: One of many configurations that could be discarded immediately because no electronic hydrogen-metal interaction exists. The isolated water tetrahedron is an unstable structure, both in vacuo and metal-bound. Therefore, ab initio calculations are done with smaller subunits. In DFT calculations, we can represent the tetrahedral structures shown in I and II by using all relevant water structures consisting of (A) a water monomer, representing one of the symmetry axes of the tetrahedron; (B) a water dimer, representing the metal-bound water with its first hydrogen-bonded neighbor (the tetrahedral center); (C) a water trimer, directly representing structure II. The bare water-metal interaction energies (excluding contributions of H2O-H2O hydrogen bonding) are calculated as Eads ) E(surf+Nmol) - Esurf - ENmol, where Nmol refers to the number of water molecules in the system, E(surf+Nmol) denotes the QM energy of the combined surface + water substructure, Esurf is the QM energy of the isolated metal surface, and ENmol is the energy of the isolated water substructure. Eads/NH2O equals 0.25 ( 0.05 eV on Ni(111) independent of NH2O ()1, 2, or 3). The complete overview of all conformations evaluated in quantum density functional calculations is given in ref 10.
Figure 2. Minimal energy phenol configuration when adsorbed on Ni(111), as found by quantum calculations.44 The green area represents the location of the nickel surface.
chosen for our MD production runs. The final potential parameters are given in Table 1. Note that even though only one point (minimized energy) per conformation is used, the fact that we use four different conformations for our potential fit, with various groups interacting simultaneously for every conformation, makes this more elaborate than a one-point parameterization. d. Neutral Phenylalanine and Tyrosine. Phenylalanine and tyrosine can, in a building block manner, be constructed by combining characteristics of the alanine and the benzene/phenol modeling. The repulsive CRβ site, present in the alanine model, is discarded in the phenylalanine and tyrosine models. Instead, because repulsive aliphatic hydrogen atoms proved to be a necessary modeling element to prevent unphysical conformations with Cβ hydrogens penetrating the surface, these were introduced (one on the CR, two on the Cβ) as virtual sites in the phenylalanine and tyrosine molecules. We use the terminology “virtual” sites because, apart from experiencing surface repulsion, they do not interact with the solvent; in the GROMOS 43a1 force field19 hydrogens connected to aliphatic CR and Cβ carbons are absorbed into united-atom potentials centered on the carbon positions. Interactions of other atoms (e.g., solvent) with these aliphatic CHn groups are described using these united atom potentials. Virtual sites in Gromacs are built from the coordinates of 2, 3, or 4 vicinal atoms; any forces on the virtual sites are spread out over these atoms at the end of each time step. The bond lengths and angular potentials of our virtual sites were taken from the all-atom OPLS force field.31 In the resulting lowest energy conformation obtained from LD simulations in vacuum, we found that both the amine nitrogen and the phenyl
ring could reach their positions of maximum interaction energy, leading to a total surface interaction energy of 1.59 eV (153 kJ/mol). This is different from the initial quantum calculation corresponding to a 1.1 eV binding conformation, where the amino and carboxyl group were chosen to point away from the surface, to minimize surface contact of the aliphatic Cβ hydrogens.9 Following our iterative multiscale modeling procedure, we tested the lowest energy LD conformation in a consecutive quantum calculation,32 and it was found to be a stable conformation (1.5 eV). One aliphatic Cβ hydrogen is close to the surface but can apparently find a stable position within a hollow site of the surface. Because the optimal LD interaction energy and conformation were close enough to the optimal conformation found by quantum calculations (within the 0.1 eV error of the quantum calculation), no further optimization of the modeling was needed. The final combination of atomsurface potentials is given in Table 1 and Figure 3. IV. Analysis of Surface Interactions at the Ni(111)/H2O Interface For benzene and phenol, an in-depth analysis can be made of the surface interaction mechanism under aqueous conditions. Because of the high symmetry of the benzene ring, its centerof-mass-surface distance was recently shown to provide a logical choice for an order parameter along which a free energy can be calculated,6 by which the surface interaction can be characterized. When extending this analysis to a phenol molecule, we meet with a slightly broken symmetry; we can, however, still take the distance between the geometric center of the ring and the surface as the order parameter of interest. For benzene and phenol, their planar structures enable us to study not only the PMF’s z-dependency but also the surface inclination dependence described by the angle between the normal of the aromatic ring and the normal of the surface, at a given z distance, as will be discussed in more detail below. Following the surface interaction mechanisms of hydrated amino acids is more problematic, not only because of their low symmetry as compared to a benzene ring but also because of their many degrees of conformational freedom and various interaction sites. Therefore, several order parameters can be
2636 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Schravendijk et al.
Figure 4. One of the four alanine-Ni(111) conformations as mentioned in Table 3. Main characteristics of this conformation (NdownCHup) are the binding of the amine nitrogen to a top site of the Ni(111) surface and the methyl group pointing away from the surface. The carbonylic oxygen is pointing slightly toward the surface. A detailed quantum analysis of this and other structures is given in ref 14.
potential of mean force based on simulations of one solute molecule in water.
Figure 3. Atom types used in Table 1. Shown are (A) benzene, (B) phenol, (C) water, (D) neutral alanine, (E) neutral phenylalanine, and (F) neutral tyrosine.Atom names refer to the atom names in Table 1. The sites CRβ and CV are specially introduced for the current modeling: CRβ is a virtual site located exactly at the center of mass of CR and Cβ. CV is a repulsive virtual site needed to prevent simultaneous adsorption of the N and carbonyl O (OC) atoms within the same molecule, as explained in the text. Note that the water hydrogens (HW) do not have any interaction with the surface. Unlabeled atoms have the general repulsion “R”.
chosen to follow the process, all of which will, however, be interdependent: at any surface distance of a given interaction site, all other interaction sites will contribute to the PMF at that particular distance. Here, we limit ourselves to presenting only the free-energy difference between the bulk-hydrated state and a selection of local free-energy minimum states with strong surface interaction, and a description of the interaction energies involved. To better understand hydration effects, we will also compare the explicit solvent MD simulations with the LD simulations of the amino acids in vacuum that were performed during our modeling procedure. The total combination of relative surface interaction energies for the various amino acids, and comparisons of surface interactions under solvated and nonsolvated conditions, can give essential information to help us understand the main chemical factors that determine metalsurface interaction of amino acids, and in this way facilitate the design of surface-binding peptides. a. Water. We used our water modeling procedure to describe the Ni(111)/water interface using different classical water models. Applied were the three-site SPC,33 SPC/E,20 and TIP3P34 models, the four-site TIP4P34 model, and the five-site TIP5P35 model. For each water model, we studied the water structure at the surface and the benzene-surface and phenol-surface
Figure 6 compares the water structure near Au(111) and Ni(111) for the SPC, SPC/E, TIP3P, TIP4P, and TIP5P water models. The water-metal interaction for Ni(111) amounts to 0.25 eV per molecule (24.1 kJ/mol) at an optimal distance of 2.4 Å; for Au(111) it amounts to 0.10 eV per molecule (9.7 kJ/mol) with an optimal distance of 3.1 Å.6 Clearly, the structure of the metal-water interface is independent of the classical water model. On Ni(111) the hydrogen density profile (dashed line) is highly structured (four peaks for z < 7Å) because hydrogen atoms belonging to water molecules in the first adsorbed water layer either correspond to OH bonds aligning the surface or OH bonds hydrogen bonded to water molecules in the second adsorbed layer. Hydrogens belonging to water molecules in the second water layer, in turn, donate hydrogen bonds to waters in the first layer and water molecules in the bulk. On Au(111) the hydrogen density profile is less structured, indicating reduced hydrogen bonding between the water layers. The peak height for the first water layer on Au(111) is in the order of 1/3rd of the peak height of the first water layer on Ni(111). The relative peak heights are, however, not a good measure of the relative water densities at the Au(111) and Ni(111) surfaces because these peaks are narrower than the molecular diameter of water and the density varies rapidly over this range. Therefore, we integrated the water (oxygen) density profile over the first peak up to the first minimum. The resulting number, expressed in units area per molecule, is given in Table 4 for all five water models. For Au(111), surface areas in the range of 9.3-9.9 Å2 per molecule were found, for Ni(111) the range was 8.4-8.9 Å2 per molecule. Note that these values are in the same order as those reported by Shelley et al.28,29 for atomistic simulations of a water slab near Hg surfaces. Different water models produce only small differences in the structure of the metal-water interface, independent of the
Hydrated Amino Acids at Metal Surfaces
Figure 5. Justification of the atomistic modeling. Solid line: average alanine-surface interaction energy (A and B) and site-surface distance (C and D) obtained from a series of 3 ns constrained LD simulations in vacuum at 300 K. Dots: optimal energy (A and B) and optimal distance (C and D) of zero-temperature quantum calculation optimized structures. In part C, the distance of the carbonylic O is shown (vertical axis) for a given constrained distance of the amino N (horizontal axis). In part D, the distance of the amino N (right vertical axis) is shown for a given constrained distance of the carbonylic O (horizontal axis). The LD data (solid lines) were obtained by constraining either the N or the O surface distance in the simulation. The distances and corresponding energies from the quantum calculations are all retrieved from the three binding configurations found in quantum calculations.14 The atomistic model samples the quantum-based conformation and energy either correctly or within a 0.05 nm distance. Parameters for this modeling are mentioned in Table 1.
strength of the water-metal interaction. This is an additional confirmation that the current modeling procedure allows for a separated treatment of molecule-surface interactions while using existing force fields for molecule-molecule interactions. b. Benzene and Phenol. Both the benzene and phenol z-dependent PMFs (Figure 7) are qualitatively very similar. The structural details of the PMFs follow the water oxygen density fluctuations (Figure 6). Upon approaching z ) 0.5 nm (corresponding to the position of the second oxygen peak in the water density profile) from larger distances, the benzene (phenol) molecule starts to expel water from the second water adsorption layer and the PMF increases rapidly. The PMF profiles correlate with the water density data from Figure 6, where we can find a first water layer around z ) 0.2-0.3 nm, a second water layer around z ) 0.5-0.6 nm, and a very weak third water layering around z ) 0.8-0.9 nm. Our previous paper6 showed that the potential of mean force (PMF) for displacing a benzene molecule from bulk water to the surface not only gives an idea about adsorption/desorption energies and intermediate energy barriers but also helps in understanding the critical steps in the adsorption/desorption mechanism.6 Any crucial effect that the choice of classical water model would have on the adsorption mechanism would therefore be visible in the corresponding PMF. Figure 7 shows that with all classical water models similar features are observed in the PMF, with a major free-energy barrier (50-60 kJ/mol high) starting almost directly at the surface (0.20 nm), ranging until ca. 0.60 nm. Additional oscillations occur further away from the surface up to ca. 1.50 nm. Differences are found in the barrier heights and the free-energy difference between the bulk-
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2637
Figure 6. Normalized water densities near Au(111) (top graph) and Ni(111) (bottom graph) for SPC,33 SPC/E,20 TIP3P,34 TIP4P,34 and TIP5P35 water models. Solid line, water oxygen; dashed line, water hydrogen.
TABLE 4: Surface Area Per Water Molecule in the First Adsorbed Layera area per molecule (Å2) water model
Au(111)
Ni(111)
SPC SPC/E TIP3P TIP4P TIP5P
9.52 9.27 9.46 9.61 9.86
8.64 8.37 8.61 8.79 8.92
a Determined by the integral of the oxygen density peak adjacent to the surface in Figure 6.
hydrated (z > 3 nm) and surface-adsorbed state (z ) 0.2 nm). These variations are partly accounted for by the error margins of these PMF calculations (between 4 and 8 kJ/mol at z ) 0.2 nm), see the caption of Figure 7, and small differences in water densities at the surface. In addition, differences in bulk water density and solvation free energies for the various water models36-40 are likely to be a source for variations in the resulting free-energy difference. In Figure 8A, two-dimensional PMFs for benzene and phenol are shown as a function of z and the inclination angle between the aromatic ring- and surface normals. A selection of snapshots of benzene at various surface distances is provided in Figure 8B. In this two-dimensional free-energy landscape, a more detailed picture is shown, and we can subdivide several zones when approaching the surface (z ) 0.2 nm) from the bulk
2638 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Schravendijk et al.
Figure 7. Solute-surface PMFs for displacing the geometrical center of the phenyl ring of (A) benzene and (B) phenol perpendicular to a Ni(111) surface, in liquid water (300 K) described with five classical water models. PMFs were obtained by integrating the average constraint force on the solute center of mass along 140 discrete points between 0.20 and 3.00 nm from the metal surface. At each point, 3 ns MD runs were performed to sample the mean constraint force. Solid black line, SPC water;33 dashed black line, SPC/E water;20 dashed/dotted black line, TIP3P water;34 solid gray line, TIP4P water;34 dashed gray line, TIP5P water.35 The errors of these PMF calculations ranged between 4 and 8 kJ/mol at z ) 0.2 nm and were estimated by calculating the block average error estimate45 for every constraint distance and taking the square root of the integral of the square of all error estimates from bulk (z ) 3.0 nm) toward the closest distance to the surface that was sampled (z ) 0.20 nm).
solvent (z > 1 nm). In the bulk zone going from z > 1 nm toward z ) 0.8 nm the sampling is distributed evenly over all angles. For benzene (less clear for phenol), it can be seen that around 0.7 nm (in between the third and second water layer; see Figure 6) already some structural effects occur and the parallel conformation (cos(θ) ) 1) is slightly less sampled. At distances corresponding to the second adsorbed water layer (between z ) 0.6 nm and z ) 0.5 nm), benzene (phenol) preferentially orients parallel to the metal surface. (The distribution of cos(θ) is narrowed down to 0.9-1.0 in Figure 8A. A snapshot of the parallel conformation can be seen in Figure 8B, at 0.49 nm.) This orientation is favored energetically because of O-H‚‚‚π hydrogen bonding involving water molecules in the first and second surface hydration layers. This type of weak water-aromatic hydrogen bonding is described properly by the force field.41 In Figure 7, this causes a shoulder to appear in the PMF at z ) 0.58 nm. A rather flat landscape with respect to the orientational degree of freedom is found in the region between the second and first adsorbed water layer (0.4 nm < z < 0.5 nm). Although, toward z ) 0.4 nm, perpendicular orientations are favored in comparison to parallel ones. In Figure 7 a shoulder is observed at these distances, and in Figure 8B one can see how a perpendicular orientation minimizes the displacement of water molecules in both the first and second adsorbed water layer. As the geometrical center of the ring approaches the first adsorbed water layer, the benzene (phenol) ring normal gets significantly tilted with respect to the surface normal and finally lines up with the surface normal, driven by an energetic stabilization of 1 eV due to the overlap of benzene(phenol)-π-orbitals with free electrons in the surface. In this process, first layer water molecules are expelled from the hydrophilic nickel surface, even before significant benzene(phenol)-surface binding interaction is present. It is clear that if the benzene would adsorb in a solvent-free environment then the z-dependent angular distribution plot would show a random orientation for distances of ca. 0.5 nm and higher and would
show a similar profile to the solvated state only for short distances below 0.3 nm, because here the orientation of the ring is governed by the presence of the surface. c. Amino Acids. To study surface-interacting amino acid conformations, we performed several MD runs in which a single interaction site-surface distance (order parameter) was constrained and a PMF as a function of that order parameter was obtained by integrating the mean force starting from z ) 3 nm downward toward the surface. Using this method, the distances with minimal free energy were determined. Snapshots of corresponding conformations as well as the average surface interaction energies and relative free energies (as compared to z ) 3 nm) at these distances are shown in Figure 9. Comparison of the free energies in Figure 9A and B as well as Figure 9C and D shows that, independent of whether the distance constraint is applied on the oxygen or nitrogen, similar values are obtained, which is an indication that sufficient sampling has been reached. The snapshots show that the amino acid conformations in Figure 9A and B as well as those in Figure 9C and D are similar. The free energies for phenylalanine viewed from the amino N (Figure 9C) and carbonylic O (Figure 9D) order parameters are lower than those for the alanine cases (about 9 kJ/mol). This is rather counterintuitive because the bulky phenyl ring did not contribute to the surface interaction energy for these order parameters, but intuitively its excluded volume can be assummed to displace more water molecules from the surface as compared to the methyl group from alanine. An explanation could be that in the case of alanine the small methyl group is able to come close to the surface, thereby weakening the interaction of amine and carbonyl groups. In contrast, the bulky phenyl group of phenylalanine will not be able to pass the hydration layers present at the surface and therefore will not hinder binding of the amine and carbonyl groups.
Hydrated Amino Acids at Metal Surfaces
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2639
Figure 8. (A) Two-dimensional PMFs for benzene (left) and phenol (right) at Ni(111) in SPC/E water (300 K). The solute (center-of-mass)-tosurface distance is plotted vertically; the cosine of the angle between the solute- and surface normal vectors is plotted horizontally. The 2D-PMF, G(z,θ) ) G(z) + G(θ|z), was calculated from constrained MD. By applying a constraint to keep the z coordinate fixed, G(z) (see Figure 7) is obtained by integrating the mean constraint force along the z coordinate. G(θ|z) ) -kBT lnP(θ|z) is obtained from the conditional distribution function, P(θ|z), sampled in the constrained MD runs. (B) Benzene snapshots at a selection of center-of-mass constraint distances from the surface.
The large positive free energy obtained by displacing the phenyl ring in phenylalanine toward the surface (Figure 9E) requires additional comments. This conformation has the lowest total solute-surface interaction energy of all structures that interact with the surface as found here, but it is the least favorable from a free-energy point of view (∆Gads ) +32.6 kJ/mol ( 19.8 kJ/mol). A reasoning to explain this observation would be to consider the fact that transferring the ring to the surface requires a displacement of more water molecules than in all other cases. In Figure 10 the surface interaction energies for water, phenylalanine, and the interaction sites of phenyla-
lanine are plotted against the distance of the constrained center of mass of the phenyl ring. The contribution of water-surface interactions decreases over the distance from bulk to optimal adsorption with 140 kJ/mol, about 70 kJ/mol more than for the other cases, where, instead of the phenyl ring, the O or N interacts with the surface (data not shown). One can get an impression of this solvent effect when comparing the phenylalanine surface interaction in explicit solvent MD (Figure 9E) with LD runs in vacuum (Figure 9F), where we find free-energy minima at the same distance (geometrical center of the ring at 0.22 nm), similar interaction energies (-134.7 kJ/mol in explicit
2640 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Schravendijk et al.
Figure 9. Snapshots of amino acids and parts of the surrounding water at the free-energy minimum distance of the constraint sites chosen, determined by the PMF method as explained in the text. Eint denotes the sum of interaction energies for all solute interaction sites with the surface. ∆G denotes the free energy of surface interaction, taken by the difference in the PMF between the distance given here, and the bulk state at z ) 3.0 nm. Error bars are calculated by calculating the block average error estimate45 for every constraint distance and taking the square root of the integral of the square of all error estimates from bulk (z ) 3.0 nm) toward the constraint distance mentioned here. ∆G for a given molecule is found to be independent of the interaction site (within the error of the calculation) and dramatically dependent on the presence of water. (A) Alanine, amino N constrained at z ) 0.24 nm. (B) Alanine, carbonyl O constrained at z ) 0.28 nm. (C) Phenylalanine, amino N constrained at z ) 0.24 nm. (D) Phenylalanine, carbonyl O constrained at z ) 0.29 nm. (E) Phenylalanine, geometrical center-of-ring constrained at 0.22 nm. (F) As in part E, but taken from an LD simulation in vacuum.
solvent MD and -137.8 kJ/mol in vacuum LD), but largely different free-energy values (+32.6 kJ/mol in explicit solvent MD vs -93.3 kJ/mol in vacuum LD). By looking at the waterand solute-surface interaction energies in Figure 10, it is apparent that a competition between water and nitrogen plays an important role: every increase in solute-surface interaction energy is accompanied by a decrease of water-surface interaction energy and vice versa. The similarities in the surface interaction mechanism of phenol and benzene (Figure 8) indicate that the tyrosine-surface interaction is similar to that of phenylalanine. V. Discussion Four factors are necessary at the level of computer simulations in order to get a complete theoretical picture of peptide adsorption on metal surfaces. First, any molecule-metal surface interaction has to be parameterized by ab initio methods and cannot be represented by generalized force fields. An efficient
procedure to obtain a correct parameterization is presented in the current paper. The second factor is related to the competitive adsorption energies of the solvent and solute, where one should take into account the number of solvent molecules that have to be displaced by the solute.6,42 It has become apparent, however, that taking these competitive effects alone does not provide the full description of the system because it does not properly account for solute hydration and the surface hydrophilicity.6 This is the third factor that should be taken into account: depending on surface hydrophilicity, dense adsorbed water layers may exist close to the surface. Intrusion of these mutually hydrogenbonded layers causes energy barriers for surface approach. Therefore, explicit solvent (atomistic) simulations are necessary, using timescales long enough to allow for solvent rearrangements. The fourth factor concerns the geometry and orientation. Especially for longer molecules like polypeptides, many conformations exist next to the surface and a correct sampling has
Hydrated Amino Acids at Metal Surfaces
Figure 10. Dependence of surface interaction energies (averaged over 3 ns simulations) of phenylalanine and water on the phenyl-ring-surface z-distance at 300 K. Shown are the overall solute-surface interaction energy (blue line), the water-surface interaction energy (green line, right vertical axis), phenyl ring-surface interaction energy (red line), and amino nitrogen-surface interaction energy (brown line). Water is expelled as the ring approaches the surface, causing a loss in watersurface interaction energy of approximately 140 kJ/mol over the distance from 0.75 to 0.20 nm. Over the same distance the solute gains a maximum of 130 kJ/mol at 0.22 nm (see snapshot E in Figure 9). The statistical noise in the water-surface and solute-surface interaction energy profiles is correlated due to a competition of water and amino nitrogen binding. The large statistical uncertainty (19.8 kJ/mol) in the free-energy calculation (snapshot E in figure Figure 9) is due to the noisy energy profile. The free-energy barrier in Figure 7 can be explained from the water-surface and solute-surface interaction energies: for 0.35 nm < z < 0.6 nm the water-surface interaction energy decreases and is not compensated by an increase of solutesurface interaction energy.
to be performed to find all possible surface-interacting conformations. Because of solvent effects, the outcome may well be non-trivial. The multiscale simulation approach presented here combines configurational and chemical information needed for engineering surface-interacting peptides. The simulations provide insight into the mechanisms of surface interactions in hydrated systems and can therefore directly be applied to support and explain observations in experimental studies. Obviously, the approach described here requires extension. QM calculations of the peptide group (CONH) and water interacting with Pt surfaces are currently being performed by us. Future work on additional modeling of amino acid residue interactions with this surface will result in a molecular construction set opening the way to the modeling of a variety of peptide-surface systems. Amino acid interactions with the nickel surface modeled in this work are transferable to surfaces of different chemical composition. The surface interaction energies of benzene and water can be arranged as a series with comparable energies for Ni, Pt, Pd, and Rh.6 Therefore, the mechanisms described here for nickel are likely to be similar on Pt, Pd, and Rh. Several preliminary generalizations comparing our modeling with experimental results can then be made. First of all, our simulations explain why phenylalanine, for which QM density functional calculations indicate that it binds strongly to a metal surface (i.e., Ni, Pt, Pd), is not found within the strong metal binders experimentally.2,3 The strongest binding configuration of a phenylalanine molecule in vacuo (the interacting aromatic ring oriented parallel to the surface) is shown by us to correspond with a highly unfavorable free energy when in the presence of solvent (Figure 9). Although in our simulations phenylalanine is a better binder than alanine (contrary to experimental findings on oligo-peptides2,3), one has to take into account that our simulations concerned only single amino acids. We will at a later point extend our simulations to peptide chains,
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2641 for which different behavior may be observed. Tyrosine is experimentally found to be a relatively strong binder among the uncharged amino acids.2,3 Because our simulations showed that the planar-ring conformation is unlikely to bind to a hydrated surface, the most reasonable explanation will be an interaction of the phenol hydroxyl group in tyrosine with the metal surface. Because of electronic polarization effects, surface defects might contribute to this interaction. Here it should also be kept in mind that in most experiments polycrystalline noble metals are used;2,3 hence, interactions with alternative crystal planes and surface defects require attention in future calculations. We note that QM density functional calculations of the adsorption of water10 and benzene43 on metal surface defects have already been performed recently. VI. Conclusions We have reported an iterative multiscale modeling procedure that uses (1) quantum density functional calculations to obtain surface interaction energies and optimum distances and (2) classical atomistic simulations to overcome energy barriers and guide localizing global conformational minima. After consistency between the atomistic modeling and the quantum calculations has been reached, the fast sampling obtainable with atomistic simulation is used to investigate interactions of biological molecules with metal surfaces in water. Following this approach, it is now possible to take into account electron delocalization effects in interactions of amino acids with metal surfaces and to fully consider the water-solute and watersurface interactions present in hydrated systems. These advancements are essential in approaching a realistic modeling of experimental peptide-surface systems. As an application, the multiscale modeling of hydrated phenol, alanine, phenylalanine, and tyrosine has been performed, and various adsorption properties (adsorption energy, free energy, structural information) have been obtained. Several general conclusions concerning the chemistry of surface adsorption can be drawn from the current study. Most importantly, it is found that quantum-based binding energies alone do not suffice to understand the thermodynamic aspects of protein-surface interactions. Instead, one should account for the competing effects of solvent and other adsorbing groups, potentially introducing free-energy barriers for surface approach. Additionally, our multiscale modeling can be applied to study geometrical effects in detail, as exemplified with our analysis of benzene and phenol, which can help in understanding the kinetics of adsorption processes. Acknowledgment. We thank Kurt Kremer, Matej Praprotnik, and Hatice Duran for helpful discussions. References and Notes (1) Sano, K.-I.; Sasaki, H.; Shiba, K. Langmuir 2005, 21, 3090. (2) Peelle, B. R.; Krauland, E. M.; Wittrup, K. D.; Belcher, A. M. Langmuir 2005, 21, 6929. (3) Willet, R. L.; Baldwin, K. W.; West, K. W.; Pfeiffer, L. N. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 7817-7822. (4) Brown, S. Nat. Biotechnol. 1997, 15, 269. (5) Whaley, S. R.; English, D. S.; Hu, E. L.; Barbara, P. F.; Belcher, A. M. Nature 2000, 405, 665. (6) Schravendijk, P.; van der Vegt, N.; Delle Site, L.; Kremer, K. ChemPhysChem 2005, 6, 1866. (7) Delle Site, L.; Abrams, C. F.; Alavi, A.; Kremer, K. Phys. ReV. Lett. 2002, 89, 156103. (8) Abrams, C. F.; Delle Site, L.; Kremer, K. Phys. ReV. E 2003, 67, 021807. (9) Delle Site, L.; Kremer, K. Int. J. Quantum Chem. 2005, 101, 733.
2642 J. Phys. Chem. C, Vol. 111, No. 6, 2007 (10) Sebastiani, D.; Delle Site, L. J. Chem. Theory Comput. 2005, 1, 78. (11) Murakhtina, T.; Delle Site, L.; Sebastiani, D. ChemPhysChem. 2006, 7, 1215. (12) Alavi, A.; Kohanoff, J.; Parrinello, M.; Frenkel, D. Phys. ReV. Lett. 1994, 73, 2599. (13) Hutter, J.; Alavi, A.; Deutsch, T.; Bernasconi, M.; Goedecker, S.; Tuckerman, M.; Parrinello, M. CPMD V. 3.4.1, 1995-1999. (14) Ghiringhelli, L.; Schravendijk, P.; Delle Site, L. Phys. ReV. B 2006, 74, 035437. (15) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701. (16) Nose´, S. Mol. Phys. 1984, 52, 255. (17) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182. (18) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463. (19) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide; vdf Hochschulverlag AG an der ETH Zu¨rich, 1996. (20) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (21) McDonald, I. R.; Rasaiah, J. C. Chem. Phys. Lett. 1975, 34, 382. (22) Thiel, P. A.; Madey, T. E. Surf. Sci. Rep. 1987, 7, 211. (23) Guidelli, R.; Schmickler, W. Electrochim. Acta 2000, 45, 2317. (24) Henderson, M. A. Surf. Sci. Rep. 2002, 46, 1-308. (25) Michot, L. J.; Vilie´ras, F.; Franc¸ ois, M.; Bihannic, I.; Pelletier, M.; Cases, J.-M. Geoscience 2002, 334, 611. (26) Taylor, C. D.; Neurock, M. Curr. Opin. Solid State Mater. Sci. 2005, 9, 49.
Schravendijk et al. (27) Verdaguer, A.; Sacha, G. M.; Bluhm, H.; Salmeron, M. Chem. ReV. 2006, 106, 1478. (28) Shelley, J. C.; Patey, G. N.; Be´rard, D. R.; Torrie, G. M. J. Chem. Phys. 1997, 107, 2122. (29) Shelley, J. C.; Be´rard, D. R. ReV. Comput. Chem. 1998, 12, 137. (30) Ghiringhelli, L.; Delle Site, L. In preparation. (31) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (32) The setup for the DFT calculation of this phenylalanine conformation is the same as that specified in ref 14, except that a 4 × 4 Ni(111) surface was always used. (33) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel Publishing Company: Dordrecht, The Netherlands, 1981; pp 331-342. (34) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (35) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (36) Guillot, B. J. Mol. Liq. 2002, 101, 219. (37) Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005, 122, 114507. (38) Hess, B.; van der Vegt, N. F. A. J. Phys. Chem. B 2006, 110, 17616. (39) Shirts, M. R.; Pande, V. S. J. Chem. Phys. 2005, 122, 134508. (40) Raschke, T. M.; Levitt, M. J. Phys. Chem. B 2004, 108, 13492. (41) Schravendijk, P.; van der Vegt, N. F. A. J. Chem. Theory Comput. 2005, 1, 643. (42) Cooper, T. G.; de Leeuw, N. H. Langmuir 2004, 20, 3984. (43) Delle Site, L.; Sebastiani, D. Phys. ReV. B 2004, 70, 115401. (44) Delle Site, L.; Alavi, A.; Abrams, C. F. Phys. ReV. B 2003, 67, 193406. (45) Hess, B. J. Chem. Phys. 2002, 116, 209.