Liquid-Phase Effects on Adsorption Processes in Heterogeneous Catalysis.

Aqueous solvation free energies of adsorption have recently been measured for phenol adsorption on Pt(111). Endergonic solvent effects of ∼1 eV suggest solvents dramatically influence a metal catalyst's activity with significant implications for the catalyst design. However, measurements are indirect and involve adsorption isotherm models, which potentially reduces the reliability of the extracted energy values. Computational, implicit solvation models predict exergonic solvation effects for phenol adsorption, failing to agree with measurements even qualitatively. In this study, an explicit, hybrid quantum mechanical/molecular mechanical approach for computing solvation free energies of adsorption is developed, solvation free energies of phenol adsorption are computed, and experimental data for solvation free energies of phenol adsorption are reanalyzed using multiple adsorption isotherm models. Explicit solvation calculations predict an endergonic solvation free energy for phenol adsorption that agrees well with measurements to within the experimental and force field uncertainties. Computed adsorption free energies of solvation of carbon monoxide, ethylene glycol, benzene, and phenol over the (111) facet of Pt and Cu suggest that liquid water destabilizes all adsorbed species, with the largest impact on the largest adsorbates.

Introduction

A heterogeneously catalyzed process starts with the adsorption of at least one of the reactants and ends with the desorption of the products.[1] Activation barriers have also been correlated to adsorption energies for about 90 years through, e.g., the well-known Brønsted–Evans–Polanyi (BEP) relationship that connects thermodynamics with kinetics.[2] Thus, understanding adsorption and desorption processes is critical for heterogeneous catalysis, and adsorption processes at gas–solid interfaces have been studied for many decades.[3−7] For solid–liquid interfaces, less progress has been made, largely due to the added complexity of the solvent molecules that can both directly and indirectly interact with the adsorbate through, e.g., hydrogen bonding and by changing the electronic structure of the adsorbent. Nevertheless, it is well-known that solvents can affect important figures of merit in catalysis, such as activity and selectivity, that go beyond reduced mass transfer rates and nonideality effects that determine a reactant’s activity/fugacity at a given liquid concentration.[8−11] These intrinsic solvent effects originate from one or a combination of the following: nonrandom distribution of solvent and solute molecules at interfaces and the nature of the solvent (pH, polarity, H-bonding ability, etc.)[12−14] leading to changes in the relative stabilization of adsorbed species,[15−17] participation of the solvent in the catalytic cycle or reaction coordinate,[18] and competitive adsorption between the solvent and adsorbate molecules.[19−21] Consequently, to design catalysts for enhanced activity, selectivity, and stability, it is critical to understand the root causes of solvent effects on adsorption processes in heterogeneous catalysis. Such an understanding becomes more important the larger the solvation effects are and the more the solvation effects vary for the various species and transition states in the catalytic system (changes in trends). Recently, aqueous solvation effects have been measured on the free energy of adsorption for phenol, benzyl alcohol, benzaldehyde, and cyclohexanol over a Pt catalyst at 298 K[22] and for phenol over Pt and Rh catalysts at four temperatures (283, 288, 298, and 314 K).[23,24] Importantly, for phenol adsorption on supported Pt catalysts, it has become possible to deconvolute the adsorption isotherm data for the different surface facets, which enabled the measurement of the aqueous solvent effect on the phenol adsorption free energy on the Pt(111) facet, which again facilitates a correlation to computational studies. Experimentally measured aqueous solvation effects are large (>1.0 eV at 298 K) and endergonic for phenol on Pt(111), suggesting that the kinetic properties of Pt nanoparticles for phenol catalysis are a strong function of solvent properties. In addition, recent simulations suggest that solvation effects are not only a strong function of adsorbate but also of the nature of the catalyst surface (e.g., Pt vs Cu).[25] Solvation effects do not cancel out in a descriptor-based catalyst design strategy, and solvation effects need to be considered explicitly. Differences in 1 eV in adsorption free energy are just too large as they could be neglected in any descriptor-based catalyst design strategy. We highlight that at 298 and 500 K, a change in a rate-controlling free-energy barrier of 0.1 eV can change the overall rate by factors of 50 and 10, respectively. However, measurements are indirect and involve adsorption isotherm models, which potentially reduces the reliability of the extracted energy value. Also, computational, implicit solvation models predict significantly smaller and exergonic solvation effects on the free energy of phenol adsorption on Pt(111), failing to even qualitatively agree with the experimental measurements.[23]

Despite the availability of fairly established methods for computing the energetics of adsorption reactions at gas–solid interfaces, methods capable of quantifying the adsorption (free) energies at solid–liquid interfaces are less developed. This is partly due to the limited availability of experimental[21,26−28] and theoretical studies[29−41] and partly due to inherent intricacy of a reaction system comprised of both a complex heterogeneous catalyst and a condensed phase. Although “on the fly” electronic structure calculations in ab initio molecular dynamics (AIMD)[42−44] simulations are well-suited for this task, as a result of the large computational cost associated with quantum mechanical calculations and configuration space sampling, AIMD simulations for computing adsorption free energies from a condensed phase are currently impractical (with small confidence interval from configuration space sampling). Typically, the high computational cost of AIMD constrains both the size of the simulation system to a few hundred atoms and the time scale of simulations to a few picoseconds.[45−47] Alternatively, implicit solvation schemes[48,49] and classical force field simulations offer a practical approach that is computationally affordable. However, force field simulations require a reliable potential for all fluid components with the adsorption site/metal surface, which, except for perhaps water, is hardly available due to a lack of reliable experimental data. Similarly, performance reliability is largely unknown for adsorption on metal surfaces with implicit solvation models due to the need to parameterize the implicit solvation models against experimental data that are again hardly available or possess unknown error bars. We highlight that for homogeneous systems, implicit solvation models have been exceptionally successful despite their limitations in capturing anisotropic site-specific interactions (e.g., hydrogen bonding).[50,51] The success of the implicit solvation models for homogeneous catalysis studies (relative to heterogeneous catalysis) can be understood by the metal site being surrounded by ligands in homogeneous catalysts, i.e., the metal site is not directly exposed to solvent molecules, while in heterogeneous catalysts, there are no ligands shielding the metal site from the solvent molecules.

To overcome these challenges for heterogeneous catalysis, a number of research groups have recently developed combined quantum mechanical/molecular mechanical (QM/MM) approaches for predicting free-energy changes for processes on solid–liquid interfaces.[34,36,37,52−54] In this class of simulations, the essential region of interest, that can often not be properly described by classical force fields, is described quantum mechanically, while the remainder of the system is described at the MM level of theory, enabling a balance of computational cost and accuracy. For example, in a metal-catalyzed adsorption reaction, adsorbate, active site, and its immediate metal neighbors are treated from first principles (called QM system/subsystem in this study), while the bulk of the solvent molecules and the nonreactive part of the simulation system are treated using classical force fields (called MM system/subsystem in this study). In this regard, we previously developed a computationally affordable and reliable method for computing free-energy changes in the presence of solvents using a hybrid QM/MM approach with electrostatic embedding, named Explicit Solvation for Metal Surfaces (eSMS)[55] and applied it to several surface reactions over various transition-metal surfaces.[25,56]

A limitation of this previously developed eSMS has been that it can only be applied for computing free-energy differences of surface reactions; that is, it can practically not be employed to compute free energies of adsorption processes from a liquid-phase environment due to (i) the long reaction coordinate in adsorption processes that increases the computational expense, and (ii) the need to be able to compute the free energy of the quantum mechanically described part of the simulation system within the harmonic approximation, i.e., only the classically described part of the simulation system is sampled extensively within eSMS. In this study, we identify a reaction coordinate for adsorption processes that enables our eSMS approach to predict adsorption free energies at a solid–liquid interface. Next, we validate our novel computational scheme against previous eSMS calculations of the aqueous solvation effect on the activation free-energy barrier of O–H bond cleavage of ethylene glycol, EG (C2H6O2), over a Cu(111) catalyst surface model[25] by computing the adsorption free energy of both the reactant and transition state. Next, we compute the aqueous solvation effect on the free energy of phenol adsorption on the Pt(111) surface, and we reanalyze the experimental data for this solvation free energy using multiple adsorption isotherm models. Computational predictions at 298 K are in good agreement with the experimental data with remaining errors likely originating from errors in the metal–water interaction potential. Thus, these results confirm a large endergonic aqueous solvation free energy for phenol adsorption with significant implications for the catalyst design for the conversion of aromatic molecules. To elucidate the significance of these large, computed and measured, solvation effects for the Pt–phenol system for various other catalytic processes, we computed the aqueous solvation effect on the adsorption free energies of benzene over Pt(111) and carbon monoxide (CO), ethylene glycol (EG), and phenol (Ph) over Pt(111) and Cu(111) surfaces.

