Aromatic Molecule Solvation in Liquid Water with Coupled Cluster Accuracy: The Balance of Pi-Interactions and Hydrophobicity

Abstract

Aromatic organic solutes in water exhibit a delicate balance between hydrophobic solvation and directional O-H$\cdots π$ hydrogen bonds, yet widely used force fields and state-of-the-art density functional approaches struggle to provide a consistent picture of these pivotal interactions. We introduce a data-efficient upfitting strategy to train a machine learning interatomic potential (MLIP) based on the graph atomic cluster expansion for aqueous aromatic molecules with CCSD(T) accuracy for condensed phase simulations, using only finite molecular clusters. We apply our method to aqueous toluene (C$_6$H$_5$CH$_3$). The resulting CCSD(T)-quality MLIP reproduces coupled cluster energies and forces in bulk and reveals that commonly employed methods do not capture the crucial balance between hydrophilic and hydrophobic solvation, distorting the interactions of aromatic molecules with their environment. Representative biomolecular force fields substantially understructure the hydrophobic solvation shell and misorient interfacial water, while overestimating $π$-contacts, yielding an inconsistent solvation balance. Even hybrid DFT and MP2 overestimate barriers to breaking of water-$π$ hydrogen bonds. Our workflow provides a practical, general route to CCSD(T)-quality condensed-phase simulations of aqueous solutions, and thus constructed interaction potentials now open the door to consistent, highly accurate benchmark studies of $π$-contacts and hydrophobic effects in biomolecular contexts such as solvation of proteins and DNA in aqueous environments.