On the challenge of obtaining an accurate solvation energy estimate in simulations of electrocatalysis

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.