# Accelerating QM/MM simulations of electrochemical interfaces through machine
learning of electronic charge densities

Andrea Grisafi<EMAIL_ADDRESS>Institut Sciences du
Calcul et des Données, ISCD, Sorbonne Université, F-75005 Paris, France
Mathieu Salanne<EMAIL_ADDRESS>Physicochimie des
Électrolytes et Nanosystèmes Interfaciaux, Sorbonne Université, CNRS, F-75005
Paris, France  Institut Universitaire de France (IUF), F-75231 Paris, France

###### Abstract

A crucial aspect in the simulation of electrochemical interfaces consists in
treating the distribution of electronic charge of electrode materials that are
put in contact with an electrolyte solution. Recently, it has been shown how a
machine-learning method that specifically targets the electronic charge
density, also known as SALTED, can be used to predict the long-range response
of metal electrodes in model electrochemical cells. In this work, we provide a
full integration of SALTED with MetalWalls, a program for performing classical
simulations of electrochemical systems. We do so by deriving a spherical
harmonics extension of the Ewald summation method, which allows us to
efficiently compute the electric field originated by the predicted electrode
charge distribution. We show how to use this method to drive the molecular
dynamics of an aqueous electrolyte solution under the quantum electric field
of a gold electrode, which is matched to the accuracy of density-functional
theory. Notably, we find that the resulting atomic forces present a small
error of the order of 1 meV/Å, demonstrating the great effectiveness of
adopting an electron-density path in predicting the electrostatics of the
system. Upon running the data-driven dynamics over about 3 ns, we observe
qualitative differences in the interfacial distribution of the electrolyte
with respect to the results of a classical simulation. By greatly accelerating
quantum-mechanics/molecular-mechanics approaches applied to electrochemical
systems, our method opens the door to nanoseconds timescales in the accurate
atomistic description of the electrical double layer.

##  I Introduction

Understanding the interface between a metallic electrode and a liquid
electrolyte has become an important objective in chemical physics. Dawlaty
_et al._ (2023) This interest is driven by the many underlying applications in
electrochemistry, electrocatalysis, and energy storage, but it is also a
challenging problem from a fundamental point of view. Most of the time, the
conventional simulation methods cannot easily be adopted due to the highly
heterogeneous nature of the system, the difficulties associated with the
application of a voltage between two electrodes, and the large time and length
scales required to simulate a representative portion of the electrochemical
interface. Schwarz and Sundararaman (2020); Groß (2023)

In the absence of chemical reactions, intermolecular interactions in bulk
electrolytes are often well handled using molecular mechanics (MM) approaches.
Goloviznina  _et al._ (2019) The main difficulty is therefore to account for
the electronic structure of the metal. For this purpose, performing quantum
mechanics (QM) calculations, most notably at the density functional theory
(DFT) level, appears as the most accurate available approach. Although it is
possible to use DFT to calculate atomic forces within a molecular dynamics
(MD) simulation, the large computational cost limits its application to
relatively simple systems. Li  _et al._ (2019); Khatib  _et al._ (2021); Chen
_et al._ (2022) Typically, having an electrolyte including a salt at finite
concentration remains a hurdle for such simulations. A possible solution is to
use a purely classical approach, which leads to the family of constant-
potential MD simulations. Jeanmairet, Rotenberg, and Salanne (2022) In this
case, the polarization of the metal is treated using fluctuating atomic
Gaussian charges, which depend on both the electrolyte positions and the
applied potential. This approach was refined to effectively include some
electronic structure effects, Scalfi  _et al._ (2020a); Goloviznina  _et al._
(2024) but the electrode representation remains highly simplified when
compared to an explicit QM-based method.

An alternative is to mix the QM/MM approaches by partitioning the system. In
fact, by performing the QM calculation on the electrode subsystem only, it is
possible to greatly decrease the computational burden. Such calculations were
for example used to study supercapacitors, Elliott  _et al._ (2023) the
electrowetting of graphene, Wei  _et al._ (2024) or the electrocatalytic
reactivity of MoS2-based materials. Abidi and Steinmann (2023) There are
however several limitations. Firstly, performing a QM calculation at every
timestep of a MD simulation ($\approx$ 1fs) is not very efficient, so that
multi-stepping techniques have to be introduced by updating the QM system
every few thousands steps. Secondly, while the energy of the system can be
recovered by adequately summing up the energies of the sub-systems, the same
is not always true for the forces. In particular, calculating the
electrostatic forces acting on the MM subsystem due to the QM subsystem is not
straightforward. The problem is generally tackled by introducing equilibrated
partial charges on the electrode atoms, Takahashi, Nakano, and Sato (2022)
which are computed following a quadratic functional form of the electrostatic
energy. Beside the assumed quadratic approximation, this type of charge
partitioning is not rigorously defined, so that some important physical
effects associated with the local anisotropy of the charge distribution may be
lost in the process.

In this work, we introduce a data-driven method that allows us to improve both
the efficiency and the accuracy of QM/MM simulations of electrochemical
interfaces. This method is grounded on recent developments in a machine-
learning (ML) approach capable of predicting the electronic charge density of
metal electrodes, including its non-local response to far-field perturbations.
Grisafi  _et al._ (2023a) In particular, we implement an interface between the
SALTED ML program Grisafi and Lewis (2023) and the MetalWalls MD code for
electrochemical simulations. Coretti  _et al._ (2022) As a key aspect of this
interface, we provide an efficient calculation of the electric field generated
by the predicted electronic charge distribution, which we derive following a
spherical harmonics generalization of the Ewald summation method. In turn,
this allows us to compute accurate electrostatic forces acting on the
electrolyte atoms, positioning our method in stark contrast with approaches
that are based on partial charges. By shifting the largest share of the
computational costs on the learning stage, simulation times of the order of
nanoseconds can be achieved with a relatively small overhead compared to
conventional full MM simulations. As an application of our method, we study
the interface between a gold electrode and an aqueous solution of NaCl at 5 M
concentration under applied electric potentials of 0 and $1$ V. The structure
of the liquid phase as well as interfacial properties such as the capacitance
are calculated and compared to previous approaches.

##  II General method

We adopt a finite-field simulation setup, where a single electrode is put in
contact with an electrolyte solution under three-dimensional periodic boundary
conditions. This setup can be conveniently used to simulate a pair of
equivalent electrodes that are kept at a constant potential difference. Dufils
_et al._ (2019) In particular, a uniform and constant electric field,
$\varepsilon_{z}$, is applied along the direction normal to the electrode
surface, such that $\Delta V=-\varepsilon_{z}\,L_{z}$, with $L_{z}$ the length
of the simulation box, corresponds to a physical electrochemical cell
potential. The representation of the electrode/electrolyte interface follows
the approach already presented in Ref. 15. Specifically, we adopt a quantum-
mechanics/molecular-mechanics (QM/MM) strategy, where an ab initio
representation of the electrode is used in combination with a classical
representation of the electrolyte defined via fixed Gaussian charges. We aim
to rely on SALTED predictions of the electrode charge density in order to
compute accurate electrostatic forces that can be used to drive the dynamics
of the system.

###  II.1 Representation of the electrode charge density

Within SALTED, the electronic charge density of the electrode is represented
as a linear spherical harmonics expansion that follows density-fitting (DF)
approximations commonly adopted in electronic-structure methods: Grisafi  _et
al._ (2023b); Lewis  _et al._ (2021)

$n_{e}(\boldsymbol{r})=\sum_{inlm}c^{nlm}_{i}\sum_{\boldsymbol{u}}R_{nl}(\left|\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}\right|)Y^{l}_{m}(\widehat{\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}})\,.$
(1)

Here, the index $i$ identifies the electrode atoms, $nlm$ the basis function
indexes, while $\boldsymbol{u}$ represent the three-dimensional translation
vectors over the periodic images of the system. We assume that radial
auxiliary functions are of Gaussian-type, i.e., $R_{nl}\propto
r^{l}e^{-r^{2}/(2\sigma^{2}_{n})}$, while $Y^{l}_{m}$ are defined as
orthonormalized real spherical harmonics.

Reference calculations of $n_{e}$ in the presence of the classical electrolyte
field are performed using the CP2K simulation program. Kühne  _et al._ (2020)
Specifically, quantum-mechanical calculations are performed at the Kohn-Sham
DFT level using the PBE functional Perdew, Burke, and Ernzerhof (1996) with
DZVP-MOLOPT-SR basis sets VandeVondele and Hutter (2007) and GTH pseudo-
potentials. Goedecker, Teter, and Hutter (1996) To extract the expansion
coefficients, $c^{nlm}_{i}$, we adopt a density-fitting approximation that
makes use of a truncated Coulomb metric, Guidon, Hutter, and VandeVondele
(2009); Bussy, Schütt, and Hutter (2023) a choice that is especially suitable
to represent the electrostatic properties of the system.

The total electrode charge density is obtained by adding to $n_{e}$ the
effective nuclear-charge density associated with the classical local part of
the GTH pseudopotential. For each electrode atom $i$, this is simply defined
by an isotropic Gaussian charge, i.e.,

$n_{\text{ion}}(r)=\frac{Z_{\text{eff}}}{\sqrt{8}\,\pi^{3/2}\,r^{3}_{\text{loc}}}\exp\left[-\frac{1}{2}\left(\frac{r}{r_{\text{loc}}}\right)^{2}\right]\,,$
(2)

where $Z_{\text{eff}}$ is the effective nuclear charge of the electrode atoms
and $r_{\text{loc}}$ is the local pseudopotential cutoff.

###  II.2 SALTED electrode model

In analogy with the prescription adopted in Ref. 15, we will assume that the
positions of the electrode atoms are kept fixed, so that the learning can be
focused on the sole electrolyte-induced variations of the electrode charge
density. Moreover, as the applied voltage $\Delta V$ is typically of the order
of a few volts over a cell that covers several nanometers along $z$, we can
further assume that the external electric field, defined as
$\varepsilon_{z}\equiv-\Delta V/L_{z}$, introduces a small linear perturbation
to $n_{e}$. Under these assumptions, we can conveniently define the SALTED
prediction target as the following density difference,

$\Delta
n_{e}(\boldsymbol{r})=n_{e}(\boldsymbol{r};\varepsilon_{z})-n^{0}_{e}(\boldsymbol{r};\varepsilon_{z})\,,$
(3)

where $n^{0}_{e}$ is the electron density of the isolated electrode under
$\varepsilon_{z}$, used as a constant baseline for the regression procedure.
Upon training, $n^{0}_{e}$ is eventually added back to the predicted $\Delta
n_{e}$ in order to recover the total electron density $n_{e}$. Notably, the
possibility of treating the external field as a small perturbation implies
that $\Delta n_{e}$ can be considered independent from $\varepsilon_{z}$. From
a ML point of view, this represents a great advantage, as it allows us to
simulate a continuum spectrum of cell potentials without the need to retrain a
new SALTED model whenever a different value of $\varepsilon_{z}$ is applied.

Consistently with the density-fitting metric used to generate the reference
quantum-mechanical data, SALTED models can be trained by minimizing a loss
function that makes use of a (truncated) Coulomb metric. Briling, Fabrizio,
and Corminboeuf (2021) In compact notation, this reads as follows,

$\begin{split}\mathcal{L}(\boldsymbol{w})&=\sum_{I=1}^{N}\bra{n^{\text{ML}}_{e}(\boldsymbol{w})-n^{\text{DF}}_{e}}\frac{1}{\boldsymbol{r}-\boldsymbol{r}^{\prime}}\ket{n^{\text{ML}}_{e}(\boldsymbol{w})-n^{\text{DF}}_{e}}_{I}\\\
&=\sum_{I=1}^{N}\left[\boldsymbol{c}^{\text{ML}}_{I}(\boldsymbol{w})-\boldsymbol{c}^{\text{DF}}_{I}\right]^{T}\mathbf{J}_{I}\left[\boldsymbol{c}^{\text{ML}}_{I}(\boldsymbol{w})-\boldsymbol{c}^{\text{DF}}_{I}\right]\,,\end{split}$
(4)

where $N$ is the number of training configurations. In the second line of this
equation, $\mathbf{J}$ is the matrix of 2-centers Coulomb integrals between
the auxiliary functions, i.e.,
$\bra{i^{\prime}n^{\prime}l^{\prime}m^{\prime}}r^{-1}\ket{inlm}$, as computed
by CP2K, $\boldsymbol{c}^{\text{DF}}$ and $\boldsymbol{c}^{\text{ML}}$ are the
density-fitting and machine-learning vectors of density coefficients
associated with the expansion of $\Delta n_{e}$, respectively, and
$\boldsymbol{w}$ is the vector of regression weights. In practice, a
regularization term is also added to the loss function to avoid overfitting;
we refer to Ref. 19 for an in-depth discussion of the SALTED method.

Following Ref. 15, the symmetry-adapted kernels used for the approximation of
the $\boldsymbol{c}^{\text{ML}}$ density coefficients can be defined to
include long-range information. In particular, kernel functions between atomic
environments are constructed from long-distance equivariant (LODE) structural
descriptors, Grisafi and Ceriotti (2019); Grisafi, Nigam, and Ceriotti (2021);
Huguenin-Dumittan  _et al._ (2023) as implemented in the “rascaline” atomistic
representations package. Fraux  _et al._ This feature of the method is
fundamental to capture the non-local variations of the electron density
induced by electrolyte charges that are located arbitrarily far from the
electrode surface. Grisafi  _et al._ (2023a) We remark that, unlike charge
equilibration schemes, this is made possible without relying on any self-
consistent procedure. In fact, the use of long-range descriptors allows us to
learn directly the self-consistent electronic polarization of the system, thus
bypassing the need to incorporate self-consistency at the ML level. This would
still hold true in the scenario where the electrolyte is also undergoing an
electronic polarization.

##  III Electric field calculation

From the predicted electrode charge distribution, a nontrivial problem is that
of performing an efficient calculation of the electric field needed to drive
the dynamics of the electrolyte. To achieve this goal, we provide in this
section an original analytical derivation of the Ewald summation method,
properly extended to account for the spherical harmonics expansion of the
electron density.

Following the representation of Eq. (1), we start by introducing a fictitious
electron density, $n^{\text{ewald}}_{e}$, obtained by combining the physical
density coefficients, $c^{nlm}_{i}$, with a new set of smooth radial
functions, $R_{nl}^{\text{ewald}}$. Crucially, these functions are defined
from a Gaussian width $\sigma_{\text{ewald}}$ that is larger than those
entering the definition of $R_{nl}$. Following the conventional Ewald
summation method, we then separate from the electron density a charge-neutral
(screened) contribution, defined as

$n^{\text{screen}}_{e}(\boldsymbol{r})\equiv
n_{e}(\boldsymbol{r})-n^{\text{ewald}}_{e}(\boldsymbol{r})\,.$ (5)

By construction, the electric field generated by $n^{\text{screen}}_{e}$ is
short-ranged, and can be efficiently computed in real space. Conversely, the
remaining (unscreened) $n^{\text{ewald}}_{e}$ density will contribute to a
long-range electric field that can be efficiently computed in reciprocal
space. In what follows, we outline the derivation of both electronic and
nuclear contributions to the electric field, as newly implemented in the
MetalWalls simulation program.

###  III.1 Short-range electronic contribution

Starting from $n^{\text{screen}}_{e}$, the short-range part of the Hartree
electronic potential can be analytically computed in real space by relying on
the Laplace expansion of the Coulomb operator. After some calculations,
detailed in the Supplementary Material, we obtain

$\begin{split}&\phi^{\text{SR}}_{\text{H}}(\boldsymbol{r})=\int_{\mathbf{R}^{3}}d\boldsymbol{r}^{\prime}\,\frac{n^{\text{screen}}_{e}(\boldsymbol{r}^{\prime})}{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}\\\
&=4\pi\sum_{i=1}^{N_{\text{at}}}\sum_{nlm}\frac{c_{i}^{nlm}}{2l+1}\sum_{\boldsymbol{u}}Y_{lm}(\widehat{\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}})\\\
&\times\left[\frac{I^{<}_{nl}(|\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}|)}{|\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}|^{l+1}}+|\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}|^{l}I^{>}_{nl}(|\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}|)\right]\,.\end{split}$
(6)

Here, $I_{nl}$ indicate radial integrals performed within ($<$) and outside
($>$) a sphere of radius
$r\equiv|\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}|$:

$\begin{split}&I_{ln}^{<}(r)=\int_{0}^{r}ds\,s^{2+l}\left[R_{nl}(s)-R^{\text{ewald}}_{nl}(s)\right]\,,\\\
&I_{ln}^{>}(r)=\int_{r}^{\infty}ds\,s^{1-l}\left[R_{nl}(s)-R^{\text{ewald}}_{nl}(s)\right]\,.\end{split}$
(7)

