%% Supplementary Information for: "MTP–LES: Long-Range Electrostatics and Born Effective Charges from Invariant Moment Tensor Potentials"
%% Compiled separately into SI.pdf.

\documentclass{article}
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage{johd}
\usepackage{float}
\usepackage{subcaption}
\usepackage{graphicx}
\usepackage[percent]{overpic}
\usepackage{amsmath}
\usepackage{xr}
\externaldocument{main}

\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thesection}{S\arabic{section}}

% Roadmap callout: rendered as a plain paragraph with a bold header.
\newenvironment{roadmap}{%
  \par\medskip\noindent\textbf{How the SI maps to the main text.}\par\nobreak
}{%
  \par\medskip%
}

\title{Supplementary Information for:\\
``{Long-Range Electrostatics and Born Effective Charges from Invariant Moment Tensor Potentials with Latent Ewald Summation (MTP-LES}''}

% \author{
% Abibat Adekoya-Olowofela$^{a,b}$,
% Sukriti Manna$^{a,b}$,
% Shola Folarin$^{a,b}$\\[0.2em]
% Aditya Koneru$^{a,b}$,
% Subramanian Sankaranarayanan$^{a,b,*}$\\[0.45em]
% {\fontsize{9}{9.5}\selectfont
% \begin{tabular}{@{}c@{}}
% $^{a}$Department of Mechanical and Industrial Engineering,
% University of Illinois Chicago, Chicago, IL 60607, USA\\
% $^{b}$Center for Nanoscale Materials,
% Argonne National Laboratory, Lemont, IL 60439, USA\\[0.2em]
% $^{*}$Corresponding author: \texttt{skrssank@uic.edu}
% \end{tabular}}
% }
% }
\date{}

\begin{document}
\maketitle

\noindent This document collects supporting material for the main paper. Sections are referenced as ``Supplementary Note~SX'' and figures, tables, and equations carry an ``S'' prefix throughout.

\begin{roadmap}
\begin{itemize}
\item Supplementary Note~S1 gives the full mathematical specification of the rotationally invariant MTP descriptors summarized in the main-text Methods.
\item Supplementary Note~S2 provides extended quantitative detail supporting the net-charge comparison in Fig.~1 of the main text.
\item Supplementary Note~S3 gives the full \textsc{vasp} input specification for the HfO$_2$ dataset.
\item Supplementary Note~S4 provides quantitative analysis of the liquid-water infrared spectrum.
\item Supplementary Note~S5 reports per-species BEC RMSE values for water and dipeptides.
\item Supplementary Note~S6 reports the dipeptide dipole-loss ablation.
\item Supplementary Note~S7 provides the raw timing data behind the runtime benchmark.
% \item Supplementary Note~S8 collects auxiliary figures referenced from the main text and SI.
\item Supplementary Note~S8 documents the periodic-image neighbour-list implementation detail required to reproduce periodic MTP--LES training.
\end{itemize}
\end{roadmap}

\tableofcontents

\section{Full specification of the invariant MTP descriptors}
\label{si:mtp_descriptors}

This note specifies the rotationally invariant moment-tensor descriptors used in
this work, following the Moment Tensor Potential construction of
Shapeev~\cite{shapeev2016,Novikov_2021}. We employ a compact set of invariant
contractions of moment tensors, given explicitly below together with its
body-order and level truncation; this set is sufficient to reach the energy and
force accuracy reported for the systems studied.


\subsection{Moment tensors}

For each atom $i$, moment tensors of rank $\nu \in \{0, 1, 2\}$ are constructed from the local neighborhood within a cutoff radius $r_\mathrm{cut}$:
\begin{align}
M^{(0)}_{\mu,i} &= \sum_{j \in \mathcal{N}(i)} R_\mu(r_{ij}, s_i, s_j), \\
M^{(1)}_{\mu,i,\alpha} &= \sum_{j \in \mathcal{N}(i)} R_\mu(r_{ij}, s_i, s_j)\, r_{ij,\alpha}, \\
M^{(2)}_{\mu,i,\alpha\beta} &= \sum_{j \in \mathcal{N}(i)} R_\mu(r_{ij}, s_i, s_j)\, r_{ij,\alpha}\, r_{ij,\beta},
\end{align}
where $r_{ij} = |\mathbf{r}_j - \mathbf{r}_i|$, $r_{ij,\alpha}$ is the $\alpha$-th Cartesian component of $\mathbf{r}_j - \mathbf{r}_i$, $s_i, s_j$ are species indices, and $\mathcal{N}(i)$ is the set of neighbours of $i$ within $r_\mathrm{cut}$.


\subsection{Chebyshev radial basis}

The radial functions are expanded using the Chebyshev polynomial basis employed
in the MLIP-2 framework~\cite{Novikov_2021}. For each pair of atomic species
$(s_i,s_j)$, the radial basis is written as

\begin{equation}
R_\mu(r,s_i,s_j)
=
\sum_{n=1}^{n_{\mathrm{rad}}}
c_{\mu,s_i,s_j,n}\,
T_n\!\left(\xi(r)\right)
f_{\mathrm{cut}}(r),
\end{equation}

where $c_{\mu,s_i,s_j,n}$ are learnable species-pair-dependent coefficients and
$T_n$ denotes the $n$th Chebyshev polynomial of the first kind, generated
recursively as

\begin{equation}
T_0=1,\qquad
T_1=\xi,\qquad
T_n=2\xi T_{n-1}-T_{n-2}.
\end{equation}

The scaled radial coordinate

\begin{equation}
\xi(r)
=
\frac{2r-r_{\mathrm{min}}-r_{\mathrm{cut}}}
{r_{\mathrm{cut}}-r_{\mathrm{min}}}
\end{equation}

maps the interval
$[r_{\mathrm{min}},r_{\mathrm{cut}}]$
onto
$[-1,1]$.
Each radial basis function is multiplied by the smooth cutoff envelope

\begin{equation}
f_{\mathrm{cut}}(r)=
\begin{cases}
\displaystyle
\left(
\frac{r_{\mathrm{cut}}-r}
{r_{\mathrm{cut}}-r_{\mathrm{min}}}
\right)^2,
&
r<r_{\mathrm{cut}},
\\[8pt]
0,
&
r\ge r_{\mathrm{cut}}.
\end{cases}
\end{equation}

This formulation is equivalent to the Chebyshev radial basis used in
MLIP-2, differing only by the constant normalization factor
$(r_{\mathrm{cut}}-r_{\mathrm{min}})^{-2}$,
which is absorbed into the learnable coefficients
$c_{\mu,s_i,s_j,n}$.


\subsection{Rotational invariants by tensor contraction}

Following the Moment Tensor Potential (MTP) formalism, each moment tensor is
assigned a level
\begin{equation}
\mathrm{lev}\!\left(M_\mu^{(\nu)}\right)=2\mu+4\nu,
\end{equation}
which determines the complexity of the corresponding basis function. The level
of a product of moment tensors is defined as the sum of the levels of the
individual tensors. Rotationally invariant basis functions are constructed by
contracting moment tensors whose total level does not exceed
$\ell_{\mathrm{max}}$ and whose tensor rank is limited by
$\nu_{\mathrm{max}}$. Throughout this work, we use
$\ell_{\mathrm{max}}=24$ and $\nu_{\mathrm{max}}=2$.

Rather than employing the complete MTP invariant basis, we retain a reduced
subset that provides sufficient expressiveness while keeping the descriptor
compact. In particular, invariants such as
$\mathrm{tr}\,M_\mu^{(2)}$ and higher-order contractions involving
$M^{(2)}$ are omitted. The rotational invariants used in this work are
\begin{align}
&\text{zeroth-order moment:}
&
B
&=
M^{(0)}_{\mu},
\\
&\text{scalar products:}
&
B
&=
M^{(0)}_{\mu_0}M^{(0)}_{\mu_1},
\qquad
\left(M^{(0)}_{\mu_0}\right)^2,
\\
&\text{vector contraction:}
&
B
&=
\sum_{\alpha}
M^{(1)}_{\mu_1,\alpha}
M^{(1)}_{\mu_2,\alpha},
\\
&\text{Frobenius contraction:}
&
B
&=
\sum_{\alpha\beta}
M^{(2)}_{\mu_1,\alpha\beta}
M^{(2)}_{\mu_2,\alpha\beta},
\\
&\text{scalar-weighted vector contraction:}
&
B
&=
M^{(0)}_{\mu_0}
\sum_{\alpha}
M^{(1)}_{\mu_1,\alpha}
M^{(1)}_{\mu_2,\alpha},
\\
&\text{vector--matrix--vector contraction:}
&
B
&=
\sum_{\alpha\beta}
M^{(1)}_{\mu_1,\alpha}
M^{(2)}_{\mu_2,\alpha\beta}
M^{(1)}_{\mu_3,\beta}.
\end{align}

The scalar-weighted vector-contraction family includes the special case
$\mu_1=\mu_2$, namely
\begin{equation}
M^{(0)}_{\mu_0}
\left|\mathbf{M}^{(1)}_{\mu_1}\right|^2.
\end{equation}
These contractions generate rotationally invariant combinations of the scalar,
vector, and second-order tensor moments while retaining information about the
local atomic environment.

Finally, a one-hot encoding of the central atomic species is concatenated with
the invariant vector to form the per-atom descriptor
\begin{equation}
\mathbf{d}_i \in \mathbb{R}^{n_{\mathrm{feat}}}.
\end{equation}
The resulting descriptor is used as the input to both the short-range energy
network and the independent latent-charge prediction network introduced below.
\section{Charge-head construction: supplementary detail}
\label{si:charge_head}

The main text adopts a decoupled and neutrality-enforced charge head over the
native LES construction (see Results and Fig.~1 of the main text). This note
provides supporting details for that choice on PbTiO$_3$: the two constructions
compared, the per-configuration net cell charge that distinguishes them, and the
energy, force, and BEC metrics that explain why the residual charge is invisible
to standard MLIP error checks.

\paragraph{The two constructions.} Both heads sit on the same short-range MTP
backbone and are trained on identical data with identical descriptors, optimizer,
and loss weights; they differ only in how the latent charges are produced. The
native head predicts the latent charges directly from the shared MTP
descriptors, as in existing LES integrations. The decoupled head uses an
independent charge branch with a per-species bias and projects the mean charge to
zero per configuration, enforcing $\sum_i q_i = 0$ exactly.

\paragraph{Net cell charge.} Across the 45 BEC configurations, the native head
develops a large residual cell charge, $|\sum_i q_i| \approx 12\,e$ per cell,
whereas the decoupled head is neutral to within $10^{-5}\,e$. The latent charges
in LES are fixed only by the requirement that they reproduce energies and forces,
which for a periodic system leaves the total cell charge unconstrained; the
reciprocal-space (tinfoil) treatment then absorbs the residual as a uniform
compensating background. The result is energetically admissible but inconsistent
with the charge-neutral physical system, which is why an explicit neutrality
constraint is required rather than optional.

\paragraph{Energies, forces, and BECs are not enough.} The residual charge leaves
no signature in the standard metrics. Both heads reach essentially identical
Born-effective-charge correlation against the PBE-DFPT references
($R = 0.995$ native, $R = 0.994$ decoupled), and neither energy nor force
accuracy is degraded by the residual charge. The non-neutrality is therefore
exposed only by direct evaluation of the latent charges, not by the parity
metrics that would normally be used to validate the model---which is the central
reason we report $\sum_i q_i$ explicitly and adopt the decoupled head as the
production model.



\section{Full VASP input specification for HfO$_2$}
\label{si:vasp_inputs}

\begin{table}[H]
\centering
\caption{VASP input parameters for the HfO$_2$ dataset.}
\label{tab:si_vasp}
\begin{tabular}{lll}
\hline
Parameter & Value & Description \\
\hline
\textsc{VASP} version & 6.4.3 & Code version \\
Functional & PBE & GGA exchange--correlation \\
PAW potentials & \texttt{PAW\_PBE Hf} (20 Jan 2003) & Released with VASP 6.4.3 \\
& \texttt{PAW\_PBE O} (08 Apr 2002) & \\
\texttt{ENCUT} & 600 eV & Plane-wave cutoff \\
\texttt{EDIFF} & $10^{-8}$ eV & Electronic convergence \\
\texttt{ALGO} & Normal & Blocked Davidson \\
\texttt{ISMEAR} & 0 & Gaussian smearing \\
\texttt{SIGMA} & 0.1 eV & Smearing width, extrapolated to 0 \\
\texttt{ISPIN} & 2 & Spin polarization on \\
\texttt{LASPH} & .TRUE. & Non-spherical PAW contributions \\
\texttt{ISYM} & 0 & Symmetry disabled \\
$k$-point density & 0.040 \AA$^{-1}$ per recip.\ vector & $3\times2\times2$ for 96-atom cell \\
$k$-point grid & $\Gamma$-centred Monkhorst--Pack & \\
\texttt{NSW} & 0 & Static evaluations \\
\texttt{IBRION} & $-1$ & No ionic relaxation \\
\texttt{ISIF} & 2 & Positions \& cell fixed \\
\texttt{LWAVE} & .FALSE. & Wavefunctions not retained \\
\texttt{LCHARG} & .FALSE. & Charge density not retained \\
\hline
\end{tabular}
\end{table}

\section{Quantitative IR spectrum analysis}
\label{si:ir_quantitative}

Peak positions are determined from the smoothed spectrum maxima; integrated intensities are the area within $\pm 200$~cm$^{-1}$ of each peak after area-normalisation over 0--4000~cm$^{-1}$.

\begin{table}[H]
\centering
\caption{IR peak positions and band intensities for liquid water: MTP--LES vs experiment~\cite{Bertie:96}.}
\label{tab:si_ir}
\begin{tabular}{lcccc}
\hline
Band & MTP--LES position (cm$^{-1}$) & Exp.\ position (cm$^{-1}$) & MTP--LES intensity & Exp.\ intensity \\
\hline
Libration    & 601  & 676  & 0.196 & 0.173 \\
H--O--H bend & 1624 & 1648 & 0.080 & 0.046 \\
O--H stretch & 3437 & 3414 & 0.541 & 0.632 \\
\hline
\end{tabular}
\vspace{2pt}
{\small Intensities are band areas (dimensionless after area normalisation over 0--4000~cm$^{-1}$). Model spectrum from the dipole-current autocorrelation without quantum correction, consistent with the main-text figure.}
\end{table}
The librational position is within 11\%, the bending mode within 1.5\%, and the
O--H stretching position within 1\% of experiment. The lower-frequency band
intensities are overestimated and the O--H stretching intensity underestimated
relative to experiment, consistent with the absence of nuclear quantum effects
in classical MD~\cite{PhysRevLett.101.017801, Marsalek2017}.

\section{Per-species BEC RMSE}
\label{si:per_species_rmse}


\begin{table}[H]
\centering
\caption{Per-species diagonal and off-diagonal BEC RMSE for water and dipeptides, in elementary charge units.}
\label{tab:si_per_species_rmse}
\begin{tabular}{lcc}
\hline
System and species & Diagonal RMSE ($e$) & Off-diagonal RMSE ($e$) \\
\hline
Water --- O & 0.103 & 0.081 \\
Water --- H & 0.087 & 0.064 \\
\hline
Dipeptides --- H & 0.092 & 0.087 \\
Dipeptides --- C & 0.187 & 0.198 \\
Dipeptides --- N & 0.256 & 0.258 \\
Dipeptides --- O & 0.147 & 0.212 \\
\hline
\end{tabular}
\end{table}


\section{Dipeptide dipole-loss ablation}
\label{si:dipeptide_ablation}

\begin{figure}[H]
\centering
\includegraphics[width=1.0\textwidth]{images/dipep_validation_train_test_no_dipole.pdf}
\caption{Validation of MTP--LES on dipeptides \emph{without} dipole loss ($w_\mu = 0$). Panels and color coding follow main-text Fig.~7. The increased scatter in the dipole parity (\textbf{c}) illustrates the benefit
of the optional dipole loss for this finite-molecule dataset.}
\label{fig:si_dipep_no_dipole}
\end{figure}

\begin{table}[H]
\centering
\caption{Effect of the optional dipole loss on the dipeptide test set.}
\label{tab:si_dipeptide_ablation}
\begin{tabular}{lcc}
\hline
Metric & $w_\mu=0$ & $w_\mu=0.1$ \\
\hline
Dipole RMSE ($e$\AA{}) & 0.307 & 0.075 \\
Diagonal BEC $R^2$ & 0.899 & 0.931 \\
Diagonal BEC RMSE ($e$) & 0.187 & 0.154 \\
Off-diagonal BEC RMSE ($e$) & 0.181 & 0.164 \\
\hline
\end{tabular}
\end{table}


\section{Wall-time benchmark}
\label{si:timing_data}
\begin{table}[H]
\centering
\caption{Wall time per MD step for MTP-only and MTP--LES on liquid water, single A100-PCIE-40GB GPU, float32, mean over 1\,000 force evaluations.}
\label{tab:si_timing}
\begin{tabular}{cccc}
\hline
$N$ (atoms) & MTP-only (ms) & MTP--LES (ms) & Overhead (\%) \\
\hline
192    & 76  & 83   & 9.0 \\
384    & 104 & 108  & 3.8  \\
768    & 120 & 122  & 1.3  \\
1\,536  & 168 & 169  & 0.5  \\
2\,304  & 217 & 222  & 2.5  \\
3\,456  & 282 & 289  & 2.6  \\
5\,184  & 397 & 411  & 3.6  \\
6\,912  & 484 & 501  & 3.4  \\
9\,216  & 613 & 635  & 3.5  \\
12\,288 & 777 & 826  & 6.3  \\
15\,360 & 996 & 1\,059 & 6.3  \\
\hline
\end{tabular}
\end{table}

\section{Implementation note: neighbour-list offsets in periodic training}
\label{si:nl_offsets}

One implementation detail is important for reproducing periodic MTP--LES force
training. For periodic systems, the periodic-image offsets in the neighbour list
must be preserved as non-differentiable constants during backpropagation. If
these offsets are recomputed from positions at each step, the resulting forces
can become inconsistent between training and validation because the offsets carry
gradients that are not part of the intended force model. We therefore precompute
the periodic-image offsets from the reference positions and attach them to the
neighbour list as non-differentiable constants. With this treatment, training and validation forces
agree to within numerical precision for the periodic benchmarks reported in the
main text.




\bibliographystyle{unsrt}
\bibliography{bib}

\end{document}