# Learning intermolecular forces at liquid-vapor interfaces

Samuel P. Niblett  Department of Chemistry, University of California, Berkeley
CA 94609  Materials Science Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94609  Mirza Galib  Department of Chemistry, University of
California, Berkeley CA 94609  David T. Limmer<EMAIL_ADDRESS>Department of Chemistry, University of California, Berkeley CA 94609
Materials Science Division, Lawrence Berkeley National Laboratory, Berkeley,
CA 94609  Chemical Science Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94609  Kavli Energy NanoScience Institute, Berkeley, CA 94609

###### Abstract

By adopting a perspective informed by contemporary liquid state theory, we
consider how to train an artificial neural network potential to describe
inhomogeneous, disordered systems. Specifically, we find that neural network
potentials based on local representations of atomic environments are capable
of describing bulk liquid properties, but can fail at liquid-vapor interfaces.
This failure results from the unbalanced attractions that build up in the
presence of broken translation symmetry from long-ranged interactions, but
cancel in the translational invariant bulk. By incorporating explicit models
of the slowly-varying long-ranged interactions and training neural networks
only on the short ranged components, we can arrive at potentials that robustly
recover interfacial properties. We find that local neural network models can
sometimes approximate a local molecular field potential to correct for the
truncated interactions, but this behavior is variable and hard to learn.
Generally, we find that models with explicit electrostatics are easier to
train and have higher accuracy. We demonstrate this perspective in a simple
model of an asymmetric dipolar fluid where the exact long-ranged interaction
is known, and in an _ab initio_ water model where it is approximated.

## Introduction

Machine learning based forcefields have significantly broadened the scope of
molecular simulations.Behler (2016); Noé  _et al._ (2020); Bartók  _et al._
(2017) Supervised learning techniques have been used to refine the
parameterization of traditional physically motivated potentials,Cisneros  _et
al._ (2016); Li  _et al._ (2017); Unke and Meuwly (2019); Cubuk  _et al._
(2017) and enabled those employing artificial neural network (ANN)
representations.Zhang  _et al._ (2018a); Behler and Parrinello (2007); Behler
(2011); Wang  _et al._ (2018) ANN potentials in particular have demonstrated
their utility in the prediction of material properties,Csányi  _et al._
(2017); Marchand  _et al._ (2020); Bischak  _et al._ (2020) and in the
accurate calculation of thermodynamicBehler  _et al._ (2008); Cheng  _et al._
(2018); Morawietz  _et al._ (2016); Bonati and Parrinello (2018) and kinetic
properties.Yang  _et al._ (2021); Galib and Limmer (2021); Käser, Unke, and
Meuwly (2020); Liu  _et al._ (2018) While there have been significant advances
in the chemical complexity of systems studied with ANN potentials as they are
flexible enough describe reactivity,Behler (2017) the overwhelming majority of
studies have employed relatively simple bulk environments, with a few notable
exceptions.Natarajan and Behler (2016); Andrade  _et al._ (2020); Wohlfahrt,
Dellago, and Sega (2020); Galib and Limmer (2021); Yue  _et al._ (2021) A
barrier preventing the extension of ANN potentials into more complex
environments is that they are truncated at a finite interaction distance,
which foundational work in liquid state theory suggests can lead to
inaccuracies in systems lacking translational invariance.Weeks, Chandler, and
Andersen (1971); Lum, Chandler, and Weeks (1999); Weeks, Selinger, and
Broughton (1995) Here we employ a perspective based on the local molecular
field theory to rationalize the failings of truncated ANN potentials in
nonuniform systems and to develop a representation and training procedure that
admits their extension to liquid-vapor interfaces.

Local molecular field (LMF) theory was introduced by Weeks and collaborators
to understand the structure and thermodynamic properties of nonuniform,
disordered systems.Weeks, Selinger, and Broughton (1995); Weeks (2002) The LMF
theory uses the Yvonne-Born Green hierarchyHansen and McDonald (1990) to
demonstrate that a system with slowly varying long-ranged interactions, like
those associated with electrostatic or dispersion forces, can be mapped to a
system with short-ranged intermolecular interactions in an effective external
field. In a bulk environment, translational invariance dictates that the
effective external field is a constant and therefore inconsequential. At an
interface, this external field corrects for the unbalanced interactions in the
truncated effective model. This perspective has been used analytically in
theories of solvation, and computationally as a means to model electrostatic
interactions without traditional Ewald sums.Rodgers and Weeks (2008); Remsing,
Rodgers, and Weeks (2011); Rodgers  _et al._ (2006); Rodgers and Weeks (2008);
Remsing  _et al._ (2014); Remsing and Weeks (2013, 2016); Cox (2020)

ANN potentials are typically formulated with mathematical representations of a
system, so-called symmetry functions, that are short-ranged,Behler and
Parrinello (2007) yet aim to approximate intermolecular potentials that are
often long-ranged. As such, it is natural to use LMF theory as a lens to
understand the limitations and successes of truncated ANN potentials. This is
the perspective we adopt here. Consistent with LMF theory, we find that bulk
properties are insensitive to the description of long-ranged interactions,
while some interfacial properties, like the interfacial polarization, are
particularly sensitive. LMF theory also suggests a means of ameliorating these
mixed successes by incorporating the slowly-varying part of the electrostatic
interactions in extended ANN potentials.

The limitations of ANN potentials employing only short ranged interactions
have been discussed extensively, and a number of procedures have been proposed
to circumvent them.Grisafi and Ceriotti (2019); Deng  _et al._ (2019); Yue
_et al._ (2021); Bartók  _et al._ (2010); Artrith, Morawietz, and Behler
(2011); Yao  _et al._ (2018) The short-ranged character of ANN potentials
results from the symmetry functions using a cutoff distance, so the most
direct route to incorporate long-ranged interactions is to include non-local
features in those molecular representations.Grisafi and Ceriotti (2019)
However, a more common approach involves explicitly evaluating additional
long-ranged forces.Deng  _et al._ (2019); Yue  _et al._ (2021); Bartók  _et
al._ (2010); Artrith, Morawietz, and Behler (2011); Yao  _et al._ (2018) In
this framework, the neural network is trained to reproduce a truncated version
of the target potential, for example in a classical model the ANN might be
trained only on the van der Waals component of the energy. During simulation,
the interactions predicted by the ANN are combined with an explicit
calculation of the remaining long-ranged components. The advantage of this
method is that long-ranged interactions are included directly, while the
disadvantages are the cost and the ambiguity in selecting the models used to
evaluate the long-ranged forces. We adopt this latter approach, employing an
explicit description of the electrostatic interactions both in systems where
the true description is known and in cases where it is not and we must
approximate it. The addition of explicit long-ranged forces into the ANN
models increases the overall accuracy of the models, decreases the difficulty
training them, and consistently recovers the interfacial properties of the
reference model. Further, we find the ANNs are not very sensitive to the
specific auxiliary long range force model employed. In the following these
findings are explored within the context of a simple model of an asymmetric
dipolar fluid, and an _ab initio_ model of liquid water.

## Reference Models

In order to understand the fundamental limitations of truncated ANN potentials
and how to move beyond them, we study ANN potentials fitted to both a fixed
charge forcefield model and an _ab initio_ model. For the forcefield, we can
compute the reference values of quantities precisely in order to validate
collective properties of the ANN system. Further, for the forcefield the form
of the long-range force is known by construction so that we can assess its
impact in a controlled manner. We also study an _ab initio_ model based on
density functional theory, analogous to most of the current efforts at
parameterizing ANN potentials. In such a setting, the exact model for the
long-ranged interactions is not known, allowing us to test the sensitivity to
our choice of interaction model. We consider a relatively simple heterogeneous
environment, a liquid in contact with its vapor using a standard slab
geometry, shown in Fig. 1. The width of the slab, equal to the length between
the two Gibbs dividing surfaces, is defined as $2\ell$ and is used throughout
to scale the $z$ axis perpendicular to the interface to compare systems with
slightly different densities. All ANN calculations are done in LAMMPS.Plimpton
(1995)

Figure 1: Characteristic snapshot of the slab geometries of the dipole fluid
(top) and water (bottom), indicating the lengthscale $\ell$, that indicates
the half width of the liquid slab.

### Asymmetric dipolar fluid

The forcefield we fit a ANN to reproduce is a single-component dipolar fluid
constructed from a flexibly bonded dimer, adapted from a model of CO,Daub,
Åstrand, and Bresme (2014) made flexible to avoid difficulties in leaning the
rigid body constraints. The bonding potential is taken to be harmonic, with
the form

$U_{\mathrm{B}}(r_{AB})=\frac{1}{2}k_{\mathrm{bond}}\left(r_{AB}-r_{\mathrm{bond}}\right)^{2}$
(1)

where $r_{AB}$ is the displacement between the $A$ and $B$ sites of the dimer,
held at an equilibrium displacement distance $r_{\mathrm{bond}}=1.2782$Å with
force constant $k_{\mathrm{bond}}=100\mathrm{kcal}/\mathrm{mol}$ Å2. By virtue
of the empirical forcefield, the short- and long-ranged interactions are known
explicitly. These interactions are comprised of a Lennard-Jones potential on
each site and compensating point charges. The Lennard-Jones potential for
pairs $i-j$ takes the form

$U_{\mathrm{LJ}}(r_{ij})=4\epsilon\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]$
(2)

where $r_{ij}$ is the displacement between sites $i$ and $j$, and
$\epsilon=0.08413\,\mathrm{kcal/mol}$ sets the energy scale for the potential.
The dimer is made asymmetric to break charge inversion symmetry so that the
extended interface can support a net polarization. This is done by choosing
the diameters of the $A$ and $B$ sites to be different, with
$\sigma_{AA}=2\sigma_{BB}=2.18113$Å, and $\sigma_{AB}=3\sigma_{BB}/2$. The
Lennard-Jones potential is truncated and shifted at $r_{c}=11.45$Å and thus is
considered a short ranged interaction. Finally, each site is charged and thus
interacts electrostatically with pair potential,

$U_{\mathrm{C}}(r_{ij})=\frac{q_{i}q_{j}}{4\pi\epsilon_{0}r_{ij}}$ (3)

with site charges $q_{A}=-q_{B}=0.21q_{e}$ where $q_{e}$ is the charge of the
electron and $\epsilon_{0}$ is the permittivity of free space. The Coulomb
potential is not truncated and is evaluated using a particle-particle mesh
Ewald summation. The net neutrality and asymmetry of the dimer results in
asymptotic dipole-dipole interactions at long distances, and the flexible bond
endows it with a low frequency polarizability. The two sites are also treated
as point masses with mass 28 amu.

The liquid-vapor interface is studied at 60 K or a reduced temperature of
$T^{*}=k_{\mathrm{B}}T/\epsilon=1.4$, which is below the critical point for
the rigid model, reported to be $T_{c}^{*}=2.0$. The temperature was
controlled using a Langevin thermostat with time constant of 10 ps.
Simulations were initialized from a regular lattice with 1200 molecules or
$N=2400$ atom and equilibrated for 10ns with a constant pressure simulation at
1 atm using the Nose-Hoover algorithm. The simulation box was then extended by
5 nm in the $+z$ and $-z$ directions, creating two liquid-vapor interfaces
perpendicular to the $z$ axis and a 10 nm vacuum gap. The final simulation box
dimensions were $L_{x}=L_{y}=2.157$nm, $L_{z}=15.808$nm. Training data for
fitting the ANN were selected at 50ps intervals from a 100ns simulation at
constant temperature and volume, allowing 1 ns for equilibration.

### revPBE-D3 water

We have also trained a set of ANN models to a many body interatomic potential
energy surface generated for water described by density functional theory. The
primary data sets for the training were generated from _ab initio_ molecular
dynamics simulations using the Gaussian Plane Wave implementation in
CP2K.VandeVondele  _et al._ (2005) All _ab initio_ molecular dynamics
simulations were carried out using the revised version of the PBE
functionalZhang and Yang (1998) along with empirical dispersion correction
(Grimme D3).Grimme  _et al._ (2010) We used a molopt-DZVP basis set and a
plane wave cut-off of 1200 Ry. This large cutoff was needed to converge the
virial component of the pressure. The core electrons were described with GTH
pseudopotential.Goedecker, Teter, and Hutter (1996)

The liquid-vapor interface is studied at 300 K. The temperature was controlled
using a Langevin thermostat with time constant of 1 ps. Simulations were
initialized from a regular lattice with 340 molecules or $N=1020$ atoms and
equilibrated for 20 ps with a constant pressure simulation at 1 atm using the
Nose-Hoover algorithm. The simulation box was then extended by 1.5 nm in the
$+z$ and $-z$ directions, creating two liquid-vapor interfaces perpendicular
to the $z$ axis and a 3nm vacuum gap. The final simulation box dimensions were
$L_{x}=L_{y}=2$nm, $L_{z}=5.0$nm. Training data for fitting the ANN were
selected at 0.005ps intervals from a 25ps simulation at constant temperature
and volume, allowing 20 ps for equilibration. A total of 5000 configurations
were used in the training.

## ANN training methods

We employ the DeePMD scheme to evaluate and train our ANN potentials.Wang  _et
al._ (2018); Zhang  _et al._ (2018b) As with almost all neural network
potentials, DeePMD assigns atomic energies and forces using a local
mathematical description of the environment.Zhang  _et al._ (2018a) The
initial representation of atom $i$’s environment is given by a matrix
$\bm{\mathcal{R}}^{i}\in\mathbbm{R}^{N_{i}\times 4}$ where $N_{i}$ is the
number of atoms within $r_{c}=8$Å of atom $i$. Each row of
$\bm{\mathcal{R}}^{i}$ represents one neighbor and is given by
$\mathcal{R}^{i}=\\{s(r_{ji}),\hat{x}_{ji},\hat{y}_{ji},\hat{z}_{ji}\\}$ where
$\hat{x}_{ji}=s(r_{ji})x_{ji}/r_{ji}$ is a scaled unit vector. Analogous
definitions hold for $\hat{y}_{ji}$ and $\hat{z}_{ji}$, which together convey
relative orientational information. The scaling function $s(r_{ji})$ assigns
greater weight to the nearest neighbors

$s(r_{ji})=\begin{cases}\frac{1}{r_{ji}},&r_{ji}<r_{cs}\\\
\frac{1}{r_{ji}}\\{\frac{1}{2}\cos[\pi\frac{r_{ji}-r_{cs}}{r_{c}-r_{cs}}]+\frac{1}{2}\\},&r_{cs}<r_{ji}<r_{c}\\\
0,&r_{ji}>r_{c}.\end{cases}$ (4)

The parameter $r_{cs}=7.9$Å defines a switching radius beyond which a cosine
function is used to take the scaling factor smoothly to zero at $r_{c}$. This
smoothness guarantees differentiability of the scaling function, which is
required for well-conditioned minimization and energy conserving dynamics.

Another requirement for reliable ANN potentials is that environment
representations should be invariant to global symmetries of the system
including translation, rotation and atomic permutation.Behler and Parrinello
(2007); Glielmo, Sollich, and De Vita (2017); Grisafi  _et al._ (2018) DeePMD
achieves this invariance by converting $\bm{\mathcal{R}}^{i}$ to an invariant
feature matrix $\mathcal{D}^{i}\in\mathbbm{R}^{M_{1}\times M_{2}}$ using an
encoding neural network whose parameters are optimized simultaneously with the
network representing the potential energy.Zhang  _et al._ (2018a) We employed
a 3-layer encoding network containing 25, 50 and 100 neurons in each layer,
with a fourth axis layer containing 16 neurons. For the dipolar model, the
main fitting network is also a three layer neural network containing 300
neurons per layer. For the _ab initio_ water model, the network is a three
layer neural network containing 600 neurons per layer.

### Loss function for optimization

During training, a loss function is minimized with respect to all neural
network parameters using the Adam stochastic gradient descent algorithm.Kingma
and Ba (2015) For the $n$th optimization step, the loss function is given by
an average of a set of $N_{\rm batch}$ configurations, denoted $\mathbf{R}$,

$\mathcal{L}(n)=\frac{1}{N_{\rm batch}}\sum_{\mathbf{R}}^{N_{\rm
batch}}L(\mathbf{R},n)$ (5)

with a step-dependent function

$L(\mathbf{R},n)=p_{e}(n)\Delta
U^{2}+\frac{p_{f}(n)}{3N}\sum_{i}^{N}\Delta\mathbf{F}_{i}^{2}+\frac{p_{\xi}(n)}{9N}\Delta\Xi^{2}\,.$
(6)

Here $\Delta U(\mathbf{R})=U_{\rm ANN}(\mathbf{R})-U_{\rm ref}(\mathbf{R})$ is
the prediction error of the ANN for the total configuration potential energy,
$\Delta\mathbf{F}_{i}=|\mathbf{F}_{i,{\rm ANN}}(\mathbf{R})-\mathbf{F}_{i,{\rm
ref}}(\mathbf{R})|$ is the error for the force on atom $i$, and
$\Delta\Xi=||\Xi_{{\rm ANN}}(\mathbf{R})-\Xi_{{\rm ref}}(\mathbf{R})||$ is the
elementwise error in the virial matrix,
$\mathbf{\Xi}$=$-\sum_{i}\mathbf{R}_{i}\otimes\mathbf{F}_{i}/2$. The three
prefactors $p_{e},p_{f}$ and $p_{\xi}$ are weights that vary with optimization
step number $n$ according to
$p(n)=p(n_{f})[1-r_{l}(n)/r_{l}(0)]+p(0)[r_{l}(n)/r_{l}(0)]$, where $r_{l}(n)$
is the learning step size and $p(n_{f})$ is a specified final value. The
learning rate decreases during optimization following
$r_{l}(n)=r_{l}(0)d_{r}^{(n/d_{s})}$, where $d_{r}$ and $d_{s}$ are termed the
decay rate and decay steps, respectively. In the limit $N_{\rm batch}$ equals
the whole training set, $\mathcal{L}(n)$ becomes deterministic, but in general
a subset of the total configurations are sampled randomly.

In addition to possible representability problems with employing ANNs with
truncated interactions to model interfaces, there is a practical training
problem. The proportion of atoms representing interfacial environments is
small, which can lead to under-learning of forces and higher prediction errors
for interfacial atoms relative to those in bulk-like environments. We found
that the stability of molecular dynamics and the accuracy of interfacial
properties are particularly sensitive to the maximum error of the neural
network, as well as to the absolute average error. To mitigate this
environment bias we have found it useful to introduce a weighting function to
the loss function to homogenize high-error atomic environments. Specifically,
we dress the squared force error

$\Delta\mathbf{F}_{i}^{2}\rightarrow\Delta\mathbf{F}_{i}^{2}\left[w-\delta
w\,{\rm tanh}\left(g-\Delta\mathbf{F}_{i}^{2}/\bar{f}^{2}\right)\right]$ (7)

that saturates at large forces. We have empirically found $w=502,~{}\delta
w=2,\bar{f}=0.022e\mathrm{V/}$Å and $g=14$ works well in training these
systems.

There are many tunable hyperparameters associated with the training algorithm.
We have found that parameter choice rarely makes a large difference to the
minimal loss function obtained during a training run, but it often does affect
the number of steps required to reach this value. Moreover, the variance in
training quality and model properties between ANNs trained using different
hyperparameters is no greater than the variance between training processes
with the same hyperparameters but different random seeds. Example
hyperparameters for training the dipole forcefield model are shown in Table 1.

Hyperparameter |  step 1 |  step 2 |  step 3
---|---|---|---
$N_{\rm batch}$ | 30 | 1 | 1
$l_{r}(0)$ | $10^{-3}$ | $10^{-3}$ | $10^{-4}$
$d_{r}$ | 0.999 | 0.999 | 0.999
$d_{s}$ | 10 | 10 | 1000
$p_{e}(0)$ | 0.01 | 0.01 | 0.01
$p_{e}(n_{f})$ | 1.0 | 0.01 | 0.01
$p_{f}(0)$ | 100.0 | 1.0 | 1.0
$p_{f}(n_{f})$ | 1.0 | 1.0 | 1.0
$p_{\xi}(0)$ | 1.0 | 0.0 | 0.0
$p_{\xi}(n_{f})$ | 1.0 | 0.0 | 0.0
Table 1: Example hyperparameters for training the dipole model.

### Explicit Long-ranged ANN models

In order to understand the impact of long-range forces on interfacial
structure and thermodynamics, we develop models represented by standard short-
ranged ANNs and those supplemented with an additional model of the long-ranged
interaction. The former, which we refer to as ANN-SR for short-ranged
interactions, is trained on the full potential energy function $U(\mathbf{R})$
of the reference model and its derivatives which provide information on the
forces and virial. For the dipolar fluid model, the total potential consists
of a sum over the explicit pairwise bonding, Lennard-Jones and Coulomb
potentials in Eqs. 1-3. For the revPBE-D3 water model, the total potential
energy derives from an expectation over the minimized electron density, and
the forces are derived from the Hellman-Feynman theorem.VandeVondele  _et al._
(2005)

The model with explicit long-ranged interactions, which we refer to as ANN-LR
for long-ranged interactions, is trained on only a piece of the reference
potential energy, while an additional piece is treated explicitly with an
empirical model. In contrast to previous work,Yue  _et al._ (2021) we do not
train the ANN-LR models on the difference between the total potential energy
and electrostatic potential. Rather we subtract only the long-ranged, slowly
varying interactions. This is accomplished using a construction supplied by
LMF theory by defining the short ranged potential as,Rodgers  _et al._ (2006)

$U_{\rm
SR}(\mathbf{R})=U(\mathbf{R})-\sum_{i<j}\frac{q_{i}q_{j}}{4\pi\epsilon_{0}r_{ij}}{\rm
erf}(\kappa r_{ij})$ (8)

where the error function, ${\rm erf}(\kappa r_{ij})$, has a natural
lengthscale $\kappa^{-1}$ and multiplies a pairwise Coulomb model for the
long-range interactions. The ANN-LR models are trained to reproduce only
$U_{\rm SR}(\mathbf{R})$ and its derivatives, while simulations are run using
that forcefield combined with an exact treatment of the remainder. Several
efficient methods exist for evaluating the remainder, and we employ the linear
scaling particle particle mesh Ewald sum.Hockney and Eastwood (2021) For the
dipolar fluid model, this procedure is equivalent to training an ANN on that
part of the potential energy typically evaluated in real-space, including the
bonding, Lennard-Jones, and part of the Coulomb potential, while still
explicitly evaluating the reciprocal space part of the Coulomb potential.
While in the dipolar fluid model, $U_{\mathrm{C}}(r_{ij})$ is known
explicitly, for the revPBE-D3 water model, the long-range-part of the
interactions is not known. We model it as a sum of Gaussian charge densities
centered on each atom with fixed amplitude. Below, we explore the sensitivity
of this approach to the choice of empirical amplitude, but most of our results
employ charges taken from the empirical SPC/E forcefield with hydrogen charge
$q_{H}=0.4238q_{e}$.Berendsen, Grigera, and Straatsma (1987)

Within LMF theory, it is known that $\kappa$ needs to be chosen small enough
such that in the absence of the long-ranged potential pair correlations are
still well-described.Rodgers and Weeks (2008) We find that $\kappa^{-1}=8.0$Å,
equivalent to the cutoff $r_{c}$ used in the neural network representations,
works well for both the dipole model as well as the _ab initio_ water model.
In principle, the charges and $\kappa$ could be optimized alongside the ANN,
but for simplicity we do not consider doing so here and employ these fixed
values.

Figure 2: A) Ensemble-averaged learning curves for ANN models of the dipole
fluid model. Dark lines correspond to the loss function including only the
contribution from the force error. Shaded areas indicate the standard
deviation of the set of models. B) Average errors in $z$-force prediction as a
function of atomic $z$ coordinate for ANN-SR and ANN-LR models normalized by
the root-mean squared error for each model. Configurations used in the
ensemble average were drawn from the ANN training set used for the models.
Both panels A and B are averaged over 20 training runs, or their resultant
models, initialized with different random seeds.

Our approach provides a clear separation between the short-range component of
the energy that is representable by a ANN and the long-range component that is
not. We find that this separation provides an accurate and efficient means of
modeling extended interfaces with ANNs. Figure 2A) shows typical initial
learning curves for ANN-SR and ANN-LR models of the dipole fluid,
demonstrating that the latter attains a significantly lower loss function,
leading to more accurate predictions of properties and greater consistency
between models. The converged loss function of the ANN-LR model is
approximately an order of magnitude smaller than that of the ANN-SR model.
This difference is expected, since the ANN-SR attempts to fit correlations
between short-range structure and long-range energies that are not easily
represented by the truncated symmetry functions.

Figure 2B) illustrates that this residual error in the ANN-SR models is due to
the presence of the interface by plotting the force error averaged over an
equilibrium trajectory conditioned by the $z$ value of the atom and normalized
by the mean squared error of each model. The errors in the bulk region of the
slab are essentially the same in both models and very small, which indicates
that both models are able to learn the short-ranged part of the potential
well. Both models also exhibit larger errors at the interface compared to the
bulk, which is a consequence of interfacial environments being
underrepresented in the training set for the slab geometry. However, the
errors in the ANN-SR model are much larger than the ANN-LR model at the
interface, indicating that the multiplicity of long-ranged environments
hinders the training process. Through subsequent active learning steps
described below, we are able to optimize both models to have comparable
residual force errors below 60 meV/Å. However, this procedure does not
eliminate the difficulty in using the ANN-SR to describe extended interfaces.

### Active learning steps

For both ANNs with and without long-ranged interactions and for both reference
potentials, we trained a set of 3 models using a sequence of active learning
steps. Data for each of the steps was taken from canonical ensemble
simulations of a liquid slab in contact with its vapor. The details of the
training of the ANN models depended slightly on the reference potential
employed.

For the dipole model, an initial ANN was constructed using approximately 10000
optimization steps with $N_{\rm batch}=30$. In all cases, the liquid densities
and surface tensions predicted by the resulting models were much too low, and
several models did not give a stable interface. To remedy these issues, we
retrained using an additional dataset containing a set of perturbed
configurations. To generate these configurations, we varied $\sigma_{AB}\in$
[3.0Å,3.2717Å] and $k_{\rm bond}\in$[20.0kcal/mol Å,100.0 kcal/mol Å] and ran
short molecular dynamics simulations. A total of 724 perturbed configurations
were extracted from a combined 400ns of simulation. In this step the batch
size was reduced to 1 for greater training efficiency. This training step was
followed by an explicit active learning step, where 3317 configurations with a
particularly high average or maximum force error (7eV/Å or 14eV/Å) were taken
from a simulation employing the existing ANN model. These configurations were
added to the original equilibrium dataset, and the loss function was
reoptimized.

Figure 3: Bulk properties of the dipole fluid. (Top) ANN-SR model predictions
of the pair distribution functions (left) and center-of-mass velocity
autocorrelation functions. (Bottom) ANN-LR model predictions of the pair
distribution functions (left) and velocity autocorrelation functions (right).
Solid lines denote the reference force-field calculation and shaded regions
denote the range of predictions from 3 ANN models. Figure 4: Bulk properties
of revPBE-D3 water. (Top) ANN-SR model predictions of the pair distribution
functions (left) and center-of-mass velocity autocorrelation functions.
(Bottom) ANN-LR model predictions of the pair distribution functions (left)
and velocity autocorrelation functions (right). Solid lines denote the
reference _ab initio_ calculation and shaded regions denote the range of
predictions from 3 ANN models.

For the water model, an initial ANN was trained using 5000 configurations
extracted from a 25 ps AIMD slab simulations and 5000 configurations of bulk
water. Initially the loss function was optimized to reduce the energy and
force error. This initial ANN was then retrained by optimizing the virial
error in addition to energy and force errors using 500 configurations,
followed by an active learning step where approximately 10000 configurations
were added to the previous data set and only energy and force error was
optimized. Those new 10000 configurations were extracted from five 1 ns
trajectories generated by an ensemble of 3 ANN models using query by
committee.Krogh, Vedelsby  _et al._ (1995) In addition to that, another 2000
configurations from that 5 ns trajectory were used to optimize the virial
error. These final models gave well correlated force and virial predictions.

## Validation of bulk liquid properties

The similar loss functions obtained for ANN-SR and ANN-LR at the end of the
active learning cycles ensures that models obtain similar average accuracy in
their prediction of the forces relative to the reference data. Such agreement
is a necessary but not sufficient condition for reliable predictions of the
thermodynamics or dynamics. We first focus on comparing bulk ensemble
properties of the ANN models with the reference potential. In a homogeneous
liquid, such as the interior of our slab geometry, LMF theory suggests that
slowly-varying long-ranged force contributions largely cancel and short-range
forces dominate the structure of the fluid. We thus expect both ANNs to
describe these properties well. Indeed, substantial previous work on ANN
models of bulk fluids indicate that such potentials can robustly recover both
the thermodynamics and dynamics of their reference systems.Cisneros  _et al._
(2016); Li  _et al._ (2017); Unke and Meuwly (2019); Cubuk  _et al._ (2017);
Zhang  _et al._ (2018a); Behler and Parrinello (2007); Behler (2011); Wang
_et al._ (2018); Behler  _et al._ (2008); Cheng  _et al._ (2018); Morawietz
_et al._ (2016); Bonati and Parrinello (2018); Yang  _et al._ (2021); Galib
and Limmer (2021); Käser, Unke, and Meuwly (2020); Liu  _et al._ (2018)

Figures 3 and 4 confirm that both ANN-SR and ANN-LR models are well able to
describe the equilibrium structure, as encoded by the pair distribution
functions $g(r)$ of the corresponding fluids. Further, both ANN-SR and ANN-LR
models are able to describe the equilibrium diffusive dynamics as encoded by
the molecular center of mass velocity autocorrelation functions,
$\langle\mathbf{v}(0)\cdot\mathbf{v}(t)\rangle$. Both sets of properties are
evaluated excluding atoms within 10 Å (dipole) or 5.0 Å (water) of the Gibbs
dividing interface, mitigating the influence from the regions of broken
translational symmetry where unbalanced long-ranged interactions may be
important. For models of the dipolar fluid and _ab initio_ water, we find
similar quantitative agreement between the ANN-SR, ANN-LR models and benchmark
calculations. Further all three independently-optimized ANN models for each
type, ANN-SR or ANN-LR, are consistent with little variability amongst them.

## Importance of long-ranged interactions for interfacial properties

Within the context of LMF theory, the inability of the ANN-SR models to
describe long-range interactions suggests that they should struggle to
accurately predict interfacial properties where such interactions become
unbalanced. We can anticipate for which properties in particular the ANN-SR
model is likely to struggle by considering a result of LMF theory. The
potential decomposition in Eq. 8 presupposes the important long-range, slowly
varying part of the potential is Coulombic in origin. For the simple dipole
fluid and water, this is a reasonable assumption as for both their leading
order asymptotic interactions are dipole-dipole forces. In such a case, a
truncated model can be constructed that accurately recovers the interfacial
properties of the reference system provided an additional external potential.
The form of the external potential, denoted $V_{R}(\mathbf{r})$, is
analytically known and given by an electrostatic potentialRodgers and Weeks
(2008)

$V_{R}(\mathbf{r})=\int d\mathbf{r}^{\prime}\frac{{\rm
erf}(|\mathbf{r}-\mathbf{r}^{\prime}|/\sigma)}{4\pi\epsilon_{0}|\mathbf{r}-\mathbf{r}^{\prime}|}\sum_{i}q_{i}\rho_{i}(\mathbf{r}^{\prime})$
(9)

where the slowly varying potential ${\rm
erf}(|\mathbf{r}-\mathbf{r}^{\prime}|/\sigma)/4\pi\epsilon_{0}|\mathbf{r}-\mathbf{r}^{\prime}|$
is convoluted with the charge density
$\sum_{i}q_{i}\rho_{i}(\mathbf{r}^{\prime})$. For the slab geometry,
translational invariance in the plane of the interface renders
$V_{R}(\mathbf{r})$ a function of $z$ only.Cox (2020)

Since $V_{R}(\mathbf{r})$ is an external electrostatic potential, properties
correlated with polarization fluctuations will be most affected by the absence
of long-range interactions. Properties weakly correlated with an electrostatic
potential are likely to be well described by a model with purely short-ranged
interactions, like the ANN-SR models. Figures 5 and 6 illustrate the number
density and orientational order profiles relative to the liquid-vapor
interfaces of our slab geometries. Consistent with the expectation that
density fluctuations are weakly correlated with polarization
fluctuations,Remsing and Weeks (2013) we find that both the ANN-SR and ANN-LR
are able to reasonably reproduce the sigmoidal density profiles,
$\rho(z)/\rho$, of the reference dipolar fluid and _ab initio_ water models
they have been trained on. This agreement can be understood as the weak
response of the mass density of neural molecules to an external electric
field. The bulk density $\rho$ for the dipolar fluid model is $\rho=0.036$Å-3
while that for the revPBE-D3 water model $\rho=0.96\mathrm{g/cm}^{3}$. The
value of the water density is in reasonable agreement with recent
estimates,Galib  _et al._ (2017) which found a strong dependence on the basis
set. Further, as we observed for the bulk liquid properties, there is a small
variance between the ensemble of the 3 independently trained potentials.

Figure 5: Interface properties of the dipolar fluid. (Top) ANN-SR model
predictions of the density (left) and orientational order (right). (Bottom)
ANN-LR model predictions of the A-site density (left) and orientational order
(right). Solid black lines denote the reference _ab initio_ calculation and
shaded regions denote the range of predictions from 3 ANNs. The different
colored lines indicate different ANN-SR models. Figure 6:  Interface
properties of revPBE-D3 water. (Top) ANN-SR model predictions of the density
(left) and orientational order (right). (Bottom) ANN-LR model predictions of
the density (left) and orientational order (right). Solid black lines denote
the reference _ab initio_ calculation and shaded regions denote the range of
predictions from 3 ANNs. The different colored lines indicate different ANN-SR
models.

The orientational order profiles shown in Figures 5 and 6, are computed by
evaluating the angle, $\theta$, between the molecule’s dipole and the $z$
axis. For the dipolar fluid and the _ab initio_ water model, the orientational
order profiles are generally poorly predicted by the ANN-SR models, but well
recovered by the ANN-LR models. This is rationalized within the framework of
LMF theory as the molecular orientation is responsible for polarization
fluctuations in the fluid, which are therefore sensitive to the unbalanced
electrostatic potential generated by truncating the dipolar interactions.Cox
(2020) For both the dipolar fluid and water, we find a wide variety of
behaviors in the ANN-SR models. These include models with relatively little
orientational alignment throughout the slab, and models with persistent
alignment permeating throughout even the bulk of the liquid. This persistent
alignment, since in both the ANN-SR trained on the dipolar fluid and the water
but not for any ANN-LR model, has also been observed recently in truncated
forcefield models of water.Cox (2020) However, as with the bulk properties,
the ANN-LR models exhibit little variability between the ensemble of the 3
independently trained potentials.

Aside from structural measures of the interface, we have also considered the
impact of long-ranged interactions on the liquid-vapor surface tension,
$\gamma$. We have computed the surface tension using the approach by Kirkwood
employing a difference in the diagonal components of the virial,Frenkel and
Smit (2001)

$\gamma=\left\langle\Xi_{zz}-\frac{1}{2}\left(\Xi_{xx}+\Xi_{yy}\right)\right\rangle\frac{L_{z}}{2}$
(10)

where $\Xi_{ii}$ is the $i$th component of the virial. Previous reports have
indicated that short-ranged ANN models can accurately predict the surface
tension.Wohlfahrt, Dellago, and Sega (2020) We find that this is indeed the
case, with the reference surface tension value computed by the empirical
forcefield dipolar fluid model following within the range of both ANN-SR and
ANN-LR models. The value of $\gamma$ for revPBE-D3 water is not known, since
it is too expensive to evaluate directly, but previous ANN results suggest
that it is 68 mN/m at 300 K.Wohlfahrt, Dellago, and Sega (2020) Our ANN-SR
estimates are slightly lower than this value. Unlike the orientation, which
will depend linearly on an external field generated by the truncated
interactions, the leading order dependence of the surface tension is
quadratic. This is because the slab has inversion symmetry and cannot depend
on the sign of the applied field. Thus the agreement between the reference
calculations and the ANN-SR models suggests that the majority of $\gamma$ is
determined locally, while the remainder is small enough that its quadratic
correction is negligible.

Model |  $\gamma_{\rm Dipole}$ (mN/m) |  $\gamma_{\rm H_{2}O}$ (mN/m)
---|---|---
Reference |  4.5 $\pm$ 0.3 | 68Wohlfahrt, Dellago, and Sega (2020)
ANN-SR |  5.2 $\pm$ 0.8 |  54 $\pm$ 6
ANN-LR |  5.0 $\pm$ 0.4 |  66 $\pm$ 14
Table 2: Surface tension values for different neural network models
parameterized for the dipole forcefield, $\gamma_{\rm Dipole}$ and the _ab
initio_ water, $\gamma_{\rm H_{2}O}$. The error in the reference is
statistical, and those of the ANN models are calculated from the range of 3
models, which is larger than the statistical uncertainty.

## ANN-SR approximation of the LMF potential

The reference training data includes force contributions from atoms outside
the symmetry function cutoff for the ANN-SR models. The ANN still attempts to
fit these contributions to its symmetry functions, so natural question is
whether the ANN attempts to learn the LMF potential. The orientational
profiles in Figs. 5 and 6 suggest that some of the models manage to learn some
degree of the long-range forces since not all of them generate a net
polarization in the bulk-like region, and they produce a smaller polarization
at the interface than do the truncated models without the added explicit long-
range potential.Cox (2020)

Figure 7 shows the learned force separation, which compares the implicit long-
range component of the ANN-SR potential with an analytical approximation to
the force for the exact dipole model. Specifically, we evaluate the LMF
potential by solving Eq. 9 using the charge density computed from the
reference forcefield calculations and take its numerical gradient in the $z$
direction. This is the one-shot approach to solving the self-consistent LMF
equations.Rodgers and Weeks (2008); Cox (2020) The ANN-SR residual forces,
computed by subtracting the short-range part of the reference forcefield from
the ANN-SR force, resemble the LMF potential but do not exactly recover it.
This indicates that the ANN-SR models attempt to approximate a molecular field
potential, which is surprising given the scale separation in the reference
data is implicit. However, they are not particularly successful and there is
significant variation in the scale and position of the implicit long-ranged
forces as a result. While some of the models reproduce the results of LMF
theory, most over- or under-estimate the restructuring potential
significantly. It is possible that tailoring loss functions to recover this
LMF potential may be successful for specific systems and geometries, however
we expect that the implicitly-learned long-ranged force will be generally
poorly transferrable given the LMF potential depends on the system geometry.

Figure 7: Estimates of the long-ranged force learned implicitly by ANN-SR
potentials, compared to the LMF force. For the three ANN models, the ensemble
average is performed over a long trajectory generated by the respective model,
and averaged over all A-sites within each $z$ interval.

## Sensitivity of ANN-LR to the charge model

Finally, in order to incorporate explicit long-ranged interfaces for the ANN-
LR model of _ab initio_ water, we used a simple empirical form. Specifically,
we associate the atom centers with fixed charges. Throughout we have used
charges taken from the SPC/E water model and found reasonable agreement with
the reference calculations. In the limit of perfect prediction by the ANN,
provided the ANN could approximate the LMF potential accurately, the use of
the long-range potential model is inconsequential. This is because we subtract
it in the fitting procedure and then add it back into the simulations.
However, no ANN is perfectly predictive, and as we have we observed, the ANNs
are not able to robustly learn the LMF potential. Consequently it is important
to explore the sensitivity, or insensitivity, of our results to our specific
choice of long-range model.

From our results validating the bulk liquid properties with both the ANN-SR
and ANN-LR models, we expect that these properties are insensitive to the
choice of long-range potential by virtue of those contributions canceling in
the part of the slab that obtains approximate transitional invariance. We test
this assertion by training ANN-LR models with range of different charge
models, up to $\pm 10\%$ of SPC/E values employed in the majority of the
calculations presented here. Indeed we found that he pair distribution
functions and velocity auto-correlation functions are within the spread of the
original ANN-LR models, indicating little sensitivity to the model long-ranged
forces.

Its less clear _a priori_ whether the interfacial properties are insensitive
to the long range model employed. In order to quantify the dependence to the
fixed charges, we again consider varying the magnitude of the charges by $\pm
10\%$ and study the effect on the density and orientational profiles.
Specifically, we fit the density profile to a functional form

$\rho(z)=\rho\left(\frac{1}{1+e^{-(z+\ell)/w_{\rho}}}-\frac{1}{1+e^{-(z-\ell)/w_{\rho}}}\right)$
(11)

with shape descriptors $\rho$ for the bulk density and $w_{\rho}$ for the
width of the liquid-vapor interface. Additionally we fit the orientational
profile to the functional form

$\langle\cos(\theta)\rangle=\alpha\left(e^{-(z+\ell)^{2}/w_{\theta}^{2}}-e^{-(z-\ell)^{2}/w_{\theta}^{2}}\right)$
(12)

with shape descriptors $\alpha$ for the amplitude of the polarization and
$w_{\theta}$ for the width of the polarized region at the interface. Both
functional forms provide good fits across the range of ANN-LR models
constructed. We train, simulate and fit 5 ANN-LR models using the previously
describe active leaning procedure with charges on the hydrogen
$q^{*}_{H}/q_{e}=0.3814,0.4238,0.4662$, and fit the density and orientational
profiles.

Summary statistics for $\rho,\alpha,w_{\rho}$ and $w_{\theta}$ for these
models are shown in Fig. 8. For all of the shape descriptors, the variability
within a set of models with the same charge is larger then the variation
between models with different charges. The two width parameters vary on sub-
angstrom scales, with the typical values of $w_{\rho}\approx 0.9$Å and
$w_{\theta}\approx 2.2$Å. There is no systematic trend in the amplitude of the
orientational order $\alpha$, though we do find that the bulk density is more
variable between models for the lower and high charge cases than for the SPC/E
charges employed. We are left with the conclusion that the internal
variability of equivalent ANN models with different stochastic training
histories is a much larger effect than the variation between different long-
ranged force models, and thus that the ANN-LR models are relatively
insensitive to the specific charges employed.

Figure 8: Box plots of the shape parameters for the density (top) and
orientation (bottom) profiles denoting the median, first quartile and range
over 5 independently trained ANN-LR with varying fixed charge models for a
slab of revPBE-D3 water.

## Conclusions

ANN models are capable of describing complex, many-bodied interactions,
enabling their deployment in a wide range of chemical phenomena. Currently,
ANN potentials are overwhelmingly formulated with atomic descriptors that are
local and short-ranged. Here, we have used LMF theory as a framework to
understand in what cases such locality can offer a good approximation to
equilibrium structure and dynamics, and when we should expect the neglect of
long-ranged interactions to matter. In systems with translational invariance,
long-range slowly varying forces like those between two molecules with
permanent dipoles in a bulk fluid largely cancel. As such the neglect of long-
range forces is inconsequential. However, at extended interfaces like that
between a liquid and its vapor, truncated interactions result in unbalanced
forces that can polarize and distort local molecular structure. While in
principle, ANNs may be able to uncover an approximation to an effective
external field that can correct for these unbalanced interactions, we have
found in practice this is difficult to achieve. The resultant models are
typically inaccurate, difficult to train, and predict interfaces that are
overly polarized.

To mediate the failings of truncated ANNs, we have proposed training models
only on the short-range interactions and employing a model of the long-ranged
forces. Rather than training models with a full electrostatic force model
subtracted from them, we use insights from LMF theory to subtract only the
part traditionally evaluated in reciprocal space. This appraoch leads to
models with a natural separation of length scales, which are able to reproduce
reference models based on an empirical fixed charge potential as well as one
derived from density functional theory. This separation has the added benefit
that local interactions are done completely by the neural network, while long-
ranged interactions are done efficiently through traditional Ewald summations,
avoiding the need to sum independent local interactions. We have employed
simple models for the long ranged interactions based on fixed point charges,
with the expectation that higher multipoles decay faster and are accommodated
by the short-ranged ANN. Further, the local description of the ANN employed
here is flexible enough to generate models that are largely insensitive to the
specific charge model employed. With an ability to study extended interfaces
with ANNs robustly, a wide range of phenomena is tractable including
understanding the effects of extended interfaces on chemical
reactivity.Niblett and Limmer (2021); Schile and Limmer (2019); Kattirtzi,
Limmer, and Willard (2017); Benjamin (2015); Mundy and Kuo (2006); Lewis  _et
al._ (2011) Addressing such questions will undoubtedly advance fields such as
atmospheric chemistry and catalysis where strongly inhomogeneous systems are
ubiquitous.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Acknowledgements

This material is based on work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research,
Scientific Discovery through Advanced Computing (SciDAC) program, under Award
No. DE-AC02-05CH11231. This research used resources of the National Energy
Research Scientific Computing Center (NERSC), a U.S. Department of Energy
Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

## References

## References

  * Behler (2016) J. Behler, “Perspective: Machine learning potentials for atomistic simulations,” The Journal of chemical physics 145, 170901 (2016).
  * Noé  _et al._ (2020) F. Noé, A. Tkatchenko, K.-R. Müller, and C. Clementi, “Machine learning for molecular simulation,” Annual review of physical chemistry 71, 361–390 (2020).
  * Bartók  _et al._ (2017) A. P. Bartók, S. De, C. Poelking, N. Bernstein, J. R. Kermode, G. Csányi, and M. Ceriotti, “Machine learning unifies the modeling of materials and molecules,” Science advances 3, e1701816 (2017).
  * Cisneros  _et al._ (2016) G. A. Cisneros, K. T. Wikfeldt, L. Ojamae, J. Lu, Y. Xu, H. Torabifard, A. P. Bartók, G. Csányi, V. Molinero, and F. Paesani, “Modeling molecular interactions in water: From pairwise to many-body potential energy functions,” Chemical reviews 116, 7501–7528 (2016).
  * Li  _et al._ (2017) Y. Li, H. Li, F. C. Pickard IV, B. Narayanan, F. G. Sen, M. K. Chan, S. K. Sankaranarayanan, B. R. Brooks, and B. Roux, “Machine learning force field parameters from ab initio data,” Journal of chemical theory and computation 13, 4492–4503 (2017).
  * Unke and Meuwly (2019) O. T. Unke and M. Meuwly, “Physnet: a neural network for predicting energies, forces, dipole moments, and partial charges,” Journal of chemical theory and computation 15, 3678–3693 (2019).
  * Cubuk  _et al._ (2017) E. D. Cubuk, B. D. Malone, B. Onat, A. Waterland, and E. Kaxiras, “Representations in neural network based empirical potentials,” The Journal of chemical physics 147, 024104 (2017).
  * Zhang  _et al._ (2018a) L. Zhang, J. Han, H. Wang, W. Saidi, R. Car, and W. E, “End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems,” in  _Advances in Neural Information Processing Systems_, Vol. 31, edited by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Curran Associates, Inc., 2018).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, “Generalized neural-network representation of high-dimensional potential-energy surfaces,” Physical review letters 98, 146401 (2007).
  * Behler (2011) J. Behler, “Atom-centered symmetry functions for constructing high-dimensional neural network potentials,” The Journal of chemical physics 134, 074106 (2011).
  * Wang  _et al._ (2018) H. Wang, L. Zhang, J. Han, and W. E, “Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics,” Computer Physics Communications 228, 178–184 (2018).
  * Csányi  _et al._ (2017) G. Csányi, J. R. Kermode, S. De, N. Bernstein, M. Ceriotti, A. P. Bartók, and C. Poelking, “Machine learning unifies the modeling of materials and molecules,” Science Advances 3, e1701816 (2017).
  * Marchand  _et al._ (2020) D. Marchand, A. Jain, A. Glensk, and W. Curtin, “Machine learning for metallurgy i. a neural-network potential for al-cu,” Physical Review Materials 4, 103601 (2020).
  * Bischak  _et al._ (2020) C. G. Bischak, M. Lai, Z. Fan, D. Lu, P. David, D. Dong, H. Chen, A. S. Etman, T. Lei, J. Sun, _et al._ , “Liquid-like interfaces mediate structural phase transitions in lead halide perovskites,” Matter 3, 534–545 (2020).
  * Behler  _et al._ (2008) J. Behler, R. Martoňák, D. Donadio, and M. Parrinello, “Pressure-induced phase transitions in silicon studied by neural network-based metadynamics simulations,” physica status solidi (b) 245, 2618–2629 (2008).
  * Cheng  _et al._ (2018) B. Cheng, E. A. Engel, J. Behler, C. Dellago, and M. Ceriotti, “Ab initio thermodynamics of liquid and solid water,” Proc. Nat. Acad. Sci. USA 116, 1110–1115 (2018), 1811.08630 .
  * Morawietz  _et al._ (2016) T. Morawietz, A. Singraber, C. Dellago, and J. Behler, “How van der waals interactions determine the unique properties of water,” Proceedings of the National Academy of Sciences 113, 8368–8373 (2016).
  * Bonati and Parrinello (2018) L. Bonati and M. Parrinello, “Silicon liquid structure and crystal nucleation from ab initio deep metadynamics,” Physical review letters 121, 265701 (2018).
  * Yang  _et al._ (2021) M. Yang, L. Bonati, D. Polino, and M. Parrinello, “Using metadynamics to build neural network potentials for reactive events: the case of urea decomposition in water,” Catalysis Today (2021).
  * Galib and Limmer (2021) M. Galib and D. T. Limmer, “Reactive uptake of n2o5 by atmospheric aerosol is dominated by interfacial processes,” Science 371, 921–925 (2021).
  * Käser, Unke, and Meuwly (2020) S. Käser, O. T. Unke, and M. Meuwly, “Isomerization and decomposition reactions of acetaldehyde relevant to atmospheric processes from dynamics simulations on neural network-based potential energy surfaces,” The Journal of Chemical Physics 152, 214304 (2020).
  * Liu  _et al._ (2018) Q. Liu, X. Zhou, L. Zhou, Y. Zhang, X. Luo, H. Guo, and B. Jiang, “Constructing high-dimensional neural network potential energy surfaces for gas–surface scattering and reactions,” The Journal of Physical Chemistry C 122, 1761–1769 (2018).
  * Behler (2017) J. Behler, “First principles neural network potentials for reactive simulations of large molecular and condensed systems,” Angewandte Chemie International Edition 56, 12828–12840 (2017).
  * Natarajan and Behler (2016) S. K. Natarajan and J. Behler, “Neural network molecular dynamics simulations of solid–liquid interfaces: water at low-index copper surfaces,” Physical Chemistry Chemical Physics 18, 28704–28725 (2016).
  * Andrade  _et al._ (2020) M. F. C. Andrade, H.-Y. Ko, L. Zhang, R. Car, and A. Selloni, “Free energy of proton transfer at the water–tio 2 interface from ab initio deep potential molecular dynamics,” Chemical Science 11, 2335–2341 (2020).
  * Wohlfahrt, Dellago, and Sega (2020) O. Wohlfahrt, C. Dellago, and M. Sega, “Ab initio structure and thermodynamics of the RPBE-D3 water / vapor interface by neural-network molecular dynamics,” J. Chem. Phys. 153, 144710 (2020).
  * Yue  _et al._ (2021) S. Yue, M. C. Muniz, M. F. Calegari Andrade, L. Zhang, R. Car, and A. Z. Panagiotopoulos, “When do short-range atomistic machine-learning models fall short?” J. Chem. Phys. 154, 034111 (2021).
  * Weeks, Chandler, and Andersen (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, “Role of repulsive forces in determining the equilibrium structure of simple liquids,” The Journal of chemical physics 54, 5237–5247 (1971).
  * Lum, Chandler, and Weeks (1999) K. Lum, D. Chandler, and J. D. Weeks, “Hydrophobicity at small and large length scales,” J. Phys. Chem. B 103, 4570–4577 (1999).
  * Weeks, Selinger, and Broughton (1995) J. D. Weeks, R. L. Selinger, and J. Q. Broughton, “Self-consistent treatment of repulsive and attractive forces in nonuniform liquids,” Physical review letters 75, 2694 (1995).
  * Weeks (2002) J. D. Weeks, “Connecting local structure to interface formation: A molecular scale van der waals theory of nonuniform liquids,” Annual review of physical chemistry 53, 533–562 (2002).
  * Hansen and McDonald (1990) J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids_ (Elsevier, 1990).
  * Rodgers and Weeks (2008) J. M. Rodgers and J. D. Weeks, “Local molecular field theory for the treatment of electrostatics,” Journal of Physics: Condensed Matter 20, 494206 (2008).
  * Remsing, Rodgers, and Weeks (2011) R. C. Remsing, J. M. Rodgers, and J. D. Weeks, “Deconstructing classical water models at interfaces and in bulk,” Journal of Statistical Physics 145, 313–334 (2011).
  * Rodgers  _et al._ (2006) J. M. Rodgers, C. Kaur, Y.-G. Chen, and J. D. Weeks, “Attraction between like-charged walls: Short-ranged simulations using local molecular field theory,” Physical review letters 97, 097801 (2006).
  * Remsing  _et al._ (2014) R. C. Remsing, M. D. Baer, G. K. Schenter, C. J. Mundy, and J. D. Weeks, “The role of broken symmetry in solvation of a spherical cavity in classical and quantum water models,” The journal of physical chemistry letters 5, 2767–2774 (2014).
  * Remsing and Weeks (2013) R. C. Remsing and J. D. Weeks, “Dissecting hydrophobic hydration and association,” The Journal of Physical Chemistry B 117, 15479–15491 (2013).
  * Remsing and Weeks (2016) R. C. Remsing and J. D. Weeks, “Role of local response in ion solvation: Born theory and beyond,” The Journal of Physical Chemistry B 120, 6238–6249 (2016).
  * Cox (2020) S. J. Cox, “Dielectric response with short-ranged electrostatics,” Proceedings of the National Academy of Sciences 117, 19746–19752 (2020).
  * Grisafi and Ceriotti (2019) A. Grisafi and M. Ceriotti, “Incorporating long-range physics in atomic-scale machine learning,” J. Chem. Phys. 151, 204105 (2019).
  * Deng  _et al._ (2019) Z. Deng, C. Chen, X.-G. Li, and S. P. Ong, “An electrostatic spectral neighbor analysis potential for lithium nitride,” npj Computational Materials 5, 75 (2019).
  * Bartók  _et al._ (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, “Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons,” Physical review letters 104, 136403 (2010).
  * Artrith, Morawietz, and Behler (2011) N. Artrith, T. Morawietz, and J. Behler, “High-dimensional neural-network potentials for multicomponent systems: Applications to zinc oxide,” Physical Review B 83, 153101 (2011).
  * Yao  _et al._ (2018) K. Yao, J. E. Herr, D. W. Toth, R. Mckintyre, and J. Parkhill, “The tensormol-0.1 model chemistry: a neural network augmented with long-range physics,” Chemical science 9, 2261–2269 (2018).
  * Plimpton (1995) S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” Journal of Computational Physics 117, 1–19 (1995).
  * Daub, Åstrand, and Bresme (2014) C. D. Daub, P. O. Åstrand, and F. Bresme, “Thermo-molecular orientation effects in fluids of dipolar dumbbells,” Physical Chemistry Chemical Physics 16, 22097–22106 (2014).
  * VandeVondele  _et al._ (2005) J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, “Quickstep: Fast and accurate density functional calculations using a mixed gaussian and plane waves approach,” Computer Physics Communications 167, 103–128 (2005).
  * Zhang and Yang (1998) Y. Zhang and W. Yang, “Comment on “generalized gradient approximation made simple”,” Physical Review Letters 80, 890 (1998).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu,” The Journal of chemical physics 132, 154104 (2010).
  * Goedecker, Teter, and Hutter (1996) S. Goedecker, M. Teter, and J. Hutter, “Separable dual-space gaussian pseudopotentials,” Physical Review B 54, 1703 (1996).
  * Zhang  _et al._ (2018b) L. Zhang, J. Han, H. Wang, R. Car, and E. Weinan, “Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics,” Physical review letters 120, 143001 (2018b).
  * Glielmo, Sollich, and De Vita (2017) A. Glielmo, P. Sollich, and A. De Vita, “Accurate interatomic force fields via machine learning with covariant kernels,” Physical Review B 95, 214302 (2017).
  * Grisafi  _et al._ (2018) A. Grisafi, D. M. Wilkins, G. Csányi, and M. Ceriotti, “Symmetry-adapted machine learning for tensorial properties of atomistic systems,” Physical review letters 120, 036002 (2018).
  * Kingma and Ba (2015) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in  _3rd International Conference on Learning Representations, ICLR 2015_ (2015).
  * Hockney and Eastwood (2021) R. W. Hockney and J. W. Eastwood,  _Computer simulation using particles_ (crc Press, 2021).
  * Berendsen, Grigera, and Straatsma (1987) H. Berendsen, J. Grigera, and T. Straatsma, “The missing term in effective pair potentials,” Journal of Physical Chemistry 91, 6269–6271 (1987).
  * Krogh, Vedelsby  _et al._ (1995) A. Krogh, J. Vedelsby, _et al._ , “Neural network ensembles, cross validation, and active learning,” Advances in neural information processing systems 7, 231–238 (1995).
  * Galib  _et al._ (2017) M. Galib, T. T. Duignan, Y. Misteli, M. D. Baer, G. K. Schenter, J. Hutter, and C. J. Mundy, “Mass density fluctuations in quantum and classical descriptions of liquid water,” The Journal of chemical physics 146, 244501 (2017).
  * Frenkel and Smit (2001) D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ , Vol. 1 (Elsevier, 2001).
  * Niblett and Limmer (2021) S. P. Niblett and D. T. Limmer, “Ion dissociation dynamics in an aqueous premelting layer,” The Journal of Physical Chemistry B 125, 2174–2181 (2021).
  * Schile and Limmer (2019) A. J. Schile and D. T. Limmer, “Rate constants in spatially inhomogeneous systems,” The Journal of chemical physics 150, 191102 (2019).
  * Kattirtzi, Limmer, and Willard (2017) J. A. Kattirtzi, D. T. Limmer, and A. P. Willard, “Microscopic dynamics of charge separation at the aqueous electrochemical interface,” Proceedings of the National Academy of Sciences 114, 13374–13379 (2017).
  * Benjamin (2015) I. Benjamin, “Reaction dynamics at liquid interfaces,” Annual review of physical chemistry 66, 165–188 (2015).
  * Mundy and Kuo (2006) C. J. Mundy and I.-F. W. Kuo, “First-principles approaches to the structure and reactivity of atmospherically relevant aqueous interfaces,” Chemical reviews 106, 1282–1304 (2006).
  * Lewis  _et al._ (2011) T. Lewis, B. Winter, A. C. Stern, M. D. Baer, C. J. Mundy, D. J. Tobias, and J. C. Hemminger, “Does nitric acid dissociate at the aqueous solution surface?” The Journal of Physical Chemistry C 115, 21183–21190 (2011).


# Improving International Climate Policy via Mutually Conditional Binding
Commitments (Track 2)

Jobst Heitzig
FutureLab on Game Theory and Networks of Interacting Agents
Complexity Science Department
Potsdam Institute for Climate Impact Research
<EMAIL_ADDRESS>Jörg Oechssler
Alfred Weber Institute for Economics
University of Heidelberg Christoph Pröschel
Technical University of Berlin
<EMAIL_ADDRESS>Niranjana Ragavan Yat Long Lo
Dyson Robot Learning Lab
<EMAIL_ADDRESS>
###### Abstract

The Paris Agreement, considered a significant milestone in climate
negotiations, has faced challenges in effectively addressing climate change
due to the unconditional nature of most Nationally Determined Contributions
(NDCs). This has resulted in a prevalence of free-riding behavior among major
polluters and a lack of concrete conditionality in NDCs. To address this
issue, we propose the implementation of a decentralized, bottom-up approach
called the Conditional Commitment Mechanism. This mechanism, inspired by the
National Popular Vote Interstate Compact, offers flexibility and incentives
for early adopters, aiming to formalize conditional cooperation in
international climate policy. In this paper, we provide an overview of the
mechanism, its performance in the AI4ClimateCooperation challenge, and discuss
potential real-world implementation aspects. Prior knowledge of the climate
mitigation collective action problem, basic economic principles, and game
theory concepts are assumed.

##  1 Introduction

Even though the Paris Agreement is often framed as a milestone in climate
negotiations, it still proved ineffective [10]. This is in part due to the
unconditionality of most Nationally Determined Contributions (NDCs), leading
to persistent free-riding by large polluters. Only very few NDCs make use of
the conditionality enabled by the Paris Agreement, and mostly only in an
inconcrete or incommensurable and thus largely inconsequential way, even
though the literature has demanded more conditionality [13].

We propose to formalize this type of conditional cooperation using a
decentralized, bottom-up Conditional Commitment Mechanism that was developed
earlier by one of us [3], inspired by the National Popular Vote Interstate
Compact [8], and was shown to have appealing theoretical properties such as
flexibility and incentives for first movers.

Here we summarize the mechanism, its performance in the AI4ClimateCooperation
challenge, and discuss some real-world implementation aspects. We assume
familiarity with the real-world collective action problem of climate
mitigation and some basic knowledge of economic and game-theoretic concepts as
described in the challenge’s white-paper [14].

##  2 Proposed solution: The Conditional Commitment Mechanism

The theoretical basis of our contribution is the Conditional Commitment
Mechanism introduced in [3]: When players in a one-shot “base” game with
(completely or partially) ordered action spaces and preferences that fulfill
certain monotonicity conditions bind themselves via so-called conditional
commitment functions (CCFs) rather than choosing actions directly, then the
outcomes of the strong Nash equilibria [2] of the resulting “commitment” game
are provably exactly the core outcomes of the underlying base game. In
particular, these strategic equilibria are individually rational and Pareto-
optimal. Players who can communicate upfront should thus be able to coordinate
right away on a strong equilibrium that realizes a Pareto-optimal outcome. It
was also shown in [3] that if the base game is repeated, also players who can
not communicate upfront will still converge to Pareto-optimal outcomes when
using a very simple form of non-exploitable better-response learning dynamics.
This suggests that AI agents using more sophisticated reinforcement learning
algorithms should converge to such good equilibria even more reliably.

Mathematically, the CCF $c_{i}$ that player $i$ chooses is a function that
maps each action profile $a_{-i}$ of the other players to a minimal action
$a_{i}=c_{i}(a_{-i})$ that player $i$ commits to realize if the other players
(denoted $-i$) realize at least the actions $a_{-i}$.111Even though it is
formally similar to a best-response function, a CCF does not encode best
responses in the underlying non-cooperative game but can rather be interpreted
as offers for cooperation. It is also somewhat similar to supply function
equilibria in Cournot competition [6]. The mechanism then computes the
(provably existing and unique) largest feasible action profile $a^{\ast}$ so
that $a^{\ast}_{i}\leq c_{i}(a^{\ast}_{-i})$ for all $i$, and all $i$ are then
bound to realize actions that are at least as large as $a^{\ast}_{i}$. 222The
resulting equilibria can be interpreted as a program equilibrium [12] between
a specific class of programs that represent CCFs.

Figure 1: Two-country example of the Conditional Commitment Mechanism. After
each country has submitted a conditional commitment function (CCF), their
liabilities are computed as the the largest feasible action combination.

Since CCFs are “action-type-agnostic”, they are a very flexible tool for
making conditional offers regarding many different forms of actions
(mitigation, consumption, tariffs, subsidies, contributions to funds,
technology transfer, vaccine donations, military aid, whatever) and they can
have any form as long as they are weakly increasing functions. E.g., they
could mix one or more discrete steps with some linearly or nonlinearly
increasing parts, etc., as illustrated in Fig. 1 for a case with just two
countries. For the theoretical results in [3], which generalize earlier
results in [11, 9], already CCFs of a very simple form suffice, namely “two-
step” functions where $c_{i}(a_{-i})=a^{2}_{i}$ iff
$f(a_{-i})\geq\theta^{2}_{i}$, otherwise $c_{i}(a_{-i})=a^{1}_{i}$ iff
$f(a_{-i})\geq\theta^{1}_{i}$, and otherwise $c_{i}(a_{-i})=0$. In this,
$a^{2}_{i}\geq a^{1}_{i}\geq 0$ are two offered action levels for player $i$,
$f$ is some monotonic aggregation function such as the minimum or average of
others’ actions, and $\theta^{2}_{i}\geq\theta^{1}_{i}$ are the “thresholds”
for this aggregate that define the condition for player $i$ to be willing to
realize at least action $a^{2}_{i}$ or $a^{1}_{i}$, respectively.

Motivated by these theoretical results, we designed our conditional commitment
negotiation protocol for the RICE model as follows: In each of RICE’s 5-year
periods: (i) each country (or “region”) $i$ simultaneously submits $k\geq 1$
many pairs $(a^{j}_{i},\theta^{j}_{i})$, $j=1\dots k$. (Typically, $k\leq 4$
suffices for good performance.) These encode a $k$-step CCF $c_{i}$ for this
period and player. Then (ii) the largest feasible profile $a^{\ast}$ is
calculated and is turned into all players’ “liabilities”, and finally (iii)
players choose economic control variables (mitigation, savings, etc.) that are
constrained by these liabilities $a^{\ast}_{i}$.

The components of the “action” $a^{j}_{i}$ are from a subset of the following
economic control variables that are defined or can be derived in the RICE
model: mitigation rate (0…1), consumption rate ($=1-{}$ savings rate), and
“non-tariff on clean imports”. The latter variable is a way to set import
tariffs depending on suppliers’ mitigation rates, representing a form of
carbon border adjustment mechanism (CBAM). More specifically, the “non-tariff
on clean imports” set by country $i$ results in the following constraint on
$i$’s tariff on imports from country $j$: tariff $\leq 1-{}$ “non-tariff on
clean imports”i $\times$ mitigation ratej.

The rationale for this choice of variables in the RICE model is that
mitigation and consumption tend to be costly for $i$ and beneficial for $-i$,
and the opposite for tariffs so that players’ preferences regarding these
action components approximately fulfill the monotonicity requirements of the
CCF theory and thus the players should converge to a strong Nash equilibrium
that realizes an approximately Pareto-optimal outcome.

The thresholds $\theta^{j}_{i}$ are w.r.t. a corresponding subset of the
aggregate variables global mitigation rate, global consumption rate, global
“non-tariff on clean imports”, which are suitably weighted global averages.
The rationale is that conditioning on global aggregates is less complex (which
should accelerate learning), more flexible to fulfill (which increases the
potential for cooperation), and more realistic than conditioning on individual
players’ actions.

For example, if $k=2$, a player might offer (i) to mitigate $\geq 20\%$,
consume $\geq 80\%$, and tax clean imports at rate $\leq 30\%$ in the coming
period if global mitigation is $\geq 30\%$, global consumption is $\geq 70\%$,
and globally avg. tariffs on clean imports are $\leq 40\%$, and (ii) to even
mitigate $\geq 40\%$, consume $\geq 90\%$, and tax clean imports at rate $\leq
20\%$ if global mitigation is $\geq 40\%$, global consumption is $\geq 80\%$,
and globally avg. tariffs on clean imports are $\leq 30\%$. Depending on
whether and how this matches with other cpuntries CCFs, the player’s actual
economic control variables are then constrained by either (i) or (ii) or
remain unconstrained for this period.

A possible real-world implementation of the CCF mechanism would be via
unilateral legislation whereby an individual country passes a bill that (i)
specifies a CCF encoding offers and conditions, and (ii) bindingly commits
that country to comply with the largest action profile that is feasible (in
the above sense) given all CCFs that are specified in similar bills that are
put into force by other countries.

For example, country A might start by passing a bill that (i) offers to match
every 2 GtC of global mitigation by $\geq 1$ GtC of domestic mitigation, up to
at most 10 GtC/year, and (ii) commits country A to this action as soon as
other countries have passed similar bills matching the conditions set in (i).
If then, say, countries B, C and D pass similar bills offering 8 GtC each if
global mitigation is $\geq 30$ GtC, then the action profile $(10,8,8,8)$ GtC
for (A,B,C,D) matches all four countries’ conditions and is thus immediately
becoming a binding liability without the countries having to sign any treaty
or pass further legislation. B could later amend their bill to include a
further offer to, say, contribute at least USD 10 bln to a global compensation
fund if country E reduces their tariffs on certain products to $\leq 10\%$. As
soon as E passes a matching bill, both actions become additional liabilities.
This way, countries can step-wise increase their ambition regarding various
interrelated variables with positive externalities for other countries in a
bottom-up emerging process not requiring multilateral treaties or upfront
negotiations between many governments in the usual sense.

##  3 Effectiveness

In theory, rational countries should easily realize that our mechanism allows
them to cooperate effectively from the very beginning to achieve Pareto-
optimal, individually and collectively rational outcomes. They can get to this
insight either by classical game-theoretical inspection of the mechanism or by
reinforcement learning in a simulated environment. In particular, our
mechanism does not require a slow ramping-up of ambitions or a complex
intertemporal policy involving reciprocity. To verify this, we performed
several forms of simulation experiments. Already in [3], we verified the
convergence to the core in a simple best-response dynamics in a three-player
public goods game with one control variable (mitigation), linear benefits, and
quadratic costs.

For this challenge, we repeated that exercise in the same simple environment
with the neural-network-based reinforcement learning algorithm supplied by the
competition organizers, restricted to simple “$k$-step-function” CCFs (see
above). This showed convergence to fairly good outcomes significantly better
than in the case without negotiations, but they did not constitute strong Nash
equilibria (which can be derived easily analytically in this simple game, see
Appendix). We suspect this is due to the insufficient capabilities of the
prescribed learning algorithms (Appendix).

We noted that due to exploration, agents sometimes converge to action
combinations above the Pareto-frontier, i.e., mitigating more than optimal. To
avoid this, we implemented an additional variant of the mechanism in which the
ordering of the $k$ offer-condition-pairs $(a^{j}_{i},\theta^{j}_{i})$ play a
prioritizing role, so that additional offers will not always increase ambition
but might also reduce ambition to optimal levels if prioritized correctly. Two
further variations of the mechanism, as well as a temporary reward shaping
trick during learning were introduced to help coordination towards a
beneficial equilibrium, which in our simulations indeed tended to raise the
resulting level of ambition (see Appendix for details).

Finally, we used the same shipped RL algorithm to assess performance in the
RICE-N framework using several control variables. As of the date of
submission, our results in this environment are promising but mixed and
somewhat inconclusive, but suggest similar issues with the RL algorithm as
already identified in the simpler environment, which is why we choose not to
use them as evidence in this track 2 submission even though we have submitted
them to track 1.

##  4 Feasibility

Real-world feasibility will mainly depend on the degree to which countries can
enforce the necessary actions within their borders, the legal feasibility of
passing the respective legislation that implements the mechanism in individual
countries, which depends on the political systems of the respective countries,
and on the level of bindingness of the resulting bills. If the latter is low,
the resulting compliance problem might have to be addressed using suitable
monitoring systems and additional sanctioning mechanisms, as in any form of
treaty. In many countries, citizens or NGOs could also enforce compliance with
passed bills by appealing to supreme courts. If the used CCFs are smooth
enough (rather than step-functions as in the RL simulations), noncompliance
could be deterred in a similar way as in ref. [5]. In our simulations, we
assumed full compliance.

Since in contrast to other mechanisms, ours does not require upfront
multilateral negotiations between governments (even though it might profit
from them) but relies purely on unilateral legislation that can be passed by
parliaments, it may result in an empowerment of parliaments relative to
governments and in a much larger transparency than negotiations behind closed
doors, thereby also increasing democratic legitimacy. It might also reduce the
frictions and delays of the typical two-level game of negotiation and
ratification occurring in traditional treaties [7].

Given the salience of these benefits, this should increase incentives and
possibilities for non-governmental actors (such as members of parliament,
NGOs, citizens) to push for CCF legislation. They can push for ambitious
climate policy more easily when proposing conditional rather than
unconditional action. Early moving domestic actors can argue that their
country can provide an example and set the agenda and a norm by encoding
certain concrete actions and conditions into its CCF. If early movers use
linear CCFs, thereby offering to match others’ ambitions with own
contributions, then domestic actors in other countries can point out that now
those countries’ actions would have much larger leverage if they also pass a
corresponding CCF. This is similar to how charitable fundraisers sometimes
incentivize donations by matching them with additional funds. These later
movers also have incentives to use metrics in their CCFs commensurable to
early movers’ metric in order to make their bill consequential.

##  5 Robustness

Besides the above-mentioned possible issues with non-compliance if countries
cannot bindingly commit, a general failure mode that we have not addressed is
irrational behavior such as not accepting individually rational offers. It is
also unclear how bills passed under the CCF mechanism will be constrained by
or influence other types of treaties that exist in parallel, such as bottom-up
carbon market linkages (see, e.g., [4]). A coalition of countries that have
established liabilities via our mechanism can however easily turn them into a
formal multilateral treaty later to increase stability further, albeit losing
some flexibility.

Due to its flexibility, our mechanism can be used in a way that is robust
against fluctuations in mitigation costs, economic performance etc.: If
countries use smooth CCFs rather than step functions, then a country that,
say, experiences a recession or unexpectedly large mitigation costs can
unilaterally change their CCF and all other countries’ liabilities would
automatically be adjusted accordingly, without them having to amend their
laws.

##  6 Ethics and climate justice

The suggested mechanism is procedurally fair in a formal sense since it treats
all countries the same. There is no special role for large actors. However,
countries with a faster legislative process might gain an advantage by passing
their bills earlier than others, thereby setting the agenda. Our theoretical
results suggest that countries might converge to any outcome in the core of
the base game, which might be more or less equitable. The fact that they also
have to be individually rational for all countries can be seen as a pragmatic
constraint on the level of justice that can be achieved, which is however hard
to avoid in any form of free negotiations.

To push in the direction of outcomes associated with calls for Climate
Justice, early movers should condition their CCFs on aggregate variables that
take into account historical responsibility. Non-country actors such as the
UN, the IPCC, or NGOs, may influence this towards equity by proposing suitable
aggregate metrics to condition on.

Finally, countries that already act unconditionally may be disadvantaged
because they have less to offer further, which is however also a common
feature of negotiations.

##  7 Novelty

The novelty of the CCF mechanism in the climate context is the switch from
multilateral coordinated negotiations to unilateral binding but mutually
conditional commitments, adding flexibility and speed.

##  8 Conclusion

The proposed mechanism combines the flexibility of the Paris agreement with
the bindingness and concreteness of traditional treaties. Between rational
countries, it should converge fast to a Pareto-optimal, individually and
collectively rational outcome. Next steps are the design of robust legal
language for domestic bills implementing the mechanism by key players within
varying legal frameworks, the identification of specific unilateral offers and
conditions to be included in a first wave of corresponding unilateral bills,
and talks with relevant policy-makers that could advance these bills.

### Author contributions

J.H. designed the mechanism and performed the game-theoretical analysis. J.H.
and C.P. implemented the mechanism in the challenge framework. C.P. performed
Monte-Carlo simulations of the learning process. N.R. visualized the results
and did research on continuous action spaces. Y.L.L. consulted about
reinforcement learning methodology. J.O. consulted about economic theory.
J.H., C.P. and N.R. wrote and edited the paper.

### Carbon accounting

We roughly estimated our computations to have caused carbon emissions on the
order of 1.5 metric tons of CO2 equivalents.

## References

  * [1] Maruan Al-Shedivat, Trapit Bansal, Yuri Burda, Ilya Sutskever, Igor Mordatch, and Pieter Abbeel.  Continuous adaptation via meta-learning in nonstationary and competitive environments.  arXiv preprint arXiv:1710.03641, 2017.
  * [2] Robert J Aumann.  Acceptable points in general cooperative n-person games.  Contributions to the Theory of Games, 4:287–324, 1959.
  * [3] Jobst Heitzig.  Efficient non-cooperative provision of costly positive externalities via conditional commitments.  Available at SSRN 3449004, 2019.
  * [4] Jobst Heitzig and Ulrike Kornek.  Bottom-up linking of carbon markets under far-sighted cap coordination and reversibility.  Nature Climate Change, 8(3):204–209, 2018.
  * [5] Jobst Heitzig, Kai Lessmann, and Yong Zou.  Self-enforcing strategies to deter free-riding in the climate change mitigation game and other repeated public good games.  Proceedings of the National Academy of Sciences, 108(38):15739–15744, 2011.
  * [6] Paul D Klemperer and Margaret A Meyer.  Supply function equilibria in oligopoly under uncertainty.  Econometrica: Journal of the Econometric Society, pages 1243–1277, 1989.
  * [7] Michael Lisowski.  Playing the two-level game: Us president bush’s decision to repudiate the kyoto protocol.  Environmental Politics, 11(4):101–119, 2002.
  * [8] Derek T Muller.  The compact clause and the national popular vote interstate compact.  Election Law Journal, 6(4):372–393, 2007.
  * [9] Joerg Oechssler, Andreas Reischmann, and Andis Sofianos.  The conditional contribution mechanism for repeated public goods–the general case.  Journal of Economic Theory, 203:105488, 2022.
  * [10] Kilian Raiser, Ulrike Kornek, Christian Flachsland, and William F Lamb.  Is the paris agreement effective? a systematic map of the evidence.  Environmental Research Letters, 15(8):083006, 2020.
  * [11] Andreas Reischmann and Joerg Oechssler.  The binary conditional contribution mechanism for public good provision in dynamic settings—theory and experimental evidence.  Journal of Public Economics, 159:104–115, 2018.
  * [12] Moshe Tennenholtz.  Program equilibrium.  Games and Economic Behavior, 49(2):363–373, 2004.
  * [13] David G Victor.  The case for climate clubs.  by International Centre for Trade and Sustainable Development (ICTSD) 7 …, 2015.
  * [14] Tianyu Zhang, Andrew Williams, Soham Phade, Sunil Srinivasa, Yang Zhang, Prateek Gupta, Yoshua Bengio, and Stephan Zheng.  Ai for global climate cooperation: Modeling global climate negotiations, agreements, and long-term cooperation in rice-n, 2022.
  * [15] Zihao Zhang, Stefan Zohren, and Stephen Roberts.  Deep reinforcement learning for trading.  The Journal of Financial Data Science, 2(2):25–40, 2020.

## Appendix

### Simple additional public goods game

In addition to RICE-N, we also performed simulations with a very simple public
goods game to assess the RL algorithm’s performance. Each of the $N$ regions
$i$ has constant capital and production, both normalized to 1. There is no
growth or capital accumulation or depreciation, and there is no trade (savings
rates and trade-related control variables are thus ignored). The only
meaningful control variable is thus the regional mitigation rate
$\mu_{i}\in[0,1]$, and regional emissions are then $1-\mu_{i}$ per period. We
assume quadratic mitigation costs of $\gamma\mu_{i}^{2}$ which are subtracted
from production, where $\gamma>0$ is a cost parameter. The global temperature
anomaly $T$ simply equals global emissions here, and there is no independent
temperature dynamics. Regional climate damages are linear, $2\beta
T=2\beta(N-\sum_{j}\mu_{j})$, so that regional period payoffs are
$\Pi_{i}=1-\gamma\mu_{i}^{2}-2\beta N+2\beta\sum_{j}\mu_{j}$. In the base
game’s Nash equilibrium (i.e., without negotiations), all players choose
$\mu_{i}\equiv\mu^{NE}=\beta/\gamma$, giving regional period payoffs of
$\Pi^{NE}=1-2\beta N+(2N-1)\beta^{2}/\gamma$. In contrast, the socially
optimal choice is $\mu_{i}\equiv\mu^{\ast}=N\beta/\gamma$, giving higher
regional period payoffs of $\Pi^{\ast}=1-2\beta N+N^{2}\beta^{2}/\gamma$.

When putting $N=2$, $\beta=0.09$ and $\gamma=0.3$, we get
$\mu^{NE}=0.3<\mu^{\ast}=0.6$ and $N\Pi^{NE}=1.442<N\Pi^{\ast}=1.496$. With
the conditional commitment mechanism, both players can easily achieve the
social optimum by putting $(a_{i},\theta_{i})\equiv(0.6,0.6)$, where
$a_{i}=\mu_{i}$ and $\theta_{i}=\sum_{j}\mu_{j}/N$, which is a Pareto-optimal
and individually and collectively rational Nash equilibrium of the mechanism.
Indeed, with the simple learning dynamics we had earlier performed in [3],
players converge fast to this outcome.

Still, at least in our simulations, the shipped RL algorithms from this
competition were unfortunately unable to learn the optimal equilibrium even in
this almost trivial two-player setup and rather ended up with global period
payoffs around 1.47, i.e., about halfway between the disagreement and socially
optimal outcomes.

### Variants of the mechanism to improve learnability by RL agents

#### Prioritized offers

In the basic version of the mechanism, the pairs $(a^{j}_{i},\theta^{j}_{i})$
are just used to define a $k$-step CCF and the mechanism selects the largest
action profile feasible in all CCFs. Once that is above the Pareto-frontier,
players can only reduce it again to optimal levels by removing offers and
replacing them by lower ones. In order to avoid falling back to zero
liabilities when removing too large offers, they would first need to add a
lower offer before removing the larger one. But since adding a lower offer
makes no immediate difference but only has a delayed effect as soon as the
larger offer is removed, this is hard to learn for RL agents.

We thus studied a variant of the mechanism in which the ordering of the
submitted points $(a^{j}_{i},\theta^{j}_{i})$ plays a prioritizing role. More
precisely, the mechanism would first construct all CCFs as before and find
those combinations of submitted points $(a^{j}_{i},\theta^{j}_{i})$ that are
feasible in all CCFs. Rather than selecting the largest such, it would then
for each player $i$ find that feasible profile which is the most prioritized
by that player. It would then select the supremum of all those “favourite
feasible” profiles. So if all players prioritize a point on the Pareto-
frontier over the other possible agreements above the Pareto-frontier, the
former would be selected, allowing players to get down to the Pareto-frontier
more easily when inadvertently being above it. For rational (rather than
learning) agents, this makes no difference, the strong equilibria are the same
as before, they can only be learned more easily now.

#### Borda-scoring

While the prioritizing approach improved performance somewhat in the simple
public good game with few players, it did not do so when the number of players
is too large, because even a single player’s wrong prioritization can drive
the supremum of the favourite feasible profiles above the Pareto-frontier. We
therefore added a further modification in which taking the supremum as the
last step is replaced by performing a simple form of automatic voting
mechanism using the submitted prioritization. More precisely, we used the
Borda scoring method for this.

#### Option to adjust towards the mean

Another issue for convergence is the non-uniqueness of Pareto-efficient,
individually and collectively rational possible agreement points in action
space and the resulting equilibrium selection and coordination problem. To
address this, we finally studied a variant of the mechanism that has an
additional step in each period: After all countries have submitted their
prioritized list of points $(a^{j}_{i},\theta^{j}_{i})$ for this period, each
country $i$ can choose to replace their highest-priority point
$(a^{1}_{i},\theta^{1}_{i})$ by a certain form of average of all other
players’ highest-priority points $(a^{1}_{j},\theta^{1}_{j})$. If all
countries choose to do so, they would thereby agree on this average as the
outcome.

### Suggested improvements in RL algorithms

It was somewhat unclear to us whether the RL algorithms supplied to
participants (based on the A2C algorithm) can be expected to converge to a
Nash equilibrium (NE) or some refinement of NE (such as subgame-perfect or
strong equilibrium), or at least to a low-regret state of the negotiation game
in the first place, and if so to which particular equilibrium if there are
many (as is usually the case in complex games like this).

From comparing the used RL approach with established alternative approaches to
find equilibria used in environmental economics (e.g., optimization with
welfare weights using gradient descent in action space), we suspect that the
used RL agents might not be reliably able to find Pareto-efficient and
collectively rational equilibria for at least two independent reasons.

#### No discretization of continuous action space

For one, the RL agents use a randomized policy over a discretized action space
rather than a deterministic policy over the actual continuous action space.
Learning via gradient descent thus cannot exploit the topological structure of
the action space. Rather than shifting actions smoothly around in action
space, they shift probability mass from points in action space to possibly
quite far away regions in action space. But given the largely approximately
convex public-good structure of the relevant utility functions, it is not to
be expected that Pareto-efficient NEs are mixed rather than pure strategies,
and that a classical gradient descent in action space (rather than in
probability space) should converge well. we hence suggest dropping the action
discretization and use algorithms such as the continuous-action version of A2C
[15].

#### Meta-learning MARL and anticipation

Independently of this, the fact that the learning algorithm for the RL agent
representing country $i$ is treating the RL agents representing the other
countries as part of the environment makes the environment non-stationary
since all RL agents are trained simultaneously and thus have changing
policies. Since convergence theorems usually have to assume stationarity, one
can expect weak convergence because of this non-stationarity. This is not
surprising given the longstanding insights from the early theory of learning
in games that showed that already simple best-response dynamics can lead to
oscillations rather than convergence since each player is “chasing a moving
target”. To overcome this, we suggest using state-of-the-art Multi-Agent
Reinforcement Learning (MARL) algorithms that treat the environment and other
agents differently, such as the “meta-learning” approach in [1].

We believe that particularly designs in which an RL agent anticipates the
other RL agents’ next learning steps would be very conducive to convergence in
a negotiation context. This is because if $i$ switches in a learning step from
making a bad offer to making a good offer, this will only pay off after the
other players have switched in their next learning step from not accepting the
offer to accepting the offer.

#### Reward shaping via temporary altruism

As we were not supposed to change the RL algorithms used in the competition,
the only modification of the RL algorithm that we actually tested in our
simulations was a certain temporary modification of the reward function in the
hope of helping the simultaneously learning agents converge to good rather
than bad Nash equilibria. The idea is to replace the selfish reward function
$r_{i}=\Pi_{i}$ (regional period payoff) by a partially altruistic reward
function $r_{i}^{\prime}=(1-\alpha)r_{i}+\alpha\sum_{j}r_{j}/N$ with an
altruism parameter $\alpha\in[0,1]$. If we set $\alpha\approx 1$ in the early
phase of learning, all agents try to maximize approximately the same objective
function, so there is no conflict between the learners and they should
consequently converge to this maximum. As this maximum is on the Pareto-
frontier of the actual payoff space, it is either within or close to that part
of the Pareto-frontier that is individually and collectively rational, i.e. in
the core. If in the later phase of learning $\alpha$ is then slowly reduced to
zero, our expectation was that the learners would then move from the found
point on the Pareto frontier to a Nash equilibrium of the mechanism that is
not too far away from it and still on the Pareto frontier. Since all points in
the core are indeed Nash equilibria, the hope was that thus the learners would
converge to the core rather than to worse equilibria.

Despite this reasoning, this modification did not show a significant
improvement of final reward in our experiments in which we tested several
different “schedules” for $\alpha$, all having $\alpha\equiv 0$ during the
last third of the learning process. This however can also be due to the
possible general issues with convergence to Nash equilibrium that we suspected
(see above). We, therefore, believe that this reward shaping by temporary
altruism during learning should be revisited with more suitable RL algorithms.