  The effect of solvent on the free energy of reaction intermediates adsorbed
on electrocatalyst surfaces can significantly change the thermochemical
overpotential, but accurate calculations of this are challenging. Here, we
present computational estimates of the solvation energy for reaction
intermediates in oxygen reduction reaction (ORR) on a B-doped graphene (BG)
model system where the overpotential is found to reduce by up to 0.6 V due to
solvation. BG is experimentally reported to be an active ORR catalyst but
recent computational estimates using state-of-the-art hybrid density
functionals in the absence of solvation effects have indicated low activity. To
test whether the inclusion of explicit solvation can bring the calculated
activity estimates closer to the experimental reports, up to 4 layers of water
molecules are included in the simulations reported here. The calculations are
based on classical molecular dynamics and local minimization of energy using
atomic forces evaluated from electron density functional theory. Data sets are
obtained from regular and coarse-grained dynamics, as well as local
minimization of structures resampled from dynamics simulations. The results
differ greatly depending on the method used and the solvation energy estimates
and are deemed untrustworthy. It is concluded that a significantly larger
number of water molecules is required to obtain converged results for the
solvation energy. As the present system includes up to 139 atoms, it already
strains the limits of computational feasibility, so this points to the need for
a hybrid simulation approach where efficient simulations of much larger number
of solvent molecules is carried out using a lower level of theory while
retaining the higher level of theory for the reacting molecules as well as
their near neighbors and the catalyst.