QM/MM-FEP Scheme for Adsorption Processes at a Solid–Liquid Interface

The key challenge for computing the free energy of adsorption in a liquid-phase environment with a hybrid quantum mechanical/molecular mechanical free-energy perturbation (QM/MM-FEP) scheme is identifying an efficient thermodynamic cycle. Specifically, to reduce the computational cost, the harmonic approximation must remain meaningful for the quantum mechanically described subsystem in the reactant and product states, and extensive configuration space sampling can only be necessary for the classically described part of the simulation system. Furthermore, the reaction coordinate must be relatively smooth and not involve overcoming large free-energy barriers, such that, only few (computationally expensive) free-energy perturbation (FEP) steps are required for computing the overall free-energy difference. Figure illustrates such a thermodynamic cycle that is explained below. A detailed computational methods section is presented below in the Methodssection of this article. The key idea of our thermodynamic cycle is to not compute the desorption of an adsorbate into a liquid phase but to only compute the solvent effect on the desorption process. In this way, free-energy differences for the overall process are smaller, and desorption itself is only the gas-phase desorption, which can be described using the harmonic approximation or correction formulas to the harmonic approximation.[57,58]

Figure 1

Free-energy perturbation approach for computing the difference in the adsorption/desorption free energy of any adsorbed species in the presence of liquid water (blue circle solid (g) + *(l) ↔ purple circle solid*(l)) and in a gas phase (blue circle solid(g) + *(g) ↔ blue circle solid*(g)), i.e., for desorption: purple circle solid*(l) + *(g) ↔ *(l) + blue circle solid *(g). Between state Iopt,0 and Iopt, the coordinates of the metal surface in the gas-phase system change to those metal coordinates in a system with adsorbate. There is no change in the liquid-phase system. Between state Iopt and I, the coordinates of the adsorbate atoms described quantum chemically in the liquid-phase system transition from those optimized in the liquid to those optimized in the gas phase. There is no change in the gas-phase system. Between states I and II, the (classical) electrostatic interaction between the solvent and the adsorbate is slowly removed in the liquid-phase system. No coordinates of atoms described quantum chemically change. Intermediate states, such as state Ia, are “nonphysical” states, which is typical for FEP calculations. From state II to III, the nonclassical contributions of water molecules changing the electronic structure of the quantum system transition from those of a metal–adsorbate system to those without the adsorbate, i.e., this energy difference constitutes the main difference between electrostatic and mechanical embedding of the QM/MM system. Green circle solid symbolizes the adsorbate being only present in the classical DL_POLY simulation and not in the quantum system. van der Waals interactions between the solvent and the adsorbate are slowly removed between states III and IV. No coordinates of atoms described quantum chemically change. Intermediate states, such as state IIIa, are again “nonphysical” states. From state IV to state IVopt, the surface atom coordinates in liquid transition to the optimized coordinates of a free site. There is no change in the gas-phase system. The gas-phase desorption is described between states IVopt and V. The (free) energy difference between states IVopt and Iopt,0 is the solvent effect on the stability of the adsorbate.

To facilitate understanding of our QM/MM-FEP scheme, we first summarize here the physical meaning of the different steps shown in Figure before discussing all steps in more detail in the forthcoming paragraphs. We start from a state containing two systems: a system with adsorbed species in a liquid with the QM subsystem coordinates optimized in the liquid and a second (clean) gas-phase optimized adsorbent system in vacuum (Iopt,0). The second, (clean) gas-phase optimized adsorbent system is identical in both the initial and final states (V). The final state (V) consists of the (clean) gas-phase optimized system, an adsorbate in the gas phase, and a (clean) liquid-phase system with QM system coordinates optimized in the liquid. Thus, the free-energy difference between states V and Iopt,0 is the free energy of adsorption from the gas phase and a (clean) liquid-phase system to an adsorbed state in liquid. The presence of the second (clean) gas-phase optimized system might appear confusing at the beginning; however, it will facilitate the overall free-energy calculation, specifically step II-to-III discussed below. First, we change the gas-phase optimized (clean) adsorbent coordinates in vacuum to those coordinates of the adsorbent when the adsorbate is adsorbed on the adsorbent in the gas phase (Iopt). There is no change in the liquid-phase system such that this is technically not an FEP calculation. Second, we transform the liquid-phase optimized system to a state I in liquid where the QM subsystem coordinates are those of a gas-phase optimized system. There is no change in the second gas-phase system; however, at the end of this transformation, the adsorbent atom coordinates are equivalent in both systems. This free-energy difference considers the fact that adsorbate and adsorbent coordinates are somewhat different in a liquid and gas-phase environment. Third, the classical electrostatic interactions between the adsorbate (and neighboring adsorbent) atoms and the water molecules are removed to reach to state II. Again, there is no change in the second gas-phase system. Fourth, the polarization of the adsorbent–adsorbate subsystem, induced by the presence of water molecules (point charges), is changed to a polarization of the adsorbent system in the absence of the adsorbate to reach state III. This step involves removing the adsorbate from the QM part of the QM/MM calculation for the liquid-phase system and adding the adsorbate to the second gas-phase system. The QM atom coordinates are equivalent in both systems, which facilitates this FEP calculation that involves no intermediate states. This step is different from the removal of the classical electrostatic interactions as it considers the effect of the point charges on the quantum chemically described subsystem. We call this step removal of nonclassical electrostatic effects, and hence, all electrostatic interactions are removed at state III. Fifth, the Lennard-Jones interactions are gradually removed in the liquid-phase system to reach state IV. There is no change in the second gas-phase system. In the liquid-phase system of state IV, there are no interactions between the adsorbate and the liquid molecules, and the adsorbate is exclusively in the second gas-phase system. However, the adsorbent coordinates in the liquid-phase system, that surrounded the adsorbate, are still those that were deformed by the adsorbate. Hence, in the sixth step, the adsorbent coordinates are transformed to those of the liquid-phase optimized adsorbent without adsorbate (state IVopt). Again, there is no change in the second gas-phase system that contains the adsorbate. Finally, the adsorbate is removed to its free state in the second gas-phase system (state V). There is no change in the liquid-phase system, and this is a typical gas-phase desorption calculation. A more detailed explanation of these steps for computing solvent effects for adsorption/desorption on a metal surface is provided in the following.

Step Iopt,0 to Iopt in Figure

Between states Iopt,0 and Iopt, the coordinates of the metal surface in the gas-phase system (optimized free site) change to those metal coordinates in a system with adsorbate. This is often a significant change in energy, but it does not involve the liquid-phase system and hence, no FEP calculation is involved in this change.

Step Iopt to I in Figure

The QM atom coordinates of the metal–adsorbate cluster system in liquid are optimized in liquid in state Iopt,0, while they are those optimized in the gas phase in state I. The free-energy change for this transition in QM atom coordinates involves standard eSMS free-energy calculations that involve a number of QM and MM calculations that are reasonably fast given the small change in QM atom coordinates during adsorption of typical adsorbates on a metal surface.[8,25,55,56] We note that for simplicity, only the coordinates of the adsorbate atoms (not metal atoms) described quantum chemically in the liquid-phase system transition from those optimized in the liquid to those optimized in the gas phase. The change in metal atom coordinates was found to be too small to justify the required computational effort. Specifically, the total potential energy function of our eSMS method is given in eq , where the first term constitutes the conventional gas-phase energy computed for a periodic slab model and the combination of the second to fifth terms is the nonclassical part of the electrostatic interaction energy of a mean-field liquid-phase environment with the adsorbate and metal surface described quantum mechanically by a cluster model. We have previously shown that unlike gas-phase adsorption energies on metal particles, this nonclassical part of the electrostatic interaction energy is fast converging with the metal-cluster size such that reliable energies can be computed for relatively small metal clusters of ∼50 metal atoms. The last term accounts for the classical (electrostatic and van der Waals) interaction energy within the MM subsystem and its interaction with the QM subsystem. The full derivation of eq can be found elsewhere.[59]In eq , the first term, EQMSurface (r̲QM), is evaluated for a periodic slab in vacuum, using the VASP[60,61] program package, and the second term, EQMcluster(r̲QM), is evaluated for a QM cluster in vacuum, with the help of the TURBOMOLE program package.[62−64] The combination of the third–fifth terms is a QM calculation in a mean field of MM water molecules using the periodic electrostatic embedded cluster method (PEECM)[65] with the fixed charge approximation (the fixed charge approximation has been validated for our eSMS approach in ref (59)). We note that the number 100 in the equation indicates that 100 MM conformations, selected equally spaced from 1 ns molecular dynamics (MD) simulations (10 ps apart), were used to represent the mean field of the MM water molecules. We also note that subscripts QM, MM, and QM/MM indicate QM subsystem (metal cluster and adsorbate), MM subsystem (e.g., TIP3P water molecules and the rest of the metal atoms), and the interface of the QM and MM subsystems, respectively (this also applies to eq ). A rationale for using the various periodic and nonperiodic density functional theory (DFT) codes is given in the Supporting Information (SI).

