Modelling complex proton transport phenomena -- Exploring the limits of fine-tuning and transferability of foundational machine-learned force fields

Abstract

The solid acids CsH$_2$PO$_4$ and Cs$_7$(H$_4$PO$_4$)(H$_2$PO$_4$)$_8$ pose significant challenges for the simulation of proton transport phenomena. In this work, we use the recently developed machine-learned force field (MLFF) MACE to model the proton dynamics on nanosecond time scales for these systems and compare its performance with long-term ab initio molecular dynamics (AIMD) simulations. The MACE-MP-0 foundation model shows remarkable performance for all observables derived from molecular dynamics (MD) simulations, but minor quantitative discrepancies remain compared to the AIMD reference data. However, we show that minimal fine-tuning -- fitting to as little as 1 ps of AIMD data -- leads to full quantitative agreement between the radial distribution functions of MACE force field and AIMD simulations. In addition, we show that traditional long-term AIMD simulations fail to capture the correct qualitative trends in diffusion coefficients and activation energies for these solid acids due to the limited accessible time scale. In contrast, accurate and convergent diffusion coefficients can be reliably obtained through multi-nanosecond long MD simulations using machine-learned force fields. The obtained qualitative and quantitative behavior of the converged diffusion coefficients and activation energies now matches the experimental trends for both solid acids, in contrast to previous AIMD simulations that yielded a qualitatively wrong picture.

Figures (3)

• Figure 1: a Snapshot from an AIMD simulation of CDP. Color scheme for atoms: hydrogen in blue, oxygen in red, phosphorus in green, and caesium in pink. b Comparison of the O-H radial distribution function $g(r)$ obtained from different MACE models (see main text) with $g(r)$ from an AIMD simulation. The inset highlights the peak of $g(r)$ at $d_\mathrm{OH} = 1.5$ Å, which is commonly referred to as the ”short strong hydrogen bond” or ”low/barrier hydrogen bond” bum2004limbach2009whiteley2017dereka2021. c Histogram showing the similarity coefficients SC as defined by Eq. (\ref{['eqn:sc']}) with respect to the AIMD results, i.e., $\mathrm{SC}[\mathrm{MACE};\,\mathrm{AIMD}]$ for all possible bond combinations for the three selected MACE models from b, see legend. The average SC is indicated by vertical dashed lines in the corresponding color.
• Figure 2: a Snapshot from an AIMD simulation of CPP. Anionic H$_2$PO$_4^-$ tetrahedra are highlighted in green, and cationic H$_4$PO$_4^+$ tetrahedra are highlighted in yellow. Color scheme for atoms: hydrogen in blue, oxygen in red, phosphorus in green/yellow (anionic/cationic), and caesium in pink. b Comparison of the anionic O-H radial distribution function $g(r)$ obtained from different MACE models (see main text) with $g(r)$ from an AIMD simulation. The inset highlights the peak of $g(r)$ at $d_\mathrm{OH} = 1.5$ Å, which is commonly referred to as the ”short-strong hydrogen bond” or ”low-barrier hydrogen bond”. c Same as b but for the cationic O-H radial distribution function.
• Figure 3: Visualization of the convergence behavior of the hydrogen diffusion coefficient with respect to simulation time for CDP and CPP. AIMD simulations at 513 K and fully fine-tuned mace simulations at 510 K are compared. Convergence of the diffusion coefficient is achieved only for timescales above 5 ns, which are not accessible by AIMD simulations. The overall diffusion coefficient (represented by solid and dashed lines) is computed as the average of the diffusion coefficients of all individual protons. The shaded area around each line represents the standard deviation calculated for these sets of individual diffusion coefficients.