Monte Carlo Simulation of the Adsorption of Phenol on Gold Electrodes . A Simple Model

Rodrigo S. Neves, Artur J. Motheo, Fernando M. S. Silva Fernandes* and Rui P. S. Fartaria a Instituto de Química de São Carlos, Universidade de São Paulo, Av. do Trabalhador Sancarlense, CP 780, 13560-970 São Carlos-SP Brazil Laboratory of Molecular Simulation and CECUL, Department of Chemistry and Biochemistry, Faculty of Science, University of Lisboa, Rua Ernesto de Vasconcelos, Bloco C8, 1749-016 Lisboa, Portugal


Introduction
The acquisition of information at molecular level, for the adsorption of organic molecules on electrochemical interfaces composed by metalic electrodes in contact with a liquid phase, is one of the main scopes of interfacial electrochemistry.Structural factors play a major role in the processes taking place in the interfacial region determining, for example, the resultant products of electron transfer reactions. 1,23][14][15][16][17][18][19] The combination of such kind of information with experimental results strongly contribute to the understanding of the behaviour of interfacial systems.
Phase transitions in adsorbed layers is another topic of interest concerned with the main adsorption processes on metalic electrodes. 20][23] The objective of the present work is to simulate, by canonical Monte Carlo method, the adsorption of phenol in a dilute aqueous solution on gold electrodes, at the potential of zero charge (pzc), and to qualitatively analyse how far we can go into the analysis of adsorption with a rather simple model.In an experimental investigation, Lezna et al. 24 concluded that, at pzc, the phenol molecule adopts a parallel orientation relative to electrode surface.However, when the potential across the interface is shifted to more positive values, this is followed by a reorientation of the molecule to a perpendicular orientation, the oxygen atom pointing to the surface and presenting a strong covalent interaction with the metal.Moreover, studies of the adsorption of phenol and correlated compounds on platinum show that for dilute solutions the same behaviour is observed, but for more concentrated ones, the quasiperpendicular orientation is always the preferred configuration. 25,26Thus, the final orientation of the adsorbate relative to the electrode seems to depend on the solution concentration.
Two main aspects of the process are considered in the present study: the approximation of the phenol molecule to the surface and its behaviour thereafter.To analyse the contribution of the solvent, potentials of mean force (PMF) were calculated for determining the preferential orientation of the phenol molecule when it approaches the electrode surface.In the next section, we present the essential aspects of the model, the interaction potentials and some simulation details.The results are, then, shown and discussed.The final section contains the conclusions of this work.

Model
The electrochemical cell The system was composed by a liquid phase, with 245 molecules of water and one molecule of phenol, encapsulated between two flat gold surfaces perpendicular to the z-direction.The use of one single phenol molecule corresponds to the condition of infinite dilution.The cell was of a rectangular shape, with 15.52 x 15.52 Å in the xy plane parallel to the metal surface, and 31.04Å in the zdirection.Periodic boundary conditions (PBC) were applied in the the x and y directions to replicate the cell in the directions parallel to the electrodes.PBC where not applied on the z-direction, which was confined by the electrodes.The size of the cell was chosen to reproduce the water density at the specified conditions of the simulations: 298.15K and 1 atm.
The use of flat gold surfaces does not take into account the effects of the atomistic corrugation of the electrodes.9][30] Despite being a simple approximation, we shall see in the next section that the present model gives results in accordance with experimental evidence.Moreover, it paves the way to more elaborate models and simulations.
Equilibration runs of 10,000 Monte Carlo cycles were followed by production runs of 100,000 cycles.One MC cycle consisted of random translations and rotations applied to each molecule of the system with an acceptance rate of approximately 50%.Several initial configurations were used for the equilibration runs.The production calculations started from the last configurations of the corresponding equilibration runs.

Models for the liquid phase
The interaction between the components of the liquid phase was modelled by Lennard -Jones (LJ) potentials between the sites of each molecule, combined with a coulombic interaction related to the partial charges on each site.
The TIP4P interaction model 31 was used for the water molecules.It consists of four interaction centers, related to the three atoms of the water molecule and to the electronic charge density due to the electronic pair on the oxygen atom.The later, is located at 0.15 Å from the oxygem and between the two hydrogen atoms.For the phenol, the potential proposed by Mooney et al. 32 was chosen.All the liquid phase molecules were considered as rigid species.The TIP4P model is rigid by its own definition.For phenol model, however, the internal degrees of freedom were constrained to the minimal energy molecular geometry. 32The potential parameters and site charges for each model are presented in Table 1, where ε and σ are the common Lennard-Jones parameters and q, the partial charge on each atom.The interaction site e, for the water molecule, represents the pair of electrons of the oxigem atom.The indexes indicated for the phenol atoms are related to Figure 1 The calculations were conducted by applying a cutoff radius of 7.5 Å to the LJ interactions with the usual long-range corrections.As for the electostatic components, they were evaluated over the entire simulation box without the consideration of long-range contributions.This last approximation introduces, of course, some discontinuities in the energy profiles, but, as we shall see, the simulation results seem to be, on the whole, consistent with the experimental evidence.Nonetheless, Ewald's sum, or the reaction field method, should be applied in a full quantitative study, especially when the sample has a small size.However, as mentioned before, the objective of the present work is to qualitatively analyse how far we can go with the simplest plausible model.In a more elaborate work, presently in progress, we will assess the importance of the long-range electrostatic contributions on the results, as well as their dependence on the sample size.

Models for the liquid phase in contact with the electrodes
The interactions of the liquid phase molecules with the flat electrodes were modelled by the combination of a LJ potential, to describe the interaction with the nonconductive electrons of the metalic surface, with an electrostatic interaction between each molecule in the solution and its image charge projected onto the electrodes, to describe the interaction with the conductive electrons.
In the case of the water molecules, the 9 -3 LJ successfully used by Philpott and Glosli, 19 for the solvent -electrode interaction in the simulation of benzene adsorption on gold electrodes, was adopted.For the phenol -electrode interaction, a 12 -3 LJ potential was derived by Hautman and Klein's method 21 used in the simulation of selfassembled alkyl thiols monolayers on gold electrodes.Both potentials are described by equations 1 and 2 respectively: , z is the distance between the interaction center and the surface and z 0 is a limit approach distance particular for each center.The parameters ε and σ are Lennard-Jones potential parameters.
Those potentials result from the integration of the 12 -6 LJ potential over the half-space related to the electrode surface.The original result is the 9 -3 LJ, but for some systems, such potential fails on the description of the interaction, mainly on the dispersion coefficient parameters. 21Thus, the 12 -3 LJ potential was used for the phenol -electrode interaction instead of the 9 -3 LJ.The potential parameters related to equations 1 and 2 are listed on Table 2.Both the potentials were truncated at 15 Å.
The interactions of the liquid phase molecules with the conductive electrons of the metal were modelled by a self-image charge potential method without distance cutoff.This approach, also takes into account the permissivity discontinuity at the interfacial region.In the original version, each charge in the liquid phase interacts with the projections of all image charges, reflected inside the electrode.This total-image method for the simulation of the polarization effects in the metal has been commonly used for water near metalic electrodes and the adsorption of liquid phases on surfaces, [5][6][7]18,19 however, some comments should be made. In th case of water near a polarizable surface, the induction effect due to a single molecule can be very large.Nonetheless, simulation studies of pure water in the vicinity of metalic surfaces using the total-image charge approach have shown a high degree of cancellation of interactions due to the sum over all the images and molecules.10 Thus, when only the self-image interactions are handled for the water-electrode interaction many body cancellations may not occur and an eventual overestimation of the interaction potential may take place, as pointed out by Shelley et al. 11 The net cancellation itself could also result from an overestimation of the induction power of water molecules distant from the electrode and separated by several layers of molecules.Indeed, it seems reasonable that water molecules of the nearest hydration layers of the electrode present a higher degree of induction when compared with the distant ones.
In the simulations reported here the self-image approach was applied with the reflection plane coincident with the electrode surface.The structural properties are in accordance to those obtained by Phipott and Glosli 19 in the simulation of the adsorption of benzene on gold, using the same 9 -3 LJ potential and a total-image charge potential.This is worth mentioning due to the structural similarities of benzene and phenol.

Structural properties
The reduced density profiles, for the hydrogens and oxygen of pure water, relative to one of the electrodes, are presented in Figure 2 and their main numerical characteristics are shown in Table 3.
The peaks in region I, centered at 2.7 Å from the surface, show that in the first layer the water molecules, in contact with the electrode, mainly have the molecular plane parallel to the electrode surface, with just a few of them with the hydrogen atoms pointing towards the surface (see sub-peak III).
The ordered structure of the first layer is a consequence of the equilibrium between the water-electrode potential and the hydrogen bonds between water molecules.It is observed in all simulations of water near metal surfaces either for flat surfaces 7 or structured ones. 8,11,18A second solvent layer is observed in region II, centered at 5.6-5.8Å.The relative positions of the oxygen and hydrogen peaks indicate a small inclination of the molecular plane, with the hydrogen atoms pointing towards the solution bulk, but the molecular orientation is mainly parallel to the electrode and still reflects a high degree of organization.In region IV, between 8 and 11 Å , a small perturbation on the water distribution is also detected.This density fluctuation shows a weaker influence of the electrode at those distances.For distances from the electrode surface greater than 11 Å the presence of the electrode seems to be not relevant and the density profiles approach the reduced asymptotic value around 1.10.These structural properties, obtained for pure water in contact with a flat gold electrode surface at pzc, are in accordance with results reported by Phipott et al. 19 and Hautman et al. 7 Therefore, the watermetal potential adopted in the present model seems to be a suitable choice for the simulation of the phenol adsorption in aqueous solutions.
Figure 3 displays the density profiles for water in the presence of one phenol molecule (3a) and the mean position of the center of mass of the phenol molecule (3b).The main numerical values of the water distribution are in Table 4.
The phenol molecule adsorbs at a distance of about 3.20 Å from the electrode, inside the first layer of water molecules, with an angle of 89 ± 4° between the aromatic ring plane and the normal to the electrode surface, in the   so-called horizontal configuration.This agrees with the experimental results for the adsorption of phenol on gold electrodes. 24The presence of the organic molecule causes some perturbations in the solvent layers at contact with the electrode.The analysis of the peak areas in region I of Figures 2 and 3 indicate that the first solvent layer suffers a decrease of approximately 16% in the number of water molecules directly adsorbed on the surface, when phenol is adsorbed, with a small displacement of the mean distance of the solvent, from 2.7 to 2.8 Å.This decrease in the number of water molecules adsorbed on the electrode is local, related to the region of the surface where phenol adsorption takes place, corresponding, on average, to displacement of 4 water molecules per phenol molecule adsorbed.The preferential orientation of water molecules in the first layer still remains the same as observed for pure water.Region III is slightly affected by the phenol presence, in the same way as discussed for region I.
On the other hand, in the second solvent layer, the oxygen density is weakly perturbed by the presence of phenol, although the hydrogen density is pushed to the solution direction approximately by 0.30 Å.In region IV, the presence of the adsorbate gives rise to an increase of structural order, with a clear peak observed at 12 Å from the electrode.Despite the influence of the phenol presence in this region, a more quantitative analysis is inadequate due to the small size of the present model and it was postponed to future simulations.
The final and more stable simulated configuration for phenol adsorption is the one with the aromatic ring parallel to the surface in accordance with the experimental observation of Lezna et al. 24 for the adsorption of phenol on gold polycrystalline electrodes.However, the present simulations of the evolution of the phenol tilt angle clearly indicate that the process takes place in two steps, as can be seen in Figure 4.
First, the phenol molecule adsorbs on the surface by the oxygen atom with the ring in a quasi-perpendicular orientation, oscillating around an angle of 25 º with the normal to the surface.After approximately 25,000 Monte Carlo cycles, there is a reorentation of the molecule to a horizontal configuration, oscillating thereafter around an angle of 89 ± 4°.A qualitative view of the two distinct events can be observed in the snapshots of Figure 5.
As already mentioned, there are experimental evidences that phenol, and some similar compounds, adsorbs on metalic surfaces at horizontal configurations for low concentration solutions, but the same behaviour is not observed at high concentrations.In this situation, the adsorption persists in the quasi-perpendicular orientations with the negative center, the oxygen in the case of phenol, pointing towards the surface.The picture of the adsorption process from the simulations reported here, associated to experimental information, suggests that the approach of the phenol molecule to the surface, at the pzc, is always in a quasi-perpendicular orientation.In low concentration solutions, a reorientation occurs and the molecule reaches a final horizontal configuration.At high concentration solutions, the number of adsorbate molecules near the electrode surface is large and lateral interactions will presumably hinder a further reorientation to a final horizontal configuration.

Potential of mean force profiles
The structural properties are fundamental to understand the mechanism of the processes taking place in the interfacial region.Additionally, it is very important to know the driving forces underlying the adosrption mechanism described in the last paragraph.That is, the approach of the phenol molecule to the electrode in a quasiperpendicular orientation of the aromatic ring, followed by its reorientation to a parallel configuration relative to the surface.This can be analysed through the potential of mean force (PMF) acting on phenol as a function of the distance, z, measured from the phenol centre of mass and the electrode along the perpendicular to the surface.
The potential of mean force acting on phenol, at distance z, is defined as the reversible work necessary to approach the adsorbate molecule, in a pre-defined orientation inside the solvent, from a distance z 0 , where the influence of the electrode is negligible, to a nearer distance z.That is the Helmohltz free energy change taking into account the presence of the solvent.
The total free energy change, ∆A, has two contributions: the PMF due to the solvent and the potential of the direct force of the electrode on the adsorbate.The last is the the phenol-electrode interaction energy given in the previous model section.The solvent contribution, ∆A S , was calculated by the free energy perturbation (FEP) method, 10,15,17,33,34 whose general working equation, in the present context, is: ( where U(z i ) is the potential energy of a configuration of the reference system, U(z i + δz i ) the potential energy of a configuration of the perturbed system, k the Boltzman constant and T the absolute temperature.The brackets represent the ensemble average over an assembly of configurations, generated for the reference system, by randomnly translating and rotating the water molecules, but keeping the phenol molecule at the distance z i (measured from its mass centre to the electrode) and at a fixed orientation (perpendicluar or horizontal).The interactions taken into account at this stage are only the solvent-solvent, solvent-adsorbate and solvent-electrode interactions.
The assembly of configurations for the perturbed system are the same as the preceding ones as far as the water molecules are concerned, but now with the phenol molecule, with the same orientation, in the new fixed distance z i + δz i where δz i is a small displacement applied to the phenol molecule for generating the perturbed system.
As the reference and perturbed states only differ in the position of the phenol molecule they are unambiguously identified by the coordinate z.The term U(z i + δz i )-U(z i ) in equation ( 3) is just the difference between the phenolwater potential energies in the two states, since the contributions of the solvent-solvent and solvent-electrode, although essential to generate the reference assembly, are the same int the reference and perturbed systems.
In order to get a proper convergence of the free energy, that difference must be much smaller than kT.Therefore, the distance z-z 0 is divided into a number of windows, N w , each characterized by the coordinet z i and equation ( 3) is applied to them.The potential of mean force, due to the solvent, acting on the phenol molecule at distance z, is then given by: (4) The total free energy change is obtained by adding the phenol-electrode potential to equation (4).
For each window, equilibration runs of 15,000 MC cycles were followed by production runs of 10,000 cycles.The first window was chosen at a distance z 0 = 13.5 Å between the phenol mass center and the electrode surface.The following windows were located at successive displacements of 0.22 Å toward the surface.The perturbation distance, δz i , was 0.1 Å. and a double-wide sampling was used. 33Additionally, during the production runs, the phenol molecule was allowed to move in the x and y directions, however keeping the orientation, for a better relaxation of the solvent in the vicinity.
Two orientations of the phenol molecule were tested: case 1, with the aromatic ring parallel to the electrode and case 2 with the ring perpendicular to the surface and the oxygen pointing towards the metal.These are the most interesting geometries considering the structural properties showed in the last section.
The contribution of the solvent, the interaction of the phenol with the electrode and the total free energy change for cases 1 and 2 can be seen in Figures 6 and 7, respectively, with fluctuations bellow 4% when the free energy values are approximately constant and between 6% and 9% when it changes abruptely.
Region I of both Figures, at about 5.5 Å from the electrode, is the penetration onset of the phenol molecule in the first solvation layer of the electrode.In case 1, the solvent contribution to the free-energy is approximately 10 kJ mol -1 higher than in case 2. When the phenol molecule approaches the first water layer, with the ring parallel to the electrode, it finds firstly that small energy barrier which is presumably due to the compression of the solvent layer against the electrode surface.In case 2, at approximately the same distance, the phenol molecule experiences a small energy minimum.For this geometry, when the center of mass is at 5.5 Å from the electrode, the oxygen atom is at ± 2.8 Å, inside the first water layer and physically adsorbed on the surface.It seems, therefore, that the difference of 10 kJ mol -1 , in the solvent contribution to the free-energy, favours the perpendicular approach.
Additionally, the analysis of the remarkable energy barriers to the approach of the phenol molecule, for distances less than 5 Å, shows that they have different origins.In case 1, the barrier arises from the further compression of the solvent layer in the direction of the electrode, caused by the parallel approach of the ring.This will disable the adsorption of the phenol molecule directly in the horizontal configuration.In case 2, the energy barrier results from the direct interaction of the phenol oxygen with the electrode surface.In this situation, no significant structural changes of the first solvent layer occur and the phenol will adsorb by the oxygen atom, in a quasi-perpendicular configuration, performing the first adsorption step.Then, the transition to the horizontal configuration is driven by the interaction of the aromatic ring with the electrode (see Figures 6 and 7).Indeed, the horizontal configuration has an interaction energy of -40 kJ mol -1 , whilst the interaction energy for the perpendicular configuration is -15 kJ mol -1 .The transition process is smooth and allows the first solvent layer to successively relax, as was observed by the evolution of the tilt angle on Figure 4.

Conclusions
We have presented canonical Monte Carlo simulations of the adsorption of phenol, in a dilute aqueous solution, on flat gold electrodes at the potential of zero charge.
The results clearly show that the adsorption process occurs in two distinct and successive steps.The first one is the physical adsorption of the phenol oxygen and the aromatic ring pointing towards the bulk of the solution.This is followed by a second step consisting of the reorientation of the molecule to a horizontal configuration, at a distance of approximately 3.20 Å from the surface, and a mean angle of 89 ± 4° with the axis normal to the surface.The final configuration is in agreement with the experimental studies of Lezna et al.. 24 The potentials of mean force indicate that the perpendicular approach is favorable to the penetration of the organic molecule across the first solvent layer on the electrode.A horizontal approach would result in the compression of that layer and in a considerable energy barrier obstructing the adsorption process.For the perpendicular approach, the barrier energy is mainly due to the direct interaction of the phenol oxygen with the electrode surface.
Experimental evidence, obtained from the study of the adsorption of aromatic molecules on platinum electrodes, 25,26 show that for low concentration solutions the adsorption process ends in a horizontal orientation, but at high concentrations the quasi-perpendicular orientation is the preferred one.The behaviours observed in the present work for the adsorption on gold, suggest that in the case of a greater number of phenol molecules on the surface, lateral interactions will difficult the reorientation and the perpendicular adsorption should indeed persist.
Despite the simplicity of the model, the structural properties and the consistency of the respective thermodynamic behaviour, suggest that the present model represents the main characteristics of the physical adsorption of phenol on gold, giving insight into the microscopic details.
The reported results constitute the starting point to a more refined model, presently under development, in which we consider the structural properties of the electrode surface to analyse the preferential top and hollow adsorption sites and their specific influence in the adsorption process.The introduction of intramolecular degrees of freedom and the developement of accurate ab-initio potential surfaces, mainly for the interactions between the molecules and gold, as well as the application of long range corrections for the electrostatic potential, as the Ewald sum, are also of the utmost importance regarding the full understanding of the whole process.
Finally, a current major challenge in computational studies of complex systems is to solve the multi-scale problem (in time and in space) without requiring an unrealistic time consuming computing effort.Work along these lines is in progress and will be reported soon.

Figure 2 .
Figure 2. Density profiles (ρ) of pure water/electrode as a function of the distance (z) from the electrode.

Figure 3 .
Figure 3. a) Density profiles (ρ) of water/electrode in the presence of one adsorbed phenol molecule as a function of the distance (z) from the elctrode.b) Average position of the mass center of the adsorbed phenol molecule.

Figure 4 .
Figure 4. Evolution of the tilt angle between the aromatic ring of the adsorbed phenol molecule and the axis normal to the electrode surface during the MC simulation.

Figure 5 .
Figure 5. Snapshots of the phenol molecule adsorbed on gold electrode.a) Quasi-perpendicular configuration.b) After the reorientation to the horizontal configuration.

Figure 6 .
Figure6.Free energy changes associated to the phenol approach to the surface in the horizontal orientation.

Figure 7 .
Figure 7. Free energy changes associated to the phenol approach to the surface in the perpendicular orientation.

Table 1 .
Interaction parameters for all the atoms on the liquid phase Representation of the phenol molecule.

Table 2 .
Interaction parameters for all atoms on the liquid phase with the surface α only hydrogen atoms bonded to carbon.

Table 3 .
Water density profiles characteristics

Table 4 .
Water density profiles characteristics in the presence of phenol