Step I to II in Figure : Classical Electrostatic Interaction Removal

In our QM/MM potential energy function of the reaction system, the electrostatic interactions between the adsorbed species and the metal atoms in direct vicinity of the adsorbate and the solvent molecules are governed by the magnitude and distance of the partial charges. To remove the electrostatic potential, first, partial charges on the “metal–adsorbate cluster” (see Figure S7) and “metal cluster” are computed using a suitable charge model (we used NPA[66] and DDEC6[67] in this study—see SI for a more detailed explanation). We note that the metal cluster is a “nonphysical” state in which the adsorbate is present but does not electrostatically interact with the surrounding water molecules that are treated classically, i.e., the partial charges on the adsorbate are set to zero (Qads = 0). Thus, to calculate the charges on the metal atoms in the “metal-cluster” system, the adsorbate atoms are removed from the “metal–adsorbate cluster” system during the charge calculation. Next, we insert adequate “nonphysical” intermediate states (such as state Ia in Figure ) in between by linearly reducing the magnitude of partial charges on the QM atoms going from the metal–adsorbate to the metal cluster. Table S3 summarizes the charges on the QM atoms of these two states for different adsorbates studied in this work.

Since the QM atom geometry remains the same during the removal of the electrostatic interaction, all terms in the potential energy function in eq , except the classical EMM+QM/MMelec+vdW(r̲QM,r̲MM) term, remain the same in every intermediate state and hence, they do not need to be computed in an FEP step calculation. In other words, removing the classical electrostatic interaction involves no QM calculations but only classical free-energy calculations.

Step II to III in Figure : Nonclassical Electrostatic Effects due to Changes in Polarization by the Solvent

The nonclassical contribution of the solvent molecules changing the electronic structure of the quantum system transitions from those of a metal–adsorbate system to those without the adsorbate, i.e., this energy difference constitutes the main difference between electrostatic and mechanical embedding of the QM/MM system. The energy difference is given bywhere the metal-(ads)-solvent and metal-ads-solvent superscripts represent states III and II, respectively (the first and last terms in parenthesis are performed using the TURBOMOLE program package). This step involves a QM calculation; however, the free-energy change calculation is not a strong function of solvent coordinates (only the computational description of a “nonphysical” state changes) and although the free-energy differences are significant, multiple independent calculations suggest that this free-energy change can efficiently be computed without introducing intermediate states (which would be impractical). Here, we also note that given the magnitude of the free-energy change from state II to III, electrostatic embedding is superior to mechanical embedding (that does not consider these nonclassical electrostatic contributions) for QM/MM-FEP calculations of metallic systems where electrons move freely.

Step III to IV in Figure : van der Waals Interaction Removal

Another interaction between the solvent and the adsorbate is the classical van der Waals (vdW) interaction. To remove the vdW interaction, we followed the classical approach of scaling of the distance-shifted Lennard-Jones (LJ) potentials.[68] This approach allows for a smooth transition in molecular simulations between real and dummy atoms for which all atomic interactions have vanished. The LJ potential that we implemented in the DL_POLY program package to annihilate or create an atom has the following functional formwhere λ is the coupling parameter in the FEP calculation, r is the interatomic distance between two atoms, and A and B are repulsive and attractive LJ parameters, respectively, calculated as a product of LJ parameters specific for each atom type of the interacting particles (A = 4εσ12, B = 4εσ6; see Table S1 for ε and σ values we used in this work). For λ = 0, the vdW interactions fully exist and an increase in λ leads to a smooth disappearance of the vdW interactions. Finally, at λ = 1, the interaction between the solvent molecules and the adsorbate is zero everywhere, corresponding to the final state where the adsorbate atom is transformed to a dummy particle. The shift parameter δ comes into play for λ larger than zero. It allows for a smooth transition from the original LJ potential to no interaction. As recommended by Zacharias et al.,[68] the parameter δ was chosen as the square of the vdW radius of the interacting atoms to allow the smoothest transition between an atom present on the surface and filling the cavity after molecule annihilation. Table S4 summarizes all δ parameter values for all pair interaction terms and, as an example, Figure S2 illustrates the LJ potential as a function of the distance between the oxygen atom of an adsorbed CO and the oxygen atom of a (TIP3P) water molecule for different values of λ.

To enable a reliable free-energy difference calculation, the removal of the vdW interactions between the solvent and the adsorbate is done in multiple steps by linearly increasing λ from 0 to 1. We note that these intermediate states, such as state IIIa (see Figure ), are “nonphysical” states. Analogous to the removal of the classical electrostatic interactions, in the potential energy function, eq , only the Eelec+vdW(r̲QM,r̲MM) term changes between two consecutive states and no QM calculations are required for this step, i.e., only classical free-energy calculations are performed at the MM level of theory. At this point, all classical interactions between the adsorbate and the solvent have vanished.

Step IV to IVopt in Figure

Given that the quantum atom coordinates of the metal cluster are still those in the gas-phase adsorbed state (obtained by removing the adsorbate atoms of the metal-cluster adsorbate in state I) and not those in the gas-phase optimized metal cluster (free site), the free-energy change for this transition in QM atom coordinates needs to be computed. This again involves standard eSMS free-energy calculations.[8,25,55,56] We note that for all free-energy calculations, we employ the Bennett acceptance ratio (BAR) as the free-energy estimator for all FEP steps.[69]

Step IVopt to V in Figure : Gas-phase Desorption

We use for the adsorbate as reference state, the adsorbate in a gas phase in equilibrium with the adsorbate in the liquid phase. Thus, only a typical gas-phase metal adsorption calculation at the QM level of theory is required for this step.

The free-energy change from state Iopt,0 to state V constitutes the free-energy change from an adsorbed molecule surrounded by a solvent to the adsorbent (metal surface) surrounded by a solvent and the adsorbate in a gas-phase state. Thus, the free-energy change from state Iopt,0 to state IVopt constitutes the solvent effect on a desorption process (ΔΔGAdsorbateliq→gas(desorption) = −ΔΔGAdsorbategas→liq(adsorption) and ΔΔGAdsorbategas→liq = ΔGAdsorbateg–l – ΔGAdsorbategas) that is likely not as DFT functional dependent as the gas-phase adsorption energy.[70] Also, the activity/fugacity calculation of the adsorbate in gas phase is often significantly easier than such a calculation in a solvent. In this paper, we assumed the gas phase to behave as an ideal gas. For adsorbates with very low gas pressure, where it is challenging to measure the fugacity, the gas–liquid partition coefficient (equilibrium constant) can be computed from COSMO-RS calculations that are quite reliable for this task.

Results and Discussion

Numerical Verification of the Proposed Thermodynamic Cycle

To numerically verify the proposed methodology for adsorption/desorption processes, the aqueous-phase effect on the free energy of activation (ΔΔGliq‡ = ΔGliq‡ – ΔGgas‡) of the O–H bond cleavage of ethylene glycol over a Cu(111) surface at 423 K is computed using two different approaches, the eSMS methodology for surface reactions[25] and the adsorption scheme proposed here. In our proposed adsorption scheme, the solvent effect on the activation barrier (ΔΔGliq‡) is given by eq where ΔΔGTSgas→liq and ΔΔGReactantgas→liq indicate the solvent effect on the adsorption of the transition state (TS) and reactant state (RS), respectively.

Figure a illustrates the free-energy profile of the solvation effects on the desorption of the RS and TS (ΔΔGAdsorbateliq→gas(desorption) = −ΔΔGAdsorbategas→liq (adsorption)). For both the RS and TS of the O–H bond cleavage of EG on the Cu(111) surface at 423 K, the transition of the QM atom coordinates (Iopt,0-to-I) leads only to a small free-energy change relative to the overall solvent effect. The transition state is significantly stabilized (ΔΔGTSgas→liq = −0.48 ± 0.03 eV, ΔΔGRSgas→liq = 0.14 ± 0.01 eV; see Table S5), while the reactant state is destabilized in liquid water (we note that free-energy profiles in Figure are plotted for the desorption process) leading to a significantly reduced activation free-energy barrier. Both the RS and TS are similarly destabilized through the creation of a solvent cavity on the surface as shown in the change in free energy during the removal of the van der Waals/Lennard-Jones interactions (see Table S5). In contrast, the classical electrostatic interactions stabilize the RS and TS in solution (ΔΔGEL-Rmv, RSliq→gas = 0.49 eV, ΔΔGEL-Rmv, TSliq→gas = 1.25 eV). The difference in classical stabilization can be attributed to the significant change in partial charges of the QM system going from RS to TS (see Table S3a). In particular, the charge of the cleaved hydrogen atom has transitioned from 0.47 e in RS to −0.16 e in TS, explaining the 0.76 eV difference in solvent stabilization due to classical electrostatic interactions. Overall, the electrostatic interaction effect is, however, altered by the nonclassical electrostatic interactions from section “II-to-III”. The latter effect is ignored in molecular embedding calculations (ΔΔGII-to-III, RSliq→gas = −0.05 eV, ΔΔGII-to-III, TSliq→gas = −0.21 eV). Therefore, we conclude that electrostatic embedding of the QM/MM system is crucial in any hybrid QM/MM approach for metallic reaction systems.

Figure 2

(Average) free-energy profiles for aqueous-phase effects on the desorption (ΔΔGliq→gas = −ΔΔGAdsorbategas→liq) of (a) reactant and transition states of O–H bond cleavage of ethylene glycol over a Cu(111) surface at 423 K, (b) CO molecule over Cu(111) and Pt(111) surfaces at 298 K, (c) ethylene glycol molecule over Cu(111) and Pt(111) surfaces at 423 K, and (d) phenol molecule over Cu(111) and Pt(111) surfaces at 298 K. ΔΔG (ΔΔGAdsorbateliq→gas) is the aqueous-phase effect on the low coverage desorption of an adsorbate. See Figure for more information about the physical meaning of states Iopt,0, I, II, III, IV, IVopt. “Remove electrostatic” and “remove Lennard-Jones” represent the removal of classical electrostatic and Lennard-Jones interactions between an adsorbate and the water molecules, respectively. For the magnitude of each contribution, see Tables and S5. We note that the jump at the beginning of the free-energy profiles corresponds to the transition of state Iopt,0 to state Iopt in Figure . There is no change in the liquid-phase system but only in the gas-phase system and hence, this is technically not an FEP step.

Table 2

Contributions to Aqueous-Phase Effect on the Free Energy (in eV) of Low Coverage Desorption (ΔΔGAdsorbateliq→gas = –ΔΔGAdsorbategas→liq) of Ethylene Glycol at 423 K (See Figure c), Carbon Monoxide at 298 K (See Figure b), Phenol at 298 K (See Figure d) over the (111) Facets of Pt and Cu Catalysts, and of Benzene at 298 K over the Pt(111) Surface (See Figure S6)a

surface	adsorbate	Iopt,0-to-Iopt	Iopt-to-I (to OPT adsorbate in liquid)	I-to-II (remove electrostatic potential)	II-to-III (electrostatic embedding)	III-to-IV (remove Lennard-Jones potential)	IV-to-IVopt (to OPT site in gas)
Cu(111)	ethylene glycol@423 K	0.04	0.07 ± 0.00	0.49 ± 0.00	–0.05 ± 0.01	–0.64 ± 0.00	–0.05 ± 0.00
	carbon monoxide@298 K	0.04	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	–0.30 ± 0.02	–0.04 ± 0.00
	phenol@298 K	0.02	–0.02 ± 0.00	0.12 ± 0.00	–0.04 ± 0.01	–0.89 ± 0.01	–0.01 ± 0.00
Pt(111)	ethylene glycol@423 K	0.01	–0.01 ± 0.00	0.46 ± 0.00	–0.02 ± 0.01	–0.72 ± 0.01	–0.14 ± 0.00
	carbon monoxide@298 K	0.16	0.00 ± 0.00	0.04 ± 0.00	–0.01 ± 0.00	–0.37 ± 0.02	–0.16 ± 0.01
	phenol@298 K	0.37	–0.01 ± 0.00	0.40 ± 0.00	–0.07 ± 0.02	–1.02 ± 0.01	–0.45 ± 0.00
	benzene@298 K	0.37	0.00 ± 0.00	0.28 ± 0.00	–0.06 ± 0.02	–0.99 ± 0.06	–0.43 ± 0.00

All numbers contain a 95% confidence interval based on limited water sampling obtained by performing multiple independent simulations.

Finally, when considering the vibrational contributions, the overall aqueous solvent effect on the adsorption free energy (ΔΔGAdsorbategas→liq) is 0.14 ± 0.01 eV for the RS and −0.48 ± 0.03 eV for the TS, such that the solvent effect on the activation barrier (ΔΔGliq‡) is −0.62 ± 0.03 eV. This agrees very well with our previous study, where we used eSMS for surface reactions and found ΔΔGliq‡,eSMS = −0.56 ± 0.04 eV.[25]

Phenol Adsorption on Pt(111) at 298 K in Liquid Water

Using the proposed QM/MM-FEP scheme, we computed the solvent effect on the adsorption free energy of phenol at low coverage (ΔΔGPhgas→liq (θPh = 0)) and at 298 K to be 0.78 ± 0.01 eV. To compare our result to the experimental work from Singh et al.,[23,71] we used the adsorption model from Nitta et al.[72,73] that considers the ability of adsorbates to cover multiple active/metal sites and that can be derived from statistical mechanics. According to the experimental work, only hydrogen and phenol are present in the aqueous phase and can adsorb on the Pt surface, hencewhere θPh, θH, and θ* are the surface coverages of phenol, hydrogen, and free site (covered by water), respectively, and n indicates the number of sites that a phenol molecule can occupy when adsorbed on the surface. We note that we have defined the “surface coverage of species i” asi.e., the “real” surface coverage is nθ. Using statistical mechanics and performing some algebra, we reach the following equation for the adsorption isotherm (see the Supporting Information for the full derivation)where Kgas (θPh = 0) is gas-phase adsorption equilibrium constant of phenol at low coverage (at 298 K: ΔHPhgas = −200 kJ/mol,[74] ΔSPhgas = −14.7R[75]), KH,Ph is Henry’s constant for phenol in water (at 298 K: KH,Ph = 5 × 10–4 bar/M[76]), CPh is the concentration of phenol in the aqueous phase, γH is a constant taking into account a hydrogen coverage effect (the γH parameter represents the hydrogen (or any other species) surface coverage relative to the free site coverage, e.g., γH = 0 is the zero H coverage approximation and γH = 10 means the H coverage is 10 times as large as the free site coverage), αPh,Phliq is a constant associated with adsorbed phenol–phenol lateral interactions in liquid, R is the gas constant, T is the system temperature, and finally ΔΔGPhgas→liq(θPh = 0) is the solvent effect on the low coverage phenol adsorption defined asWe note here that ΔGads,Phl and ΔGads,Phg–l are the adsorption free energies of phenol in the aqueous phase from its solvated state in solution and from its free state in the gas phase, respectively, and ΔGPhsolvation is the solvation free energy of phenol, which can be obtained from Henry’s constant at relevant temperatures.[76]

Next, we fitted eq to the experimental data at different phenol concentrations at 298 K to extract ΔΔGPhgas→liq(θPh = 0) and αPh,Phliq. Figure shows the best fits for different n and γH that are experimentally not known in liquid water. A change of γH from 0 to 1 only changes ΔΔGPhgas→liq(θPh = 0) by less than 0.10 eV; however, for a larger γH (high hydrogen coverage) such as 10 (10 times higher hydrogen coverage relative to the free site coverage), this change is about 0.25 eV. Experimental work claims γH < 1.0 since the partial pressure of hydrogen is quite low. Furthermore, the plots in Figure indicate that the fits with lower n predict the higher concentration data better, while a higher n leads to better fits at a low phenol concentration (see also Figure S8). We attribute this observation and the fact that, for a fixed n value, we fit a negative lateral phenol–phenol interaction parameter αPh,Phliq to different phenol orientations on the surface at low and high coverages. At low coverage, phenol adsorbs horizontally parallel to the surface, while increasing the phenol coverage forces the adsorbed phenol molecules to slant upward to occupy less space on the surface for additional adsorbed phenol molecules (see Figure S9). Fortunately, ΔΔGPhgas→liq(θPh = 0) is not very sensitive to the n value and we conclude that independent of these intricacies, the experimental ΔΔGPhgas→liq(θPh = 0) value (1.14 eV) is in good agreement with our computational predictions (0.78 eV) of a strong endergonic solvent effect. Aqueous solvation effects are indeed large for phenol adsorption on Pt(111) with significant consequences for the metal catalyst design for conversion of similar reactants. In this context, we note that using a stronger adsorbing force field for the Pt–water interaction will lead to a reduced difference between the computational prediction and experimental measurement. Interestingly, the computationally predicted low coverage free energy of phenol adsorption in liquid water, ΔGadsg–l ≈ −0.8 eV, is even closer to the experimentally measured one (between −0.5 and −0.6 eV), as PBE-D3 with harmonic approximation slightly underpredicts the gas-phase adsorption free energy relative to the best experimental prediction.

Figure 3

Fits of eq (adsorption isotherm) to the experimental data (green circles) of the aqueous-phase phenol adsorption on a Pt(111) catalyst surface at 298 K from Singh et al.[23,71] (see Figures S10–S12 for similar plots at different temperatures). The γH parameter represents the hydrogen surface coverage relative to the free site coverage (γH = 0 shows zero H coverage approximation and γH = 10 means the H coverage is 10 times as large as the free site coverage), n represents the number of sites a phenol molecule occupies when adsorbed on the surface, α (αPh,Phliq) is a constant associated to the phenol–phenol lateral interactions on the surface in water, and ΔΔG (ΔΔGPhgas→liq (θPh = 0)) is the aqueous-phase effect on the low coverage phenol adsorption on the surface. SSR is the sum of squared residuals, which shows the error associated with each fit.

Benzene Adsorption on Pt(111) at 298 K in Liquid Water

In the gas phase, benzene adsorbs on Pt terrace sites with its molecular plane parallel to the surface, and the most energetically favored adsorption site was found to be the bridge 30 conformation (Figure S7h), which agrees well with previous experimental[77] and theoretical studies.[78] Singh et al.[22] measured aqueous-phase benzene adsorption over a Pt(111) electrode surface[79] and concluded that at ∼50% coverage of benzene (i.e., 0.125 ML or 0.1 ML benzene coverage, assuming each benzene occupies four or five sites, respectively), the free energy of adsorption of benzene in the aqueous phase from its solvated state (C6H6(l) + *(l) ↔ C6H6*(l)) is approximately −14 kJ/mol. To calculate the adsorption free energy of benzene in the aqueous phase from its free state in the gas phase (C6H6(g) + *(l) ↔ C6H6*(l)), we can utilize Henry’s law constant for benzene at 298 K (KH = 5 bar M–1),[22] ΔGsolvation = −RT ln KH = +4 kJ/mol. Hencewhere ΔGadg–l is the free energy of adsorption of benzene in the aqueous phase from its gas-phase reference state. To compute the solvent effect on the benzene adsorption free energy, we also need to estimate the free energy of benzene adsorption in the gas phase (C6H6(g) + *(g) ↔ C6H6*(g)) such thatwhere ΔΔGgas→liq is the solvent effect on benzene adsorption. Campbell et al.[77] studied the gas-phase heat of adsorption of benzene over the Pt(111) surface at 300 K using single crystal adsorption calorimetry (SCAC) and reported that the heat of adsorption of benzene follows a quadratic expression in coverage up to 0.14 ML coverage, ΔHadsg = 197 – 314θ – 3546θ2 kJ/mol. Hence, at 0.125 ML coverage, the heat of adsorption of benzene in the gas phase is computed to be −102.34 kJ/mol. At 0.1 ML coverage, the gas-phase heat of adsorption becomes −130.14 kJ/mol. The entropy loss of an adsorbed molecule can also be calculated utilizing Campbell’s expressions,[80] ΔSads = −0.3Sgas – 3.3R = −0.11 kJ/(mol·K), where the standard entropy of gas-phase benzene at 298 K is Sgas = −269.01 J/(mol·K).[81] This leads to a free energy of benzene adsorption (assuming 0.125 ML coverage) in the gas phase of ΔGadsg = ΔHads – TΔSads = −70.12 kJ/mol. At 0.1 ML coverage, the adsorption free energy is computed to be −97.91 kJ/mol. Finally, we can employ eq to compute the solvent effect on the adsorption free energy of benzene at 0.125 ML coverageAt 0.1 ML coverage, this becomesUsing our QM/MM-FEP scheme, we computed the solvent effect on the adsorption free energy of benzene at 298 K to be 0.82 ± 0.06 eV, which is again in very good agreement with the aforementioned experimental prediction. We note that we used the TIP4P water model[82] for benzene calculations and the TIP3P water model[83] for phenol calculations. Given that both water models generally predict very similar results for biological systems, it is instructive to observe that also in our simulations, good agreements with the experiment are observed at room temperature independent of the specific water model (that admittedly use an identical water geometry).

Aqueous-Phase Effect across Adsorbates and Metals

To elucidate the significance of these large, computed and measured, solvation effects for the Pt–phenol system for various other catalytic processes, we computed the aqueous solvation effect on the adsorption free energies (ΔΔGAdsorbategas→liq = ΔGAdsorbateg–l – ΔGAdsorbategas) of four different adsorbates, ethylene glycol (EG), carbon monoxide (CO), benzene, and phenol (Ph), over the (111) surface facet of Cu and Pt. In addition to our QM/MM-FEP calculations, we also used two different well-known implicit solvation methods, iSMS[84] and VASPsol[85,86] (see Table ). Positive solvent effects on the adsorption free energies indicate that liquid water destabilizes the adsorbate. In this regard, the calculations at 298 K indicate that an aqueous phase destabilizes all adsorbed species and that the impact is most pronounced for phenol and benzene and smallest for CO, irrespective of the metal surface (ΔΔGPhgas→liq ≈ ΔΔGBenzenegas→liq > ΔΔGCOgas→liq). We attribute the destabilization of the adsorbates on the surface primarily to the solvent cavity formation energy and the desire of water molecules to adsorb on a metal surface.

Table 1

Aqueous-Phase Effect on the Free Energy of Low Coverage Adsorption (ΔΔGAdsorbategas→liq = ΔGAdsorbateg–l – ΔGAdsorbategas) in eV for Ethylene Glycol over Cu(111) and Pt(111) at 423 K, for CO on Cu(111) and Pt(111) at 298 K, for Phenol over Cu(111) and Pt(111) at 298 K, and for Benzene over Pt(111) at 298 Ka

			ΔΔGAdsorbategas→liq
adsorbate/surface/T (K)	ΔGAdsorbategas	ΔGAdsorbateg–l,QM/MM-FEP	QM/MM-FEP	iSMS	VASPsol
ethylene glycol/Cu(111)/423	0.45	0.59 ± 0.01	0.14 ± 0.01	0.02	–0.15
ethylene glycol/Pt(111)/423	0.35	0.80 ± 0.02	0.45 ± 0.02	–0.19	–0.21
carbon monoxide/Cu(111)/298	–0.22	0.07 ± 0.01	0.29 ± 0.01	–0.03	–0.02
carbon monoxide/Pt(111)/298	–1.24	–0.89 ± 0.03	0.35 ± 0.03	0.01	–0.05
phenol/Cu(111)/298	–0.59	0.26 ± 0.01	0.85 ± 0.01	–0.19	–0.22
phenol/Pt(111)/298	–1.61	–0.83 ± 0.01	0.78 ± 0.02	–0.41	–0.47
benzene/Pt(111)/298	–1.55	–0.73 ± 0.06	0.82 ± 0.06	–0.42	–0.41

The aqueous phase results are based on the proposed scheme in this work (QM/MM-FEP) and two implicit solvation schemes: iSMS and VASPsol. For the QM/MM-FEP calculations, the 95% confidence interval (based on limited water sampling and multiple independent simulations) is also given.

Furthermore, the metal identity also plays a key role in the solvation effect, primarily through an indirect influence of the solvent environment on the electronic structure of the metal atoms that in turn affects the stability of the adsorbed species.[25] To elucidate the root cause for the different solvent effects across the metals and the adsorbates, we developed a linear model for describing the variability of the solvation free energy of adsorption (s is the standard deviation and the bar sign shows mean value) in which f1 and f2 are the number of hydrogen bonds (H bond) between the adsorbate and water molecules (we employed a geometric definition in which a hydrogen bond exists if the distance between the donor oxygen (Od) and the acceptor oxygen (Oa), ROO, is less than 3.2 Å and the angle ∠HOdOa is smaller than 20°; see Figure S13) and number of surface sites occupied by the adsorbate (#Site), respectively (see Table S6). Based on the best fit parameters (α1 = −0.44, α2 = 0.86), we conclude that different solvent effects arise primarily from the size of the molecule and partially from the different H-bonding between the adsorbate and water molecules. We note that we also tested the net charge on the adsorbate, number of waters displaced by the adsorbate, and gas-phase adsorption energy as possible descriptors, but we found the best fit using the H-bond and #Site. The number of displaced water molecules has previously been proposed as a meaningful descriptor,[22] and it is related to the size of the molecule or the #Site descriptor. We show in Table S7 the number of displaced waters by adsorption of CO, phenol, and benzene over Pt(111) at 298 K, 1.88 ± 0.71, 5.46 ± 0.97, and 4.43 ± 0.49, respectively, which indicates a correlation of the number of displaced water molecules and the destabilization in adsorption free energies in solution. Possibly, the correlation with the number of displaced water molecules is worse than the #Site descriptor because of the larger uncertainties in our predicted number of displaced water molecules.

Analysis of the Specific Contributions to the Solvation Free Energy of Adsorption

Figure illustrates the free-energy profiles of the aqueous solvent effect on the desorption of the various adsorbates. The profiles are comprised of five subsections: Iopt,0-to-I, I-to-II (remove electrostatic potential), II-to-III (electrostatic embedding), III-to-IV (remove Lennard-Jones potential), and IV–IVopt (To OPT site in gas), whose individual values are listed in Tables and S5.

Iopt,0-to-I

This includes steps Iopt,0-to-Iopt and Iopt-to-I. The former depends on how much the adsorbate distorts the coordinates of the surface metal atoms upon gas-phase adsorption. It means this contribution is larger for adsorbates that bind to the surface more strongly. The contribution of step Iopt-to-I is often very small since the gas and liquid optimized coordinates are often similar. For EG, CO, benzene, and phenol on Pt(111), the solvent effect of this contribution is essentially zero (within the error bars) for all species.

Classical Electrostatic Contribution to the Solvent Effect (I-to-II)

Removal of the classical electrostatic interactions between the adsorbate and the liquid water molecules generally increases the total free energy of the system, i.e., the classical electrostatic interactions stabilize all adsorbed species. As expected, this effect is most dominant for ethylene glycol, somewhat smaller for phenol, and least relevant for CO, which is attributed to the fact that attractive forces between charged species decrease the energy of the system, while repulsive forces between charged species have the opposite effect. In this context, perturbing charges of the adsorbate atoms from their original magnitude to zero decreases the attractive forces between the adsorbate and the water molecules and in turn increases the energy of the total system.

Nonclassical Electrostatic Effects due to Changes in Polarization by the Solvent (II-to-III)

As described above, removal of the nonclassical electrostatic contribution of the water molecules changing the electronic structure of the quantum system usually decreases the free energy of the system in the desorption process shown in Figures and 2, increasing the adsorbate stabilization from the classical electrostatic contributions.

Lennard-Jones Contribution to the Solvent Effect (III-to-IV)

Removing the LJ interactions between the adsorbed species and the water molecules generally decreases the free energy of the system as shown in Table and Figure . Slow removal of the LJ interaction energy leads to a reduction of the water cavity for the adsorbates. Cavity formation is generally an endergonic process and given that any adsorbate requires a similar-sized water cavity on any metal surface, we found this LJ removal contribution to be almost independent of the metal surface (see Table ). This again explains why larger adsorbates that require larger water cavities are often also more destabilized on a metal surface.

Transition of Metal–Cluster Coordinates (IV–IVopt) Contribution to the Solvent Effect

This contribution is significant for adsorption of molecules that bind strongly to the surface in the gas phase and thus distort the coordinates of the surface metal atoms upon adsorption. For EG, CO, benzene, and phenol on Pt(111), the zero-point corrected adsorption energies, −0.35, −1.73, −2.23, and −2.29 eV, respectively, correlate with this contribution to the overall solvent effect on adsorption of 0.14, 0.16, 0.43, and 0.45 eV, respectively. We note that although the magnitude of this contribution is similar to that of Iopt,0-to-Iopt, both their physical meaning and computations involved are different. This transition occurs in the presence of water and requires FEP steps, while the Iopt,0-to-Iopt contribution is only the difference between the self-consistent field (SCF) energy of the gas-phase optimized free site (metal slab) and the single-point SCF energy of the metal slab, which is obtained by removing the adsorbate from the optimized metal–adsorbate slab.

Implicit (iSMS and VASPsol) vs Explicit (QM/MM-FEP) Solvation Models

Finally, we compare the results of the solvent effect on the adsorption free energies of EG, CO, benzene, and phenol over Pt(111) and Cu(111) computed with our QM/MM-FEP scheme against two different implicit solvation schemes: VASPsol and our iSMS methodology. Table shows that both implicit solvation methods fail to capture the true solvent effect, probably because they underestimate the endergonic cavity formation energy. In contrast, our explicit solvation scheme agrees quite well with experimental data probably because the metal–water force field is specifically optimized for this free-energy contribution.[87] Interestingly, both implicit solvation schemes predict similar solvation free energies at 298 K, the temperature for which the VASPsol parameters have been optimized. While iSMS permits a straightforward change in the system temperature, we could only adjust the dielectric constant for water with temperature in VASPsol and otherwise used the default VASPsol parameters. Thus, at 423 K, the solvation free-energy predictions deviate significantly between the two implicit solvation methods, and iSMS predicts an endergonic solvent effect of 0.02 eV for EG adsorption on Cu(111) at 423 K, which agrees at least qualitatively with the value of 0.14 ± 0.01 eV predicted using our QM/MM scheme. To improve implicit solvation models for metallic surface systems, we recommend focusing on the cavity formation energy description.

Enthalpic and Entropic Contributions: Van’t Hoff Plots

To compute the enthalpic and entropic contributions to the solvent effect on the adsorption free energies, we repeated our explicit solvation calculations (QM/MM-FEP) for CO and phenol on Pt(111) at five different temperatures within the temperature interval of 285–353 K and constructed van’t Hoff plots shown in Figure . The results indicate that the enthalpic and entropic contributions to the solvation free energy of adsorption are comparable for phenol and that phenol adsorption on Pt(111) in liquid water is destabilized both enthalpically and entropically (ΔΔHPhgas→liq = 0.45 ± 0.10 eV, −TΔΔSPhgas→liq@298 K = 0.33 ± 0.08 eV). For CO adsorption, our error bars for the free energy of solvation are too large to unambiguously determine the enthalpic and entropic contributions with a van’t Hoff plot (ΔΔHCOgas→liq = 0.13 ± 0.15 eV, −TΔΔSCOgas→liq@298 K = 0.22 ± 0.15 eV). Still, it appears that the entropic solvation effect is not negligible relative to the enthalpic contribution. Heenen et al.[88] also investigated the aqueous-phase CO adsorption over Pt(111) using AIMD simulations and reported ΔΔHCOgas→liq = −0.15 ± 0.12 eV. In another study, Bodenschatz et al.[89] investigated the aqueous-phase CO adsorption over Pt(111) surface using DFT and MD simulations and reported ΔΔHCOgas→liq = −0.01 ± 0.09 eV. Although the enthalpic contribution in CO adsorption in these two studies is negative, while ours is positive, we emphasize that the error bars in all of these simulations are too large to conclusively determine the sign of the solvation effect on the enthalpy. Here, we add that AIMD-based methods have significant uncertainties originating from the limited configuration space sampling that are difficult to estimate.

Figure 4

Temperature dependence of the aqueous-phase effect on the free energy of low coverage adsorption of phenol and carbon monoxide over a Pt(111) surface (ΔΔG = ΔΔGAdsorbategas→liq(θAdsorbate = 0)). The error bars indicate 95% confidence intervals based on limited water sampling and multiple (three) independent simulations.

Phenol on Pt(111): Experimental vs Computational

Singh et al.[22,24,71] recently also studied adsorption of aqueous phenol at four different temperatures (283, 288, 298, 314 K) to acquire experimental enthalpic and entropic solvation contributions. We also fitted eq to their measured data (see Figures and S10–S12) and Table shows the extracted ΔΔGPhgas→liq(θPh = 0) values at different temperatures. Overall, the agreement between the computed (QM/MM-FEP) and experimental solvation free energies of adsorption is good for the low coverage adsorption free energy of phenol on Pt(111) at all temperatures. However, the experimental solvation enthalpies and entropies of phenol adsorption in water disagree with our computational data. Large enthalpic solvation effects are predicted that are reduced by entropic contributions (ΔΔHPhgas→liq = 1.83 eV, ΔΔSPhgas→liq = 2.30 meV/K). In our opinion, the solvation effect on the enthalpy and entropy of adsorption is too large and an artifact of the small experimental temperature interval in the van’t Hoff plot together with sizable errors of ∼0.1 eV of the experimental solvation effects on the Gibbs free energies of adsorption. Specifically, in the small temperature interval from 283 to 314 K, the experimentally extracted ΔΔGPhgas→liq value differs by only 0.08 eV (see Table ), which is smaller than the error bar of ∼0.1 eV, making it impractical to even predict the sign of the solvation effect on the entropy of adsorption. A consequence of this small temperature window in the van’t Hoff plot is that the apparent experimentally extracted enthalpy of adsorption from the gas phase to an adsorbed phenol in liquid water (ΔHads,Phg–l) is barely exothermic, while the entropy of adsorption (ΔSads,Phg–l) appears to increase. In contrast, when using either DFT computed or experimentally measured gas-phase enthalpies and entropies of adsorption together with our QM/MM predicted enthalpies and entropies of solvation, we predict sizable phenol enthalpies of adsorption of −1.62 to −1.86 eV and (as expected) negative entropies of adsorption of −2.37 to −3.42 meV/K.

Table 3

Temperature Dependence of the Aqueous-Phase Effect on Free Energy (ΔΔGPhgas→liq), Enthalpy (ΔΔHPhgas→liq), and Entropy (ΔΔSPhgas→liq) of Low Coverage Phenol Adsorption on the Pt(111) Catalyst Surfacea

	T (K)	ΔΔGPhgas→liq (eV)	ΔΔHPhgas→liq (eV)	ΔΔSPhgas→liq (meV/K)	ΔHads,Phg–l (eV)	ΔSads,Phg–l (meV/K)
experimental	283	1.16 (1.44)	1.83 (1.82)	2.30 (1.28)	–0.25 (−0.26)	1.03 (0.01)
	288	1.15 (1.43)
	298	1.14 (1.44)
	314	1.08 (1.39)
QM/MM-FEP	285	0.77 ± 0.01	0.45 ± 0.10	–1.10 ± 0.28	–1.86 ± 0.10 [−1.62 ± 0.10]	–3.42 ± 0.28 [−2.37 ± 0.28]
	298	0.78 ± 0.02
	315	0.80 ± 0.01
	335	0.82 ± 0.02
	353	0.85 ± 0.01

Experimental energy and entropy data obtained by fitting eq (n = 5 and γH = 10; the numbers in parentheses in the experimental row obtained using n = 5 and γH = 0.1 in fitting) to the raw adsorption data from Singh et al.[22,23,71] The γH parameter represents the hydrogen surface coverage relative to the free site coverage (γH = 0 shows zero H coverage approximation and γH = 10 means the H coverage is 10 times as large as the free site coverage), n represents the number of sites a phenol molecule occupies when adsorbed on the surface. Computational data obtained by the QM/MM-FEP scheme proposed in this study. ΔHads,Phg–l and ΔSads,Phg–l are the enthalpy and entropy of phenol adsorption in liquid water from its free state in the gas phase (Ph(g)+ n*(liq) ↔ Ph* (liq)), respectively (ΔΔZPhgas→liq = ΔZPhg–l – ΔZPhgas, Z ≡ G, H, S). Numbers in bracket [ ] in the QM/MM-FEP row were calculated using the experimental gas-phase enthalpy and entropy of phenol adsorption on the Pt(111) catalyst surface (ΔHPhgas = −200 kJ/mol,[74] ΔSPhgas = −14.7R[75]).

Conclusions

In this study, we have developed an explicit solvation scheme, based on hybrid QM/MM-FEP calculations, for quantifying solvation effects on adsorption free energies for heterogeneous catalysis applications. To numerically verify the proposed scheme, the solvent effect on the activation free energy of O–H bond cleavage in EG over Cu(111) has been computed through two adsorption free energy of solvation calculations (reactant and transition states) and by performing conventional eSMS calculations for surface reactions. Using the adsorption calculations, we predict a solvation effect on the activation free energy of −0.62 ± 0.03 eV, which agrees very well with our conventional eSMS calculations (ΔΔGliq‡ = −0.56 ± 0.04 eV) and which highlights that the sampling error is similarly small in the new scheme that has the added benefit of being able to describe adsorption events.

Next, we computed solvation free energies of phenol and benzene adsorption on Pt(111) and reanalyzed experimental data for solvation free energies of phenol and benzene adsorption on this surface using multiple adsorption isotherm models. Our explicit solvation calculations predict an endergonic solvation free energy of 0.78 ± 0.02 eV for phenol and 0.82 ± 0.06 eV for benzene, which are in very good agreement with measurements to within the experimental and force field uncertainties. Implicit solvation schemes such as VASPsol and iSMS fail to capture the true solvent effect, probably because they underestimate the endergonic cavity formation energy. To elucidate the significance of these large, computed and measured, solvation effects for the Pt–phenol and the Pt–benzene system for various other catalytic processes, we computed the aqueous solvation effect on the adsorption free energies of carbon monoxide, ethylene glycol, and phenol over the Pt(111) and Cu(111) surfaces. Irrespective of the metal surface, the water environment destabilized all three adsorbed species with the largest impact on phenol and benzene and smallest for CO (ΔΔGPhgas→liq ≈ ΔΔGBenzenegas→liq > ΔΔGCOgas→liq), which is attributed primarily to the size of the molecule (number of sites occupied) and the corresponding solvent cavity formation and partially to different H-bonding between the adsorbate and water molecules.

Methods

To describe the potential energy surface of a metal–adsorbate system during desorption with a hybrid QM/MM scheme, adsorbate, active site, and its immediate metal neighbors are treated from first principles (called QM system/subsystem in this study), while the bulk of the solvent molecules and the “nonreactive” part of the simulation system (metal atoms far away from the adsorbate) are described by classical force fields (called MM system/subsystem in this study). Here, we explain why we need to use different software in our QM/MM scheme. First, the interaction between metal surfaces and adsorbates in vacuum can be described well by periodic planewave DFT calculations (e.g., VASP[60,61]). Electrons flow freely in metals, which leads to long-range interactions that can typically not be described by cluster model DFT calculations. In the presence of liquid water molecules, the metal–adsorbate system gets polarized by the solvent. Fortunately, this polarization is short-ranged and thus, it can be described by periodic electrostatic embedded cluster method (PEECM[65]) calculations,[20] as implemented in the TURBOMOLE[62−64] program package, of a small metal-cluster–adsorbate system.[62−64] To avoid overcounting of interactions, another nonperiodic DFT cluster calculation has to be performed in the absence of the (periodic) solvent with the TURBOMOLE program package. In principle, the polarization of QM atoms by the water point charges can also be computed with periodic DFT codes (although not with VASP); however, due to the long-range water–water interactions, very large metal slabs containing thousands of metal atoms would be necessary, which is very slow and therefore not done here. Next, to explain the classical solvent interactions (electrostatic and van der Waals) within the solvent and with the QM subsystem, classical molecular dynamics simulations must be performed (e.g., DL_POLY[90]). The following sections detail the calculation methods employed in each of the aforementioned program packages.

Planewave DFT Calculations

Gas-phase DFT calculations were carried out employing periodic boundary conditions and using the Vienna Ab Initio Simulation Package (VASP 5.4.4).[60,61] The frozen-core, all-electron projector augmented-wave (PAW)[91] method was utilized to avoid the singularities of Kohn–Sham wavefunctions at the nuclear positions. The exchange-correlation energy is calculated within the generalized gradient approximation (GGA)[92] using the Perdew–Burke–Ernzerhof (PBE)[93,94] functional. For phenol adsorption, we used Grimme’s DFT-D3 methodology[95] as implemented in VASP to describe the van der Waals interactions. Brillouin zone integrations have been performed with a 4 × 4 × 1 Monkhorst–Pack[96]k-point grid, and electronic wavefunctions at each k-point were expanded using a discrete planewave basis set with kinetic energies limited to 400 eV. For phenol adsorption, we increased the energy cutoff to 450 eV although for all practical purposes, energies are already converged with a cutoff of 400 eV. Fractional occupancies of bands were allowed within a window of 0.10 eV using a first-order Methfessel–Paxton[97] smearing method, which permitted us to calculate the entropic contributions due to accurate smearing. The modified version of the Makov–Payne[98] method with Harris corrections was applied to the stress tensor and forces to calculate dipole and quadrupole corrections (along the surface normal) to the total energy. Using the supercell approach, a 4 × 4 unit cell with four layers of metal atoms, in which the bottom two layers were fixed, while the top two layers were relaxed, has been constructed employing a 15 Å vacuum on top of the surface, which restricted the interaction between the periodic images along the surface normal. Self-consistent field (SCF) calculations were converged to 1.0 × 10–7 eV. A force criterion of 0.02 eV Å–1 was used on relaxed atoms for geometry optimizations. The transition-state structure of the O–H bond cleavage of ethylene glycol (EG) over a Cu(111) surface model was located using a combination of climbing-image nudged elastic band[99,100] and dimer[101,102] methods. Finally, the minima and the first-order saddle points were validated by computing the Hessian matrix and vibrational spectra.

Periodic Electrostatic Embedded Cluster and Nonperiodic Cluster Calculations

Cluster model DFT calculations in a water point charge field and in vacuum have been carried out using the TURBOMOLE 7.2 program package.[62−64] We used a cluster model consisting of two layers with a total number of 51 metal atoms. Convergence of the total QM/MM system energy with respect to the lateral size and depth of the cluster geometry can be found elsewhere.[59] To represent the adsorbate atoms and the metal atoms, an improved version of the default TURBOMOLE basis sets (def-bases) with split valence and polarization functions (def2-SVP)[103,104] was used. Pt atoms were represented using scalar relativistic effective core potentials (ECPs) in conjunction with split valence basis sets augmented by polarization functions.[104,105] Exchange-correlation effects were calculated using the PBE functional[93,94] (for phenol adsorption, Grimme’s DFT-D3 methodology[95] was used). Finally, the RI-J approximation with auxiliary basis sets was used to approximate the Coulomb integrals.[106,107] Single-point energy calculations were performed with an SCF energy convergence criterion of 1.0 × 10–7 Hartree, and a Gauss–Chebyshev type spherical grid, m4, was employed to perform the numerical integrations.[63]

Molecular Dynamics (MD) Simulations

The DL_POLY 4.03 molecular simulation program package[90] was utilized for the MD simulations. Using the supercell approach, the initial 4 × 4 unit cell for each clean surface metal was augmented laterally to a 16 × 20 surface with further vacuum added in the Z-direction resulting in a simulation box containing 1280 metal atoms. As a result, a simulation box size of 44.98 Å × 48.69 Å × 49.01 Å for a Pt(111) surface and 41.04 Å × 44.43 Å × 49.01 Å for a Cu(111) surface were achieved. The QM cluster atoms were then substituted into the constructed surface. The simulation box height was selected based on the work from Behler et al.[108] who found that simulations of metal–water interfaces should contain a water layer of ∼40 Å height. Subsequently, the experimental saturated liquid water density of ∼1.0 g/cm3 at 298 K for CO, benzene, and phenol adsorption calculations on Pt(111) and Cu(111) surfaces and of ∼0.92 g/cm3 at 423 K for O–H bond cleavage’s calculations of ethylene glycol (EG) on Cu(111) as well as EG adsorption on Pt(111) were achieved by packing the simulation box of the Pt(111) and Cu(111) surfaces with ∼3000 and ∼2400 water molecules, respectively. All metal and adsorbate atoms were kept fixed, while the geometry of water molecules was constricted to that of TIP3P[83] geometry with the RATTLE algorithm,[109] a velocity version of the SHAKE algorithm,[110] in conjunction with the velocity Verlet (VV) integrator[111] to solve Newton’s equations of motion. The force field parameters of liquid water were obtained from the TIP3P model (we used the TIP4P/2005 water model[82] for the benzene calculations as we are currently in the process of studying solvation effects on benzene adsorption at elevated temperatures where the TIP4P/2005 water model is better at reproducing the water phase behavior—at room temperature we do not expect significant deviations in simulation results between water models as the metal–water interaction potential has been developed independent of the specific water model and both water models are commonly used for the simulation of room-temperature phenomena). The Lennard-Jones parameters were obtained from the OPLS force field[112−114] for EG, benzene, and phenol and from Straub et al.[115] for CO. In addition to the OPLS parameters, the Lennard-Jones parameters from the combined B3LYP/6-31_G*/AMBER Potential[116] were used for the hydrogen atoms of EG. Lennard-Jones parameters for hydrogen atoms are important in QM/MM optimizations that permit hydrogen atoms to approach water molecules and leave the protective environment of a neighboring carbon or oxygen atom. The metal–water interaction was represented by a Lennard-Jones metal potential.[87] The LJ cross-term of the intermolecular parameters are calculated by Lorentz–Berthelot mixing rules through equations and . All Lennard-Jones parameters are included in Table S1. For O–H bond breaking of EG on the Cu(111) surface, and EG adsorption on the Pt(111) surface, the charges for the QM atoms were estimated using the natural population analysis (NPA),[66] while for CO, benzene, and phenol adsorptions on the Pt(111) and Cu(111) surfaces, the DDEC6[67] (a refinement of the density-derived electrostatic and chemical (DDEC) approach) charge model was used (to compare the performance of NPA and DDEC6 in our scheme, we also repeated the calculations of CO adsorption on Pt(111) with the NPA charge model. See page S3 of the Supporting Information for a discussion on the selection of the charge model).

To describe the interaction of the TIP3P (TIP4P/2005 for the benzene calculations) water point charges with the quantum chemically described cluster model, we employed the periodic electrostatic embedded cluster method (PEECM)[65] as implemented in TURBOMOLE. Simulations were carried out in a canonical ensemble (NVT) with the Nosé–Hoover thermostat[117,118] with a 1 ps relaxation time constant. We selected the NVT ensemble over the NPT ensemble since the system is asymmetric and the dimensions in two directions are constant and controlled by the metal slab. We note that all ensembles lead to the same results in the limit of large systems, and we previously showed that our simulation system is sufficiently large for obtaining converged results.[44] Electrostatic interactions were accounted for using the Smoothed Particle Mesh Ewald (SPME) method[119] with automatic parameter optimization with default SPME precision. A 12 Å cutoff radius was adopted for the van der Waals interactions and the transition between short- and long-range electrostatic interactions. If not specified differently, all systems were equilibrated for 250 ps and sampled for 1000 ps (1 ns) using a 1 fs timestep to obtain 1000 MM conformations, which are 1 ps apart, for each experiment. To optimize structures in an aqueous-phase environment, we utilized the fixed-size ensemble approximation with 5000 MM conformations (5 ns sampling) recorded every 1 ps.