While the individual terms of $I_{ln}^{>}$ are quickly vanishing for
$r\to\infty$, the requirement that $I_{ln}^{<}$ vanishes sufficiently fast for
increasing values of $r$ implies that the normalization of the Ewald radial
functions must be chosen so that to satisfy the following equality:

$\int_{0}^{\infty}ds\,s^{2+l}R^{\text{ewald}}_{nl}(s)=\int_{0}^{\infty}ds\,s^{2+l}R_{nl}(s)\,.$
(8)

This choice guarantees that both types of radial integrals in Eq. (7) tend to
zero with superexponential rapidity, thus dominating over the algebraic terms
of the order of $1/r^{l+1}$ and $r^{l}$ in Eq. (LABEL:eq:ewald-pot-real).
Crucially, this screening effect allows us to compute the short-range part of
the Hartree potential by only collecting contributions of electrode atoms that
fall within a finite cutoff distance. In particular, we select a cutoff of
$r_{\text{cut}}^{l}=4\sqrt{2+l}\,\sigma_{\text{Ewald}}$, which is defined to
account for the spatial range covered by the Ewald Gaussian-type functions.

To compute the electric field as
$\boldsymbol{E}^{\text{SR}}_{\text{H}}=-\nabla\phi^{\text{SR}}_{\text{H}}$, it
is convenient to adopt the substitution
$\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}\equiv\left(r,\theta,\phi\right)$
in Eq. (LABEL:eq:ewald-pot-real), so that to make use of the gradient in
spherical coordinates, i.e.,

$\tilde{\nabla}\equiv\left(\frac{\partial}{\partial
r},\frac{1}{r}\frac{\partial}{\partial\theta},\frac{1}{r\sin\theta}\frac{\partial}{\partial\phi}\right)\,.$

Upon this change of variables, the contribution to the short-range part of the
electric field brought by each electrode atom can be analytically computed. In
compact notation, we have

$\begin{split}\boldsymbol{E}^{\text{SR}}_{\text{H}}(\boldsymbol{r})&=4\pi\sum_{i=1}^{N_{\text{at}}}\sum_{nlm}\frac{c_{i}^{nlm}}{2l+1}\,\sum_{\boldsymbol{u}}\mathbf{U}(\widehat{\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}})\\\
&\cdot\tilde{\nabla}\left[\left(\frac{1}{r^{l+1}}\,I_{ln}^{<}(r)+r^{l}\,I_{ln}^{>}(r)\right)Y^{l}_{m}(\theta,\phi)\right]\,,\end{split}$
(9)

where

$\mathbf{U}=\begin{pmatrix}\sin\theta\cos\phi&\cos\theta\cos\phi&-\sin\phi\\\
\sin\theta\sin\phi&\cos\theta\sin\phi&\cos\phi\\\
\cos\theta&-\sin\theta&0\end{pmatrix}\,,$ (10)

represents the unitary transformation matrix between the Cartesian and
spherical unit vectors
$(\hat{x},\hat{y},\hat{z})\to(\hat{r},\hat{\theta},\hat{\phi})$, such that
$\nabla=\mathbf{U}\cdot\tilde{\nabla}$ recovers the Cartesian gradient.
Explicit analytical expressions of radial integrals and spherical harmonics
derivatives are reported in the Supplementary Material.

###  III.2 Long-range electronic contribution

From a smooth choice of $n^{\text{ewald}}_{e}$, the long-range part of the
Hartree potential can be conveniently written as a plane-waves expansion
compatible with the three-dimensional periodicity of the system:

$\begin{split}\phi^{\text{LR}}_{\text{H}}(\boldsymbol{r})&=\sum_{\boldsymbol{k}\neq\boldsymbol{0}}\frac{4\pi}{k^{2}}\,\tilde{n}^{\text{ewald}}_{e}(\boldsymbol{k})\,e^{i\boldsymbol{k}\cdot\boldsymbol{r}}\,.\end{split}$
(11)

Note that we excluded the $\boldsymbol{k}=0$ term because of the assumed
charge neutrality of the overall system, i.e., electrons and nuclei. From the
expansion of the plane wave in spherical harmonics, the Fourier components of
the electron density are found to get the form of

$\begin{split}\tilde{n}^{\text{ewald}}_{e}(\boldsymbol{k})&=\frac{1}{\Omega}\int_{\Omega}d\boldsymbol{r}\,e^{-i\boldsymbol{k}\cdot\boldsymbol{r}}n^{\text{ewald}}_{e}(\boldsymbol{r})\\\
&=\frac{4\pi}{\Omega}\sum_{i=1}^{N_{\text{at}}}e^{-i\boldsymbol{k}\cdot\boldsymbol{r}_{i}}\sum_{nlm}c^{nlm}_{i}(-i)^{l}Y^{l}_{m}(\hat{\boldsymbol{k}})\,\hat{I}_{nl}(k)\,,\end{split}$
(12)

where $\Omega$ is the cell volume. Here, $\hat{I}_{nl}$ are defined as the
following radial integrals,

$\begin{split}\hat{I}_{nl}(k)&=\int_{0}^{\infty}dr\,r^{2}j_{l}(kr)R^{\text{ewald}}_{nl}(r)\,,\end{split}$
(13)

with $j_{l}$ a spherical Bessel function of the first kind. Upon plugging Eq.
(12) into Eq. (11), the final analytical expression for the long-range Hartree
potential is

$\begin{split}\phi^{\text{LR}}_{\text{H}}(\boldsymbol{r})&=\frac{32\pi^{2}}{\Omega}\sum^{\text{half}}_{\boldsymbol{k}\neq\boldsymbol{0}}\frac{1}{k^{2}}\sum_{i=1}^{N_{\text{at}}}\sum_{nlm}c^{nlm}_{i}\,\hat{I}_{nl}(k)\,Y^{l}_{m}(\hat{\boldsymbol{k}})\\\
&\times\begin{cases}(-1)^{\frac{l}{2}}\cos\left[\boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{r}_{i})\right]\,,\
l\%2=0\\\
(-1)^{\frac{l+3}{2}}\sin\left[\boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{r}_{i})\right]\,,\
l\%2\neq 0\end{cases}\,.\end{split}$ (14)

Note that we exploited the real nature of the electron density to sum
$\boldsymbol{k}$-vectors only over a semisphere in reciprocal space of radius
$k_{\text{cut}}$. In practice, $k_{\text{cut}}$ is indirectly defined by the
chosen value of $\sigma_{\text{Ewald}}$, following the numerical convergence
criterion implemented in MetalWalls. Marin-Laflèche  _et al._ (2020)

Upon Cartesian differentiation of Eq. (14), the long-range part of the
electronic electric field is finally given by

$\begin{split}\boldsymbol{E}^{\text{LR}}_{\text{H}}(\boldsymbol{r})=&\frac{32\pi^{2}}{\Omega}\sum^{\text{half}}_{\boldsymbol{k}\neq\boldsymbol{0}}\frac{\boldsymbol{k}}{k^{2}}\sum_{i=1}^{N_{\text{at}}}\sum_{nlm}c^{nlm}_{i}\,\hat{I}_{nl}(k)\,Y^{l}_{m}(\hat{\boldsymbol{k}})\\\
&\times\begin{cases}(-1)^{\frac{l}{2}}\,\sin\left[\boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{r}_{i})\right]\,,\
l\%2=0\\\
(-1)^{\frac{l+1}{2}}\,\cos\left[\boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{r}_{i})\right]\,,\
l\%2\neq 0\end{cases}\,.\end{split}$ (15)

An explicit expression of the radial integrals, together with additional
implementation details, are reported in the Supplementary Material.

###  III.3 Pseudopotential nuclear contribution

From Eq. (2), computing the nuclear contribution to the electric field is well
established, Marx and Hutter (2009) and it reduces to the problem of treating
a set of classical Gaussian charges. Following the conventional Ewald method,
this amounts to compute a screened short-range potential as

$\phi^{\text{SR}}_{\text{ion}}(r)=\frac{Z_{\text{eff}}}{r}\left[\text{erf}\left(\frac{r}{\sqrt{2}r_{\text{loc}}}\right)-\text{erf}\left(\frac{r}{\sqrt{2}\sigma_{\text{ewald}}}\right)\right]\,,$
(16)

and a long-range contribution in reciprocal space, whose Fourier components
are given by

$\tilde{\phi}^{\text{LR}}_{\text{ion}}(k)=\frac{4\pi}{\Omega}\frac{Z_{\text{eff}}}{k^{2}}\exp\left(-\frac{1}{2}k^{2}\sigma^{2}_{\text{ewald}}\right)\
.$ (17)

For consistency, we select a value of $\sigma_{\text{ewald}}$ equivalent to
that adopted for the electronic field calculation.

The short-range part of the electric field can be computed by relying on the
same transformation matrix between Cartesian and spherical unit vectors
already introduced in Eq. (10). Considering that Eq. (16) has no angular
components, we can write

$\begin{split}\boldsymbol{E}^{\text{SR}}_{\text{ion}}(\boldsymbol{r})=\sum_{i=1}^{N_{\text{at}}}\sum_{\boldsymbol{u}}\mathbf{U}(\widehat{\boldsymbol{r}-\boldsymbol{r}_{i}-\boldsymbol{u}})\cdot\left[-\frac{\partial}{\partial
r}\phi^{\text{SR}}_{\text{ion}}(r),0,0\right]\,.\end{split}$ (18)

The long-range part of the nuclear electric field is finally computed as
follows,

$\begin{split}\boldsymbol{E}^{\text{LR}}_{\text{ion}}(\boldsymbol{r})&=2\sum^{\text{half}}_{\boldsymbol{k}\neq\boldsymbol{0}}\boldsymbol{k}\,\tilde{\phi}^{\text{LR}}_{\text{ion}}(k)\sum_{i=1}^{N_{\text{at}}}\sin\left[\boldsymbol{k}\cdot(\boldsymbol{r}-\boldsymbol{r}_{i})\right]\,.\end{split}$
(19)

##  IV SALTED/MetalWalls interface

Upon implementing the electric field just derived into the MetalWalls
molecular dynamics engine, the SALTED/MetalWalls simulation can be easily run
thanks to a Python-based interface between the two programs. In what follows,
we report a description of the simulation workflow, including pseudocode of
the relevant steps.

  1. 1.

Import SALTED and MetalWalls modules:

import salted, metalwalls

  2. 2.

Load pretrained SALTED information:

salted_info = salted.init_pred.build()

  3. 3.

Load isolated electrode information:

structure = ase.io.read("init.xyz")
coords_electrode = structure.get_positions()[640:]
coefs_0 = np.loadtxt("coefs_isolated.txt")

  4. 4.

Define MetalWalls variables.

my_system =
metalwalls.mw_system.MW_system_t()
my_algorithms =
metalwalls.mw_algorithms.MW_algorithms_t()
my_parallel =
metalwalls.mw_parallel.MW_parallel_t()

  5. 5.

Initialize electrode charge density:

coords_electrolyte = my_system.xyz_ions
coords = np.vstack((coords_electrolyte,
coords_electrode))
structure.set_positions(coords)
my_system.coefs = salted.salted_prediction
.build(salted_info,structure)
my_system.coefs += coefs_0

  6. 6.

Initialize MetalWalls simulation:

metalwalls.mw_tools.run_step(my_system, my_parallel, my_algorithms, 0, 0,
do_output, step_output_frequency)

  7. 7.

Start the dynamics:

for step in range(my_system.num_steps):

  8. 8.

Update electrolyte atomic positions:

metalwalls.mw_tools.step_setup(my_system, my_algorithms, my_parallel, step)
coords_electrolyte = my_system.xyz_ions
coords = np.vstack((coords_electrolyte,
coords_electrode))
structure.set_positions(coords)

  9. 9.

Update electrode charge density:

my_system.coefs = salted.salted_prediction
.build(salted_info,structure)
my_system.coefs += coefs_0

  10. 10.

Compute atomic forces:

metalwalls.mw_tools.step_compute_forces
(my_system, my_parallel, step)

  11. 11.

Continue from point 8.

Importantly, the electrostatic forces acting on the electrolyte atoms $j$ are
computed at step 10 as follows,

$\boldsymbol{f}_{j}=q_{j}\,\boldsymbol{E}(\boldsymbol{r}_{j})\,,$ (20)

where $q_{j}$ are the electrolyte partial atomic charges, while

$\boldsymbol{E}=\boldsymbol{E}^{\text{SR}}_{\text{H}}+\boldsymbol{E}^{\text{LR}}_{\text{H}}+\boldsymbol{E}^{\text{SR}}_{\text{ion}}+\boldsymbol{E}^{\text{LR}}_{\text{ion}}\,,$
(21)

is the total electric field generated by the predicted electrode charge
distribution, as decomposed from the derivation carried out in Sec. III.

##  V Results

To test the simulation workflow previously described, we revisit the example
of ionic capacitor introduced in Ref. 15. In particular, we consider a Au(100)
electrode made of 4 unit cell repetitions along the $xy$ plane and 7 atomic
layers along $z$, which is put in contact with an aqueous solution of NaCl 5 M
under 3D periodic boundary conditions (Fig. 1). The electrolyte is made of 200
water molecules and 20 NaCl ions pairs, which, together with the gold
electrode, occupy a box size of $L_{x/y}=11.54$ Å and $L_{z}=64.34$ Å.
Classical electrolyte interactions are defined using a TIP4P/2005 water model
Abascal and Vega (2005) with Na+ and Cl- Lennard-Jones parameters taken from
Ref. 35; fixed classical atomic charges are given by $q_{\text{M}}=-1.1128$,
$q_{\text{H}}=+0.5564$ and $q_{\text{Na/Cl}}=\pm 1.0$.

Following the discussion reported in Sec. II, we apply a uniform electric
field along the direction orthogonal to the Au(100) surface, allowing us to
mimic a pair of equivalent gold electrodes that are kept at a constant voltage
drop. In particular, we consider two different cell potentials, namely $\Delta
V=0$ V and $\Delta V=1$ V, which correspond to applied fields of
$\varepsilon_{z}=0.0$ V/Å and $\varepsilon_{z}=0.016$ V/Å, respectively.
Reference QM/MM calculations are performed for 2000 uncorrelated frames
selected from a classical MetalWalls trajectory obtained running a finite-
field simulation at a cell voltage of $1$ V. The basis set used for the gold
electron density representation is defined from an uncontracted version of the
automatically generated auxiliary functions implemented in CP2K. Stoychev,
Auer, and Neese (2017) In particular, we adopt an auxiliary basis set that
includes up to $l_{\text{max}}=4$ spherical harmonics, accounting for a total
of 182 Gaussian-type functions per gold atom.

Figure 1: Representation of an ionic capacitor made of a Au(100) gold
electrode put in contact with a water/NaCl electrolyte solution under 3D
periodic boundary conditions. According to Ref. 18, a finite electric field
$\varepsilon_{z}$ is applied along $z$ to mimic a voltage drop of $\Delta
V=-\varepsilon_{z}L_{z}$ between two parallel electrodes.

###  V.1 Validation of SALTED electrode model

From the QM/MM dataset so generated, SALTED models of the electrode charge
density are trained following the discussion of Sec. II-B. In practice, a
subset of $M=400$ gold atomic environments is selected to recast the learning
problem into a low dimensional space. Grisafi  _et al._ (2023b) Symmetry-
adapted kernel functions are computed from LODE descriptors that include both
atom-density and potential-like structural features. Grisafi  _et al._
(2023a); Grisafi, Nigam, and Ceriotti (2021) In this work, these are defined
using Gaussian widths of $\sigma=0.5$ Å and $\sigma=2$ Å, respectively,
together with a radial cutoff of $r_{\text{cut}}=10$ Å. To validate the model,
we select 1600 random configurations for training and retain the remaining 400
for testing. The prediction error is measured consistently with the definition
of the SALTED loss function reported in Eq. (4). In particular, we compute the
percentage root mean square error as

$\%\text{RMSE}=\sqrt{\frac{\mathcal{L}_{\text{test}}(\boldsymbol{w})}{\sigma^{2}_{\text{test}}}}\times
100\,,$ (22)

where $\mathcal{L}_{\text{test}}(\boldsymbol{w})$ follows from Eq. (4), while

$\sigma^{2}_{\text{test}}=\sum_{I=1}^{N_{\text{test}}}\left[\boldsymbol{c}^{\text{DF}}_{I}\right]^{T}\mathbf{J}_{I}\,\boldsymbol{c}^{\text{DF}}_{I}\,,$
(23)

represents the variance of the learning target computed over the test set.
Learning curves are reported in the Supplementary Material.

Figure 2: Electrode charge density response (a) and derived electric field (b)
induced by the classical charges of a test electrolyte configuration. Red
lines: SALTED prediction results. Black dashed lines: Kohn-Sham DFT reference.
Shaded gray area indicates the region occupied by the gold electrode.

With as few as 200 training structures, the SALTED error saturates at $\sim$3%
RMSE, demonstrating the capability of the model of capturing the non-local
fluctuations of the electronic charge induced by different spatial
arrangements of the electrolyte atoms. This result is confirmed by the
excellent agreement between the Kohn-Sham and predicted electron density
variation $\Delta n_{e}$, as depicted in Fig. 2-a for a representative test
configuration. Despite the slightly less accurate prediction of the charge
density profile with respect to what reported in Ref. 15, we find that the use
of a Coulomb metric in the SALTED loss function allows us to obtain a highly
accurate electric field. This is illustrated in Fig. 2-b, where the
longitudinal variation of the predicted electric field with respect to the
field of the isolated electrode, $\Delta E_{z}$, is reported against the
quantum-mechanical reference. Notably, we observe that the agreement with DFT
is especially good outside the electrode region, being the most relevant for
driving the dynamics of the electrolyte.

As an ultimate proof of the model accuracy, we compute the predicted Cartesian
components of the electrostatic atomic forces associated with the 400
electrolyte configurations used for testing, as compared against the quantum-
mechanical reference. Both reference and predicted forces are computed as in
Eq. (20), starting from the DFT and SALTED electron-density coefficients,
respectively. Consistently with the high level of statistical correlation
shown in Fig. 3, we obtain a remarkably small root mean square error of $1.0$
meV/Å, which amounts to about 0.7% of the standard deviation of the DFT forces
in the test set. By and large, these results highlight the intrinsic
effectiveness of a machine learning model that has the electron density as a
prediction target to compute the electrostatic properties of the system.

Figure 3: Correlation plot between the predicted and reference electrostatic
forces generated by the electric field of a gold electrode, as computed over
256’000 electrolyte atoms belonging to 400 test configurations. Red dots:
SALTED prediction. Black dashed line: KS-DFT reference. Predicted forces are
associated with a $\text{RMSE}=1.0$ meV/Å over a standard deviation of
$\sigma_{f}=141.9$ meV/Å.

###  V.2 Simulation of ionic capacitor

We now report the results of the data-driven molecular dynamics, which is run
following the simulation workflow described in Sec. IV. Finite-field
simulations associated with applied cell potentials of $\Delta V=0$ V and $-$1
V are performed at room temperature in timesteps of 2 fs using a Nosé-Hoover
thermostat. Nosé (1984) Predictions of the electrostatic atomic forces are
obtained from a SALTED model trained on $N=2000$ configurations, as validated
from the discussion previously carried out. By parallelizing the calculation
over an Intel i9-12900 CPU, we find that each step of the dynamics takes about
0.5 s, corresponding to a speedup $>10^{3}$ with respect to first-principles
QM/MM calculations. This substantial increase in computational efficiency
results to be decisive for achieving timescales of the order of nanoseconds,
thus making it possible to access the slow dynamics of the electrolyte while
preserving a highly accurate representation of electrode charge density. For
both applied voltages, in particular, we are able to reach long simulation
times of $\approx 3$ ns.

It is instructive to compare our results against those of a classical
MetalWalls simulation. In this case, the electronic polarization of the gold
electrode is approximated through a charge-equilibration scheme that makes use
of fluctuating Gaussian charges. Marin-Laflèche  _et al._ (2020) Specifically,
we represent the partial charges on the gold atoms through a Gaussian width of
$\sigma_{\text{MW}}=1.06$ Å, and simulate classical trajectories for a total
time that is, once again, of about 3 ns. We note that similar simulations are
only 50 times faster than those performed through our method. This is
remarkable if considering that, in addition to the SALTED predictions, quantum
electrostatic forces are computed from the charge density generated by 182
spherical harmonics functions per gold atom, rather than from a single
isotropic Gaussian function.

Figure 4: Longitudinal density profiles of oxygen atoms on the right hand side
of the ionic capacitor studied in this work, subject to a null cell potential
of $\Delta V=0$ V. Dashed line: classical MetalWalls simulation. Full line:
SALTED/MetalWalls simulation. Shaded gray area indicates the gold electrode.

Regardless from the applied potential, we find that the thermal distribution
of water molecules obtained through our method displays a systematic shift
with respect to the result of a classical simulation. In particular, we
observe a distancing of about 0.5 Å of the first adsorption peak, which comes
along with overall smoother density profiles. This is shown in Fig 4 in the
case of a vanishing applied voltage, indicating that a classical electrode
model tends to give a stiffer potential of mean force with respect to what
expected from a first-principles approach.

Figure 5: Longitudinal density profiles of sodium chloride (a), and hydrogen
atoms (b) at the negatively charged electrode surface of the ionic capacitor
studied in this work, subject to an applied cell potential of $\Delta V=1$ V.
Dashed lines: classical MetalWalls simulation. Full lines: SALTED/MetalWalls
simulation. Red and blue color codes refer to Na+ and Cl- ions, respectively.
Shaded blue area indicates the gold electrode.

More striking differences between the classical and data-driven simulations
emerge when looking at the results obtained under an applied potential $\Delta
V=1$ V. As an example, we report in Fig. 5-a the ionic density profiles of Na+
and Cl- associated with the screening of the negatively charged metal surface.
While the distribution of chlorine shows an overall agreement between the two
methods, a qualitative difference is found in the adsorption of sodium
cations. In particular, a double pick is produced as a result of the classical
simulation, which does not appear when using SALTED to predict the quantum-
mechanical response of the electrode. This artifact of the classical model can
be associated with the rigidity by which the distribution of the metal charge
is represented, which mostly affects the accuracy of the electric field in the
proximity of the gold electrode. In fact, it has already been shown that,
depending on the chosen value of $\sigma_{\text{MW}}$, MetalWalls simulations
can yield substantial variations in the distribution of Na+ at the Au(100)
surface. Serva  _et al._ (2021)

In addition to the ionic profiles, a clear-cut discrepancy is also observed in
the distribution of hydrogen atoms; this is illustrated in Fig. 5-b.
Interestingly, we find that a prepeak is produced by the classical simulation
at about 2 Å from the gold slab. This prepeak is associated with an occasional
reorientation of the water molecules that expose the hydrogen atoms towards
the electrode, a phenomenon that is however not predicted by our method. This
is confirmed by an analysis of the orientation of the water dipoles with
respect to the normal to the electrode surface, as reported in the
Supplementary Material.

A further benefit of driving the dynamics of the electrochemical interface
through electron-density predictions is that of giving access to the total
charge $\pm Q$ accumulated at the electrode surfaces. In a finite-field
simulation setup, $Q$ can be computed by integrating $\Delta n_{e}$ from the
middle of the metal slab to the classical electrolyte region. Grisafi  _et
al._ (2023a) In Fig. 6, we report the time evolution of $Q$, computed every
picosecond, for both null and finite cell potentials. When considering the
results at $\Delta V=0$ V, we observe a fluctuation of the electrode charge
around 0 for both classical and SALTED-based simulations. This is to be
expected from the symmetry of the isolated metallic slab, implying that no net
polarization can be found at thermal equilibrium without an externally applied
electric bias. Conversely, simulations run at $\Delta V=1$ V lead to a finite
average charge $\langle Q\rangle$. In particular, we find that the computed
classical value of $\langle Q\rangle$ is about twice as large than what
predicted through our method, i.e., 1.16$e$ versus 0.63$e$, respectively. This
difference in the system’s electric capacitance is due to a substantial
overestimation of the electronic polarization of the metallic slab obtained
through the classical model, a result that was already hinted in Ref. 15 when
comparing predictions of $\Delta n_{e}$ performed on a test classical
trajectory. Although a better agreement might be obtained by tuning the value
of the classical Gaussian width $\sigma_{\text{MW}}$, the observed discrepancy
suggests that disposing of a physically sound electrode model is preferable to
obtain reliable results.

Figure 6: Evolution of the electrode surface charge in the ionic capacitor
studied in this work. Black and blue lines refer to SALTED/MetalWalls
simulations at 0 and $1$ V, respectively. Gray and light blue lines refer to
classical MetalWalls simulations at 0 and $1$ V, respectively.

In addition to $\langle Q\rangle$, having access to the thermal fluctuations
of $Q$ allows us to provide an estimate of the differential capacitance
associated with the formation of the electrical double layer: Scalfi  _et al._
(2020b)

$C^{\text{EDL}}_{\text{diff}}=\beta\left<(Q-\langle Q\rangle)^{2}\right>\,.$
(24)

At $\Delta V=0$ V, a similar differential capacitance is obtained for both
type of simulations, resulting in a value of
$C^{\text{EDL}}_{\text{diff}}=8.9\mu$F/cm2 for the classical model, and a
slightly reduced value of $C^{\text{EDL}}_{\text{diff}}=7.4\mu$F/cm2 for the
data-driven model. This discrepancy is magnified when looking at the results
at $\Delta V=1$ V. In line with what already computed in Ref. 15, the
classical simulation is associated with a value of
$C^{\text{EDL}}_{\text{diff}}=10.1\mu$F/cm2. Conversely, our SALTED/MetalWalls
simulation comes along with much more attenuated charge fluctuations, yielding
a differential capacitance that is about twice as small, i.e.,
$C^{\text{EDL}}_{\text{diff}}=4.6\mu$F/cm2. While simulations longer than 3 ns
would be required to provide a tighter convergence of
$C^{\text{EDL}}_{\text{diff}}$, these results remark, once again, the need of
introducing an accurate representation of the electrode in the description of
the electrochemical interface.

##  VI Conclusions

The presented method represents a relevant example of how to incorporate a
machine-learning approach that addresses the prediction of electronic-
structure properties into the atomistic study of complex materials. By
providing a formally rigorous integration between the SALTED and MetalWalls
programs, we have shown how to extend the reach of applications of QM/MM
approaches that aim at preserving a quantum description of metal electrodes in
the simulation of electrochemical interfaces. Most notably, the possibility of
driving the dynamics of the electrolyte over nanoseconds timescales represents
a great leap forward towards obtaining a thermally relaxed, yet accurate,
description of the electrical double layer.

By taking an ionic capacitor as an example, we have found that disposing of a
quantum treatment of the electrode charge density represents a critical aspect
for predicting the structural properties of the electrolyte at the interface.
In particular, our results remark how a classical electrode model cannot
reproduce the complex nature of the system’s electrostatics close to the metal
surface. This is not surprising considering that strong electric fields are
known to be encountered whenever an electrolyte atom penetrates the electronic
cloud of the metal electrode, Schmickler (1996) an effect that is especially
difficult to include through the use of classical atomic charges. In this
regard, an interesting question will be that of understanding how our method
compares against data-driven charge equilibration schemes based on partial
charges, where a ML representation of atomic electronegativies, Ko  _et al._
(2021) electric dipoles, Staacke  _et al._ (2022) and/or polarizabilities,
Shao  _et al._ (2022) is adopted to predict the electrostatics of the system.
We mention, in particular, a recent ML interface with MetalWalls that can be
used to refine the atomic Gaussian charges of the electrode by relying on a
linear density-response formalism. Dufils  _et al._ (2023)

In perspective, a notable extension of our method could involve promoting
redox active molecular species to the QM region, thus enabling the simulation
of electrochemical reactions in the presence of an explicit ionic solution. In
this context, an important aspect will be that of integrating our approach
with already existing local ML potentials. This possibility would be
especially attractive both for studying the charging mechanisms in
pseudocapacitive materials, Choi  _et al._ (2020) as well as for rationalizing
the role of the electrolyte in electrocatalytic processes. Ringe  _et al._
(2019) Furthermore, disposing of a three-dimensional map of the electronic
charge density will constitute an added value of our method in view of
predicting the regioselectivity of the catalyst to electron-transfer
phenomena.

## Supplementary Material

See the Supplementary Material for a complete derivation of the electronic
electric field, together with complementary results.

## Author Declarations

The authors have no conflicts to disclose.

## Source code

Printing of density coefficients and 2-center auxiliary integrals is made
available under the official trunk version of CP2K.
111https://github.com/cp2k/cp2k The latest SALTED implementation used to
produce the results of this work can be accessed from GitHub.
222https://github.com/andreagrisafi/SALTED A version of MetalWalls suitable to
be interfaced with SALTED can be downloaded from GitLab.
333https://gitlab.com/andreagrisafi/metalwalls/-/tree/salted_interface?ref_type=heads

## Data availability

The data that support the findings of this study are openly available from
Zenodo. 444https://zenodo.org/doi/10.5281/zenodo.11175494

###### Acknowledgements.

We thank Alessandra Serva, Laura Scalfi and Augustin Bussy for useful
discussion. This work was supported by the French National Research Agency
under the France 2030 program (Grant ANR-22-PEBA-0002).

## References

  * Dawlaty  _et al._ (2023) J. M. Dawlaty, S. Perkin, M. Salanne, and A. P. Willard, The Journal of Chemical Physics 159, 150401 (2023).
  * Schwarz and Sundararaman (2020) K. Schwarz and R. Sundararaman, Surface Science Reports 75, 100492 (2020).
  * Groß (2023) A. Groß, Current Opinion in Electrochemistry 40, 101345 (2023).
  * Goloviznina  _et al._ (2019) K. Goloviznina, J. N. Canongia Lopes, M. Costa Gomes, and A. A. H. Pádua, Journal of Chemical Theory and Computation 15, 5858 (2019).
  * Li  _et al._ (2019) C.-Y. Li, J.-B. Le, Y.-H. Wang, S. Chen, Z.-L. Yang, J.-F. Li, J. Cheng, and Z.-Q. Tian, Nature Materials 18, 697 (2019).
  * Khatib  _et al._ (2021) R. Khatib, A. Kumar, S. Sanvito, M. Sulpizi, and C. S. Cucinotta, Electrochimica Acta 391, 138875 (2021).
  * Chen  _et al._ (2022) A. Chen, J.-B. Le, Y. Kuang, and J. Cheng, The Journal of Chemical Physics 157, 094702 (2022).
  * Jeanmairet, Rotenberg, and Salanne (2022) G. Jeanmairet, B. Rotenberg, and M. Salanne, Chemical Reviews 122, 10860 (2022).
  * Scalfi  _et al._ (2020a) L. Scalfi, T. Dufils, K. G. Reeves, B. Rotenberg, and M. Salanne, The Journal of Chemical Physics 153, 174704 (2020a).
  * Goloviznina  _et al._ (2024) K. Goloviznina, J. Fleischhaker, T. Binninger, B. Rotenberg, H. Ers, V. Ivanistsev, R. Meissner, A. Serva, and M. Salanne, arXiv , 204.07848 (2024).
  * Elliott  _et al._ (2023) J. D. Elliott, M. Chiricotto, A. Troisi, and P. Carbone, Carbon 207, 292 (2023).
  * Wei  _et al._ (2024) Z. Wei, J. D. Elliott, A. A. Papaderakis, R. A. Dryfe, and P. Carbone, Journal of the American Chemical Society 146, 760 (2024).
  * Abidi and Steinmann (2023) N. Abidi and S. N. Steinmann, ACS Appl. Mater. Interfaces 15, 25009 (2023).
  * Takahashi, Nakano, and Sato (2022) K. Takahashi, H. Nakano, and H. Sato, The Journal of Chemical Physics 157, 234107 (2022).
  * Grisafi  _et al._ (2023a) A. Grisafi, A. Bussy, M. Salanne, and R. Vuilleumier, Phys. Rev. Mater. 7, 125403 (2023a).
  * Grisafi and Lewis (2023) A. Grisafi and A. Lewis, “Salted: Symmetry-adapted learning of three-dimensional electron densities,”  (2023).
  * Coretti  _et al._ (2022) A. Coretti, C. Bacon, R. Berthin, A. Serva, L. Scalfi, I. Chubak, K. Goloviznina, M. Haefele, A. Marin-Laflèche, B. Rotenberg, S. Bonella, and M. Salanne, J. Chem. Phys. 157, 184801 (2022).
  * Dufils  _et al._ (2019) T. Dufils, G. Jeanmairet, B. Rotenberg, M. Sprik, and M. Salanne, Phys. Rev. Lett. 123, 195501 (2019).
  * Grisafi  _et al._ (2023b) A. Grisafi, A. M. Lewis, M. Rossi, and M. Ceriotti, Journal of Chemical Theory and Computation 19, 4451 (2023b).
  * Lewis  _et al._ (2021) A. M. Lewis, A. Grisafi, M. Ceriotti, and M. Rossi, Journal of Chemical Theory and Computation 17, 7203 (2021).
  * Kühne  _et al._ (2020) T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hutter, The Journal of Chemical Physics 152, 194103 (2020).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * VandeVondele and Hutter (2007) J. VandeVondele and J. Hutter, The Journal of Chemical Physics 127, 114105 (2007).
  * Goedecker, Teter, and Hutter (1996) S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 1703 (1996).
  * Guidon, Hutter, and VandeVondele (2009) M. Guidon, J. Hutter, and J. VandeVondele, Journal of Chemical Theory and Computation 5, 3010 (2009), pMID: 26609981.
  * Bussy, Schütt, and Hutter (2023) A. Bussy, O. Schütt, and J. Hutter, The Journal of Chemical Physics 158, 164109 (2023).
  * Briling, Fabrizio, and Corminboeuf (2021) K. R. Briling, A. Fabrizio, and C. Corminboeuf, The Journal of Chemical Physics 155, 024107 (2021).
  * Grisafi and Ceriotti (2019) A. Grisafi and M. Ceriotti, J. Chem. Phys. 151, 204105 (2019).
  * Grisafi, Nigam, and Ceriotti (2021) A. Grisafi, J. Nigam, and M. Ceriotti, Chem. Sci. 12, 2078 (2021).
  * Huguenin-Dumittan  _et al._ (2023) K. K. Huguenin-Dumittan, P. Loche, N. Haoran, and M. Ceriotti, The Journal of Physical Chemistry Letters 14, 9612 (2023).
  * (31) G. Fraux, P. Loche, S. Kliavinek, K. K. Huguenin-Dumittan, D. Tisi, and A. Goscinski, “rascaline,” .
  * Marin-Laflèche  _et al._ (2020) A. Marin-Laflèche, M. Haefele, L. Scalfi, A. Coretti, T. Dufils, G. Jeanmairet, S. K. Reed, A. Serva, R. Berthin, C. Bacon, S. Bonella, B. Rotenberg, P. A. Madden, and M. Salanne, Journal of Open Source Software 5, 2373 (2020).
  * Marx and Hutter (2009) D. Marx and J. Hutter,  _Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods_ (Cambridge University Press, 2009).
  * Abascal and Vega (2005) J. L. F. Abascal and C. Vega, The Journal of Chemical Physics 123, 234505 (2005).
  * Yagasaki, Matsumoto, and Tanaka (2020) T. Yagasaki, M. Matsumoto, and H. Tanaka, Journal of Chemical Theory and Computation 16, 2460 (2020).
  * Stoychev, Auer, and Neese (2017) G. L. Stoychev, A. A. Auer, and F. Neese, Journal of Chemical Theory and Computation 13, 554 (2017).
  * Nosé (1984) S. Nosé, J. Chem. Phys. 81, 511 (1984).
  * Serva  _et al._ (2021) A. Serva, L. Scalfi, B. Rotenberg, and M. Salanne, The Journal of Chemical Physics 155, 044703 (2021).
  * Scalfi  _et al._ (2020b) L. Scalfi, D. T. Limmer, A. Coretti, S. Bonella, P. A. Madden, M. Salanne, and B. Rotenberg, Phys. Chem. Chem. Phys. 22, 10480 (2020b).
  * Schmickler (1996) W. Schmickler, Chemical Reviews 96, 3177 (1996).
  * Ko  _et al._ (2021) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, Nature Communications 12, 398 (2021).
  * Staacke  _et al._ (2022) C. G. Staacke, S. Wengert, C. Kunkel, G. Csányi, K. Reuter, and J. T. Margraf, Machine Learning: Science and Technology 3, 015032 (2022).
  * Shao  _et al._ (2022) Y. Shao, L. Andersson, L. Knijff, and C. Zhang, Electronic Structure 4, 014012 (2022).
  * Dufils  _et al._ (2023) T. Dufils, L. Knijff, Y. Shao, and C. Zhang, Journal of Chemical Theory and Computation 19, 5199 (2023).
  * Choi  _et al._ (2020) C. Choi, D. S. Ashby, D. M. Butts, R. H. DeBlock, Q. Wei, J. Lau, and B. Dunn, Nature Reviews Materials 5, 5 (2020).
  * Ringe  _et al._ (2019) S. Ringe, E. L. Clark, J. Resasco, A. Walton, B. Seger, A. T. Bell, and K. Chan, Energy Environ. Sci. 12, 3001 (2019).
  * Note (1) https://github.com/cp2k/cp2k.
  * Note (2) https://github.com/andreagrisafi/SALTED.
  * Note (3) https://gitlab.com/andreagrisafi/metalwalls/-/tree/salted_interface?ref_type=heads.
  * Note (4) %****␣text.tex␣Line␣1125␣****https://zenodo.org/doi/10.5281/zenodo.11175494.


# Gnuastro: measuring radial profiles from images

Raúl Infante-Sainz Mohammad Akhlaghi Sepideh Eskandarlou Centro de Estudios de
Física del Cosmos de Aragón (CEFCA), Plaza San Juan, 1, E-44001, Teruel, Spain
Raúl Infante-Sainz<EMAIL_ADDRESS>
###### Abstract

Radial profiles play a crucial role in the analysis and interpretation of
astronomical data, facilitating the extraction of spatial information.
However, highly customizable (for different scenarios) measurements over each
elliptical annulus can be challenging. In response, we present astscript-
radial-profile, which is part of Gnuastro from version 0.15 and has an
extensive documentation. A convenient feature of this program is its
capability to make the measurements with different operators (mean, median,
sigma-clipping, and many more) over ellipses, very quickly and directly on the
command-line with minimal dependencies. This research note is reproducible
with Maneage, on the Git commit 104aad5.

Astronomy software (1855), Astronomical techniques (1684), Astronomy data
visualization (1968), Broad band photometry (184), Open source software (1866)

##  1 Introduction

A radial profile consists of measurements over elliptical annuli at increasing
distances/radii. In general, a radial profile can be generalized to any number
of dimensions. For example, in a 3D dataset (cube), a radial profile consists
of measurements over growing ellipsoidal shells. Radial profiles are most
commonly employed in 2D images because the third dimension in 3D astronomical
data (IFUs or radio data cubes) is spectral, not spatial. The radial profile
is a table with two columns: the distance to the center and the measured
quantity.

Even with good visualization tools and proper color bars set, quantifying the
spatial variation of an object’s brightness, and comparing it with other
objects is challenging. A radial profile is therefore a convenient way of
extracting spatial information of an astronomical target.

Astronomical sources can have irregular shapes (for example stellar streams
and Galactic cirrus). In these cases, it is possible to compute radial
profiles for a specified azimuth angle or direction. More generally,
astronomical sources, driven by gravitational forces, tend to exhibit
symmetrical, often ellipse-symmetric shapes. As a consequence, this type of
aperture is very useful to derive physical properties such as surface
brightness profiles, mass density profiles, and star formation rate profiles,
among others, as demonstrated in recent studies (see e.g., Infante-Sainz et
al., 2020; Liu et al., 2022; Mosenkov et al., 2020; Martínez-Lombilla et al.,
2023; Watkins et al., 2023).

Once the annuli are defined, selecting the appropriate operator for measuring
the properties becomes crucial. While simple mean or sum measurements are
common, other situations require more sophisticated operators. For instance,
measuring the brightness within each radius with a simple mean or sum can be
severely affected by outliers at certain radii: foreground or background
objects or star-forming regions within the galaxy itself (see Figure 1). As a
consequence, robust estimators of the measurement such as median or sigma
clipping operators are desired, but a diverse set of measurements are rarely
implemented in many implementations of radial profiles.

Figure 1:  Color image and surface brightness radial profiles of M81 galaxy.
_Left:_ Color image with gray background obtained by using the astscript-
color-faint-gray script (Infante-Sainz & Akhlaghi, 2024). _Right:_ Surface
brightness radial profiles of the rSDSS image, measured with two different
operators. The y-axis represents the surface brightness (SB), while the x-axis
represents the distance to the center of the galaxy. The red circles
correspond to the measurements obtained by a simple mean, and they exhibit
peaks due to the bright objects (stars). In contrast, the blue squares
represent the measurements obtained using a $3\sigma$ clipped mean, resulting
in a smoother profile by removing the extreme outlier values (for example
those from contaminant sources).

In this research note, we present a powerful solution to address these
problems: the Gnuastro executable program astscript-radial-profile. This
program addresses the complexities arising from the generation of radial
profiles on astronomical images. It is able to create radial profiles of
astronomical sources very fast, from the command-line, and with an enormous
variety of measurements111Measurements:
https://www.gnu.org/software/gnuastro//manual/html_node/MakeCatalog-
measurements.html. The following section provides an overview, using the
radial profile of the M81 galaxy as an illustrative example.

##  2 Creating radial profiles

Gnuastro’s astscript-radial-profile is an installed script (written in POSIX
shell), and it uses several Gnuastro programs (that are installed with it and
available on the host operating system) to generate the radial profile. It is
possible to run the script and keep the intermediate files (passed between the
programs) to inspect and optimize the execution for each particular science
case. A summary of the general steps are described below.

  * •

Aperture construction. The script uses Gnuastro’s MakeProfiles (astmkprof) to
construct elliptical annuli based on user-defined parameters, including the
center position, axis ratio, position angle, azimuthal angle, and others.

  * •

Measurement extraction. Gnuastro’s MakeCatalog (astmkcatalog) is then used to
obtain the requested measurements over the elliptical annuli generated
previously.

  * •

Output generation. The script produces a formatted table (possibly FITS)
containing the resulting radial profile.

By default, if no options are specified, the program computes the radial
profile using circular apertures centered on the input image center, employing
a simple mean as the measurement operator. To customize the radial profile,
users can modify multiple options and parameters, allowing adjustments to the
ellipse positioning, shape, and measurement operators. Beyond the large number
of possible measurements, we find particularly interesting the possibility of
changing the size of the apertures by undersampling or oversampling, as well
as the generation of radial profiles with sub-pixel precision. For a
comprehensive understanding of the different parameters and how to use this
script in more detail, we refer to the Gnuastro manual222Manual:
https://www.gnu.org/software/gnuastro/manual.

The left panel of Figure 1 depicts a color (and gray background) visualization
of the M81 galaxy captured from the J-PLUS DR3 (Cenarro et al., 2019) using
the iSDSS, rSDSS, and gSDSS filters. This color image has been obtained with
another Gnuastro’s installed scripts astscript-color-faint-gray script (for
more details, see Infante-Sainz & Akhlaghi, 2024). In this image there are
several other sources and over M81, star-forming regions produce non
elliptically-symmetric features. The right panel of Figure 1 shows the radial
profiles computed by using two different operators: simple mean and sigma-
clipped mean. The radial profiles look quite different. While the simple mean
radial profile has several peaks that are due to the presence of nearby
objects (bright stars in the surroundings of the M81 galaxy), the sigma
clipping radial profile is smoother and does not have these peaks.

##  3 Acknowledgments

The workflow uses Maneage (_Man_ aging data lin _eage_ , Akhlaghi et al.,
2021, commit 8161194). This research note is created from the Git commit
104aad5, hosted on Codeberg333Git repository of the paper (radial-profile
branch): https://codeberg.org/gnuastro/papers. which is archived on Software
Heritage444Software Heritage (SoftWare Hash IDentifier, SWHID):
swh:1:dir:d5029e066916cb64f0d95d20eb88294acc78b2b1 for longevity. Supplements
are also available on Zenodo555Zenodo:
https://doi.org/10.5281/zenodo.10124582.

The analysis uses GNU Astronomy Utilities (Gnuastro, Akhlaghi & Ichikawa,
2015; Akhlaghi, 2019) v0.21.43-3101. Gnuastro has been funded by the Japanese
MEXT scholarship and its Grant-in-Aid for Scientific Research (21244012,
24253003), ERC 339659-MUSICOS, Spanish MINECO AYA2016-76219-P, and
NextGenerationEU ICTS-MRR-2021-03-CEFCA. We acknowledge the funding by
Governments of Spain and Aragón through FITE and Science Ministry
(PGC2018-097585-B-C21, PID2021-124918NA-C43).

## References

  * Akhlaghi (2019) Akhlaghi, M. 2019, arXiv e-prints, arXiv:1909.11230, doi: 10.48550/arXiv.1909.11230
  * Akhlaghi & Ichikawa (2015) Akhlaghi, M., & Ichikawa, T. 2015, ApJS, 220, 1, doi: 10.1088/0067-0049/220/1/1
  * Akhlaghi et al. (2021) Akhlaghi, M., Infante-Sainz, R., Roukema, B. F., et al. 2021, CiSE, 23, 82, doi: 10.1109/MCSE.2021.3072860
  * Cenarro et al. (2019) Cenarro, A. J., Moles, M., Cristóbal-Hornillos, D., et al. 2019, A&A, 622, A176, doi: 10.1051/0004-6361/201833036
  * Infante-Sainz & Akhlaghi (2024) Infante-Sainz, R., & Akhlaghi, M. 2024, Research Notes of the AAS, 8, 10, doi: 10.3847/2515-5172/ad1aae
  * Infante-Sainz et al. (2020) Infante-Sainz, R., Trujillo, I., & Román, J. 2020, MNRAS, 491, 5317, doi: 10.1093/mnras/stz3111
  * Liu et al. (2022) Liu, Q., Abraham, R., Gilhuly, C., et al. 2022, ApJ, 925, 219, doi: 10.3847/1538-4357/ac32c6
  * Martínez-Lombilla et al. (2023) Martínez-Lombilla, C., Infante-Sainz, R., Jiménez-Ibarra, F., et al. 2023, A&A, 678, A62, doi: 10.1051/0004-6361/202346280
  * Mosenkov et al. (2020) Mosenkov, A. V., Smirnov, A. A., Sil’chenko, O. K., et al. 2020, MNRAS, 497, 2039, doi: 10.1093/mnras/staa1885
  * Watkins et al. (2023) Watkins, A. E., Salo, H., Kaviraj, S., et al. 2023, MNRAS, 521, 2012, doi: 10.1093/mnras/stad654

##  Appendix A Software acknowledgement

This research was done with the following free software programs and
libraries: 1.23, Bzip2 1.0.8, C compiler (Apple clang version 13.1.6
(clang-1316.0.21.2.5)), CFITSIO 4.1.0, CMake 3.24.0, Dash 0.5.11-057cd65,
Discoteq flock 0.4.0, Expat 2.4.1, File 5.42, Fontconfig 2.14.0, FreeType
2.11.0, GNU AWK 5.1.1, GNU Astronomy Utilities 0.21.43-3101 (Akhlaghi &
Ichikawa, 2015; Akhlaghi, 2019), GNU Autoconf 2.71, GNU Automake 1.16.5, GNU
Bash 5.2-rc2, GNU Bison 3.8.2, GNU Coreutils 9.1, GNU Diffutils 3.8, GNU
Findutils 4.9.0, GNU Grep 3.7, GNU Gzip 1.12, GNU Libtool 2.4.7, GNU M4
1.4.19, GNU Make 4.3, GNU Multiple Precision Arithmetic Library 6.2.1, GNU
Multiple Precision Floating-Point Reliably 4.1.0, GNU NCURSES 6.3, GNU Nano
6.4, GNU Readline 8.2-rc2, GNU Scientific Library 2.7, GNU Sed 4.8, GNU Tar
1.34, GNU Texinfo 6.8, GNU Wget 1.21.2, GNU Which 2.21, GNU gettext 0.21, GNU
gperf 3.1, GNU libiconv 1.17, GNU libunistring 1.0, GPL Ghostscript 9.56.1,
Git 2.37.1, Help2man , Less 590, Libffi 3.4.2, Libgit2 1.3.0, Libidn 1.38,
Libjpeg 9e, Libpaper 1.1.28, Libpng 1.6.37, Libtiff 4.4.0, Libxml2 2.9.12,
Lzip 1.23, OpenSSL 3.0.5, Perl 5.36.0, Python 3.10.6, WCSLIB 7.11, X11 library
1.8, XCB-proto (Xorg) 1.15, XZ Utils 5.2.5, Zlib 1.2.11, cURL 7.84.0, libICE
1.0.10, libSM 1.2.3, libXau (Xorg) 1.0.9, libXdmcp (Xorg) 1.1.3, libXext
1.3.4, libXt 1.2.1, libpthread-stubs (Xorg) 0.4, libxcb (Xorg) 1.15, pkg-
config 0.29.2, podlators 4.14, util-Linux 2.38.1, util-macros (Xorg) 1.19.3,
xorgproto 2022.1 and xtrans (Xorg) 1.4.0. The LaTeX source of the paper was
compiled to make the PDF using the following packages: courier 61719
(revision), epsf 2.7.4, etoolbox 2.5k, helvetic 61719 (revision), lineno 5.3,
pgf 3.1.10, pgfplots 1.18.1, revtex4-1 4.1s, tex 3.141592653, textcase 1.04
and ulem 53365 (revision). We are very grateful to all their creators for
freely providing this necessary infrastructure. This research (and many other
projects) would not be possible without them.