Local environment neighbor sensitivity analysis: visualization of cation effect at liquid–solid interface

Abstract

We present a method to extract the effects of environmental neighbors on local energies in heterogeneous systems from machine-learning force fields and demonstrate its ability to visualize cation-induced stabilization of hydroxyl groups on Pt in water.

---

In recent years, machine-learning force fields (MLFFs) have attracted significant attention as promising models for describing interatomic interactions, offering an alternative to conventional first-principles (FP) calculations^1,2 and classical force fields.^3,4 MLFFs are supervised machine-learning models trained on FP data for potential energies, atomic forces, and stress tensors to predict interatomic interactions.^5–15 Unlike classical force fields, which rely on fixed functional forms based on physical laws such as Coulomb and dispersion interactions, MLFFs employ flexible functions developed in the machine-learning (ML) community to represent interatomic interactions. As a result, MLFFs can describe a wide range of interactions—including ionic, dispersive, metallic, covalent, and hydrogen bonding—within a single framework, achieving near FP accuracy while being several orders of magnitude faster than FP calculations. Owing to these properties, MLFFs enable quantitative predictions for complex multi-element organic and inorganic materials that are difficult to treat accurately with classical force fields and computationally prohibitive for FP methods. Their applications have therefore expanded to phenomena such as phase transitions, thermal and ionic transports, chemical reactions and catalytic processes in gases, liquids, solids, and their interfaces.^5–12,14,15

The development of MLFFs has progressed rapidly, and a wide variety of functional forms has been proposed. Nevertheless, most MLFFs represent interatomic interactions based on two common concepts: descriptors of many-body interactions and nonlinear mapping. Because the potential energy E of a material originates from many-body interactions, it can be expressed using the following cluster expansion:^16–19

Here, R_i denotes the position vector of the i-th atom, and R^N represents the set of position vectors of all N atoms in the system. E^(ν) denotes a ν-body interaction, which is defined as a term whose (ν + 1)-th derivative with respect to atomic positions vanishes, as ∂^ν+1E^(ν)/∂R_i1⋯∂R_iν+1 = 0.²⁰ In MLFFs, descriptors representing ν-body local structures are introduced, and the potential energy E is expressed as a function of these descriptors. As the order ν of the descriptors is systematically increased, E can be represented more accurately. However, the number of descriptors grows exponentially with ν, leading to a corresponding increase in computational cost. To address this issue, many MLFFs restrict the descriptor order to ν ≈ 3 or 4 and instead represent many-body effects through higher-order terms generated by nonlinear functions of these descriptors. This strategy for efficiently incorporating high-order many-body interactions is known as the “density trick”.¹¹ Neural networks and kernel functions are commonly used as the nonlinear functions.

In this way, MLFFs can predict interatomic interactions with high accuracy and efficiency. However, a major challenge for MLFFs lies in the difficulty of physically interpreting their predictions. In classical force fields, where interactions are predefined and decomposed into contributions such as Coulomb and dispersion interactions, predicted energies can be analyzed in terms of their underlying interaction components. Similarly, FP calculations provide deep physical insight through quantities such as electron densities, wave functions, and Hamiltonian matrices. In contrast, such analyses are not directly available in MLFFs.

Here, we present a method within the MLFF framework to visualize the surrounding species and their positions that influence the energy of a specific site in a system, and apply it to interfacial electrochemical reactions. As a case study, we focus on the effects of alkali cations on OH adsorbed on a Pt(111) surface in aqueous solution. Experiments have shown that small alkali cations, such as Li⁺, stabilize OH and slow the oxygen reduction reaction (ORR), whose elementary steps include OH reduction.²² This was the first observation demonstrating cation effects in interfacial electrocatalytic reactions and subsequently highlighted the importance of cation effects in reactions such as hydrogen evolution^23,24 and CO₂ reduction.^24–26 However, assessing the stability of hydrated cations fluctuating in water at finite temperature remains challenging. This difficulty arises because interfacial solvent molecules reorganize on nanosecond time scales,²⁷ which are difficult to access using conventional FP methods^22,28,29 that rely on structural optimization at absolute zero or MD simulations limited to picosecond time scales. As a result, many aspects of cation effects on surface adsorbates remain unclear. In this work, we compute the effects of Li⁺, K⁺, and Cs⁺ on the redox potential of OH reduction using a recently developed MLFF-aided FP thermodynamic integration (TI) approach based on finite-temperature MD simulations. We show that smaller cations stabilize OH and lower the redox potential in silico. We then demonstrate, through conventional analyses of cation distributions and interfacial diffusion coefficients, that smaller cations are trapped near OH and exhibit reduced diffusivity. Finally, by applying the visualization method, we obtain further insight into the relationship between the localization of cation distributions and the stabilization of OH, suggesting that attractive interactions between alkali cations and OH play an important role.

The visualization method is constructed on the basis of a kernel-based MLFF approach^30–32 implemented in the Vienna Ab initio Simulation Package (VASP).^33,34 Similar to the Gaussian Approximation Potential (GAP) proposed by Bartók and co-workers,^17,35 this method expresses the total potential energy as a sum of atomic energies:

The atomic energy E_i is represented as a functional of the atom density distribution ρ_i surrounding atom i as E_i = F[ρ_i(r)]. In the present study, the density ρ_i is represented by the three-body angular distribution function ρ_i⁽³⁾(r,s,θ),^17,32 which gives the probability of finding one atom at a distance between r and r + dr from atom i, and another atom at a distance between s and s + ds from atom i, such that the angle ∠kij lies between θ and θ + dθ. This probability is expressed as ρ_i⁽³⁾(r,s,θ)r²s²sinθdrdsdθ. Accordingly, r denotes the rotationally invariant coordinate (r,s,θ) [see details in Section S1 of SI]. In practice, ρ_i is expanded in an orthonormal basis set {ϕ_j|j= 1,…,N_D} as

and the vector of the resulting expansion coefficients, x_i = (c_1i,…,c_NDi)^T, is used as a descriptor.³² The functional F is represented as a linear combination of kernel basis functions that measure the similarity between the descriptor x_i of atom i and those of reference atoms i_B = 1,…,N_B:

The reference atoms are selected from structures providing the training data, and the regression coefficients w_iB are determined to optimally reproduce the training data of energies, forces and stress tensor components. For the kernel basis functions, we use the smooth overlap of atomic positions (SOAP).¹⁷

In this MLFF, the contribution of each atom to the total potential energy E can be readily obtained as E_i. However, since E_i is expressed as a linear combination of nonlinear kernel basis functions of the high-dimensional descriptor x_i, it is not immediately obvious which aspects of the local environment influence E_i. Nevertheless, by exploiting the fact that E_i is a functional of the local density distribution ρ_i, and that the descriptors are expansion coefficients of ρ_i in an orthonormal basis, as described in eqn (3), the influence of the local environment can be quantified through the variation of E_i with respect to ρ_i, which can be derived as the simple expression given below (see its derivation in Section S2 in SI):

The functional derivative δE_i/δρ_i(r) characterizes the local response of the energy E_i to an infinitesimal change in ρ_i at position r relative to atom i. With this functional, the local energy is explicitly decomposed to identify which species at which positions most strongly influence it. Notably, the exact variation in eqn (5) can be derived owing to the simple form of the local energy in eqn (4) and the expansion in a complete orthonormal basis set in eqn (3). Although eqn (5) is formulated for the kernel-based representation, a similar exact expression remains applicable to the atomic cluster expansion (ACE)¹⁹ and neural-network representation^5,9,36 of F when many-body descriptors generated from orthonormal basis sets are used as inputs.^11,12,15 We refer to this approach as local environment neighbor sensitivity (LENS) analysis. Details of the computational protocol are provided in Section S3 in SI.

We applied this analysis to cation effects on a hydroxyl group (OH) adsorbed on a Pt(111) surface in water. Prior to the LENS analysis, we computed the redox potential U_OH for the OH reduction reaction, OH* + H⁺ + e⁻ → * + H₂O, where * denotes a surface site. The calculation was performed using an MLFF-aided FP-TI scheme based on finite-temperature MD calculations. Details of this method are described in Section S4 of the SI and in ref. 15,37–40. Here, we briefly outline the approach. This approach yields the reaction free energy ΔA by integrating the derivative of Hamiltonian 〈∂H/∂λ〉_λ along a coupling path λ from the reactant state at λ = 0 to the product state at λ = 1 using MLFFs. Owing to its high computational efficiency, the MLFF approach enables statistical sampling over timescales of tens of nanoseconds, which is necessary to achieve satisfactory convergence of the reaction free energy.³⁸ However, MLFFs deviate from FP potential energy surfaces (PESs) due to incomplete descriptions of short-range interactions and the neglect of long-range interactions. These deficiencies can be systematically corrected by performing an additional TI from the MLFF to the FP-PES. Although this second TI step requires computationally demanding FPMD calculations, the energetic differences between the MLFF and FP descriptions are small (see Section S6 in the SI). As a result, the free energy difference between the two methods converges within a few picoseconds. To make the FPMD calculations feasible, it was necessary to employ a relatively small interfacial model, specifically a system consisting of 48 water molecules corresponding to a concentration of 1.2 mol L⁻¹ (see Section S5 in the SI). Although this concentrated aqueous environment may limit direct quantitative comparison with experimental results obtained at 0.1 mol L⁻¹, the two-step TI scheme corrects potential errors inherent in the MLFF while accounting for nanosecond-scale solvent reorganization dynamics, thereby enabling accurate determination of FP free energy differences.

The effects of three alkali cations, Li⁺, K⁺, and Cs⁺, on U_OH were examined using TI calculations. As a diffuse limit of the cation distribution, the same calculations were also performed for a system containing a homogeneous background charge instead of explicit cations. Following the TI calculations, LENS was applied to visualize the interactions arising from the cations.

Fig. 1 shows the computed interfacial redox potential U_OH as a function of the Shannon radius of the cation, r_ion.⁴¹ For comparison, the redox potential U_OH obtained for a system containing a homogeneous positive background charge instead of an explicit cation is also shown. The redox potential U_OH follows the order Li⁺ < K⁺ < Cs⁺ ≈ background charge, indicating a decrease in OH* stability with increasing cation size. This observation matches with the experiment.²² Another notable feature is the strong correlation between U_OH and the interfacial cation diffusion coefficient D_int (see Section S7 of the SI). As D_int increases, U_OH approaches the value for a diffusive homogeneous background charge. This trend suggests that smaller cations are trapped near OH* through attractive interactions and stabilize OH*. The localization of small cations is evident from the distributions shown in Fig. 2. The radial and three-dimensional concentration distributions indicate that Li⁺ incorporates the O atom of OH* into its first hydration shell and preferentially resides near the adsorbate, whereas the corresponding peak is significantly weaker for K⁺ and Cs⁺. None of the cations directly contact the Pt(111) surface; water molecules or OH* always separate them from the surface, consistent with FP structure optimizations by Greeley et al.²²

Fig. 1 Computed redox potential U_OH of the OH reduction, OH* + H⁺ + e⁻ → * + H₂O, versus the Shannon radius of cation r_ion. Solid line indicates computed U_OH for the system involving a homogeneous background charge. Error bars are determined from the block averaging analysis.^3,4 Graphics show the snapshots of interfacial structures obtained from the MD simulations at 300 K, and the numbers are the Bader charges of cation and OH*. Small white, medium red, medium blue, medium light green, large purple and large light green spheres show H, O in H₂O, O in OH*, Li⁺, K⁺ and Cs⁺. Graphics are obtained by using VESTA.²¹

Fig. 2 Radial distribution functions between X⁺ (X = Li, K or Cs) and O in OH*. Insets show side and top views of three-dimensional density isosurfaces of the alkaline cation at 0.002 Å⁻³ for blue and 0.007 Å⁻³ for red. Graphics are obtained by using VESTA.²¹

Bader charge⁴² and electronic structure analyses show that OH* is stably negatively charged via electron donation from Pt 5d to O 2p orbitals,⁴³ while the cations remain close to +1e, indicating that the interaction is predominantly electrostatic (see Fig. 1 and Section S8 in the SI). Although the kernel-based MLFF does not explicitly include long-range interactions, their effects are effectively incorporated into short-range interactions within the cutoff radius—such as the cation-OH* distances considered in this study—and are reproduced for the same system used in training within the deviations reported in Fig. S2.

All of the above results suggest that smaller cations interact attractively with OH* and stabilize it. Additional insight for this interaction is obtained from the LENS analysis. Fig. 3(a)–(c) visualize the functional derivative δE_i/δρ_i(r). Here, we focus on the three-body angular distribution involving H atoms, the O atom in OH*, and a cation X⁺ (X = Li, K, or Cs). The relative coordinate r is defined by the distance r between H and O, the distance s between X⁺ and O, and the angle θ = ∠HOX. In Fig. 3, the angle is fixed at 110°, which is close to the direction of the sp³ orbital of the central O atom and corresponds to a region where the functional derivative takes a low value. The LENS analysis reveals a distinct, highly localized negative response of E_i to increasing ρ_i at r ≈ 1 Å and s ≈ 2 Å for Li⁺, a broader negative response centered at r ≈ 1 Å and s ≈ 3 Å for K⁺, and an more diffuse and weaker negative response around r ≈ 1 Å and s ≈ 3.5 Å for Cs⁺. Notably, the variational response appears deeper for K⁺. However, the negative peak is spatially delocalized, leading to a correspondingly delocalized distribution of K⁺ (Fig. 2 and Section S9 in the SI) and a reduced stabilizing contribution to OH*. A rough linear-response estimate of the local energy change, , shows that Li⁺ stabilizes OH* more strongly than K⁺ [Fig. 3(d)]. Although the magnitude of this local energetic response differs from the change in U_OH, which reflects the free-energy difference of the entire system, and the difference between K⁺ and Cs⁺ cannot be resolved within this approximation, the LENS analysis nevertheless captures key features of the attractive interaction between OH* and small cations.

Fig. 3 Visualization of the functional derivative δE_i/δρ_i obtained using LENS (a)–(c) and the stabilization energy ΔE_i^lin estimated from the linear-response approximation (see text) (d). Here, E_i is the energy of the O atom in OH*, and ρ_i denotes the probability density of finding an H atom at a distance r to r + dr from this O atom and a cation X⁺ (X = Li, K, or Cs) at a distance s to s + ds from the same O atom, with ∠ HOX = θ = 110°, as illustrated in the insets. LENS yields smooth functional-derivative landscapes below ca. 1 Å for Li⁺ and 2 Å for the other ions; however, these regions are not sampled in the MD trajectories. The values in these regions arise from kernel extrapolation and are therefore physically meaningless.

In summary, our MLFF-aided TI scheme combined with the interaction visualization method LENS shows that small alkaline cations are trapped by OH adsorbed on Pt(111) in water and stabilize this interfacial species via attractive interactions, consistent with experiments.²² Notably, the kernel-based representation of interatomic interactions, which does not explicitly encode physical laws, learns from FP data to capture the attractive interaction between small cations and OH without human intervention. Furthermore, LENS reveals information beyond the total-energy summation in eqn (4) by decomposing atomic energy contributions in terms of local species distributions. These results show that MLFFs combined with LENS can non-empirically identify species that stabilize or destabilize atoms or atomic groups in complex heterogeneous environments, highlighting their potential for materials design.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). The supplementary information includes derivation of eqn (5), details of models and computational procedures. The data and LENS code used for this work are provided in a zipped file to facilitate reproduction of the results. See DOI: https://doi.org/10.1039/d6cc00436a.

References

1. R. M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, Cambridge, 2004.
2. D. Marx and J. Hutter, Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods, Cambridge University Press, 2009.
3. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2017.
4. D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press, Inc., USA, 2nd edn, 2001.
5. J. Behler, J. Chem. Phys., 2016, 145, 170901.
6. V. L. Deringer, M. A. Caro and G. Csányi, Adv. Mater., 2019, 31, 1902765.
7. R. Jinnouchi, K. Miwa, F. Karsai, G. Kresse and R. Asahi, J. Phys. Chem. Lett., 2020, 11, 6946–6955.
8. T. W. Ko, J. A. Finkler, S. Goedecker and J. Behler, Nat. Commun., 2021, 12, 398.
9. J. Behler, Chem. Rev., 2021, 121, 10037–10072.
10. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Chem. Rev., 2021, 121, 10142–10186.
11. F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi and M. Ceriotti, Chem. Rev., 2021, 121, 9759–9815.
12. V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti and G. Csányi, Chem. Rev., 2021, 121, 10073–10141.
13. D. M. Anstine and O. Isayev, J. Phys. Chem. A, 2023, 127, 2417–2431.
14. I. Batatia, S. Batzner, D. P. Kovács, A. Musaelian, G. N. C. Simm, R. Drautz, C. Ortner, B. Kozinsky and G. Csányi, Nat. Mach. Intell., 2025, 7, 56–67.
15. R. Jinnouchi and S. Minami, ACS Nano, 2025, 19, 22600–22644.
16. J. E. Mayer, J. Chem. Phys., 1937, 5, 67–73.
17. A. P. Bartók, R. Kondor and G. Csányi, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 184115.
18. A. V. Shapeev, Multiscale Model. Simul., 2016, 14, 1153–1173.
19. R. Drautz, Phys. Rev. B, 2019, 99, 014104.
20. A. Glielmo, C. Zeni and A. De Vita, Phys. Rev. B, 2018, 97, 184307.
21. K. Momma and F. Izumi, J. Appl. Crystallogr., 2008, 41, 653–658.
22. D. Strmcnik, K. Kodama, D. van der Vliet, J. Greeley, V. R. Stamenkovic and N. M. Marković, Nat. Chem., 2009, 1, 466–472.
23. A. Goyal, S. Louisia, P. Moerland and M. T. M. Koper, J. Am. Chem. Soc., 2024, 146, 7305–7312.
24. M. S. Shibata, Y. Morimoto, A. Z. Weber and I. V. Zenyuk, ACS Energy Lett., 2025, 10, 3922–3930.
25. J. Resasco, L. D. Chen, E. Clark, C. Tsai, C. Hahn, T. F. Jaramillo, K. Chan and A. T. Bell, J. Am. Chem. Soc., 2017, 139, 11277–11287.
26. S. Ringe, E. L. Clark, J. Resasco, A. Walton, B. Seger, A. T. Bell and K. Chan, Energy Environ. Sci., 2019, 12, 3001–3014.
27. D. T. Limmer, A. P. Willard, P. Madden and D. Chandler, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4200–4205.
28. I. T. McCrum and M. J. Janik, J. Phys. Chem. C, 2016, 120, 457–471.
29. H. H. Kristoffersen, K. Chan, T. Vegge and H. A. Hansen, Chem. Commun., 2020, 56, 427–430.
30. R. Jinnouchi, J. Lahnsteiner, F. Karsai, G. Kresse and M. Bokdam, Phys. Rev. Lett., 2019, 122, 225701.
31. R. Jinnouchi, F. Karsai and G. Kresse, Phys. Rev. B, 2019, 100, 014105.
32. R. Jinnouchi, F. Karsai, C. Verdi, R. Asahi and G. Kresse, J. Chem. Phys., 2020, 152, 234102.
33. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
34. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
35. A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Phys. Rev. Lett., 2010, 104, 136403.
36. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 98, 146401.
37. R. Jinnouchi, F. Karsai and G. Kresse, npj Comput. Mater., 2024, 10, 107.
38. R. Jinnouchi, J. Chem. Phys., 2024, 161, 194110.
39. R. Jinnouchi, F. Karsai and G. Kresse, Chem. Sci., 2025, 16, 2335–2343.
40. R. Jinnouchi and S. Minami, J. Phys. Chem. Lett., 2025, 16, 265–273.
41. R. D. Shannon, Acta Crystallogr., Sect. A, 1976, 32, 751–767.
42. W. Tang, E. Sanville and G. Henkelman, J. Phys.: Condens. Matter, 2009, 21, 084204.
43. M. Arce, P. Quaino and E. Santos, Catal. Today, 2013, 202, 120–127.

---