# Dielectric Saturation in Water from a Long Range Machine Learning Model

Harender S. Dhattarwal  Department of Chemistry and Chemical Biology, Rutgers
University, Piscataway, NJ 08854  Ang Gao  Department of Physics, Beijing
University of Posts and Telecommunications, 100876, Beijing, China  Richard C.
Remsing<EMAIL_ADDRESS>Department of Chemistry and Chemical
Biology, Rutgers University, Piscataway, NJ 08854

###### Abstract

Machine learning-based neural network potentials have the ability to provide
ab initio-level predictions while reaching large length and time scales often
limited to empirical force fields. Traditionally, neural network potentials
rely on a local description of atomic environments to achieve this
scalability. These local descriptions result in short range models that
neglect long range interactions necessary for processes like dielectric
screening in polar liquids. Several approaches to including long range
electrostatic interactions within neural network models have appeared
recently, and here we investigate the transferability of one such model, the
self consistent neural network (SCFNN), which focuses on learning the physics
associated with long range response. By learning the essential physics, one
can expect that such a neural network model should exhibit at least partial
transferability. We illustrate this transferability by modeling dielectric
saturation in a SCFNN model of water. We show that the SCFNN model can predict
non-linear response at high electric fields, including saturation of the
dielectric constant, without training the model on these high field strengths
and the resulting liquid configurations. We then use these simulations to
examine the nuclear and electronic structure changes underlying dielectric
saturation. Our results suggest that neural network models can exhibit
transferability beyond the linear response regime and make genuine predictions
when the relevant physics is properly learned.

##  I Introduction

Ab initio molecular dynamics simulations enable computational predictions of
interatomic interactions and chemical reactivity Car and Parrinello (1985);
Sun  _et al._ (2016); Chen  _et al._ (2017); Tuckerman  _et al._ (1996); Marx
and Hutter (2009); Kresse and Hafner (1993); Szabo and Ostlund (2012);
Dhattarwal  _et al._ (2022); Groß and Sakong (2022), but the expense of
performing electronic structure calculations limits their use to small system
sizes and short time scales. Empirical interaction potentials, or force
fields, are regularly used to model large systems and long time scales Karplus
and Petsko (1990); Hansson  _et al._ (2002); Riniker (2018); Venable  _et al._
(2019); Brooks  _et al._ (2021), but it is difficult to include processes like
chemical reactions and electronic polarization in these classical models Leven
_et al._ (2021); Shi  _et al._ (2021); Jing  _et al._ (2019); Bedrov  _et al._
(2019). To bridge the gap between ab initio and force field-based simulations,
machine learning-based neural network models are being developed to achieve ab
initio accuracy at a fraction of the computational cost Behler and Parrinello
(2007); Zhang  _et al._ (2018); Hansen  _et al._ (2015); Artrith  _et al._
(2011).

Neural network models learn interatomic interactions from ab initio
calculations to enable efficient sampling of potential energy surfaces Blank
_et al._ (1995); Rowe  _et al._ (2020); Noé  _et al._ (2019); Behler (2016);
Yao  _et al._ (2022). By learning ab initio-level interactions, the resulting
neural network models are able to predict electronic polarization effects and
bond breakage and formation in large systems and on larger timescales than
accessible in current ab initio simulations Behler (2017); Xu  _et al._
(2021); van Gerwen  _et al._ (2022); Vandermause  _et al._ (2022); Cools-
Ceuppens  _et al._ (2022); Manzhos and Carrington (2021); Dhabal  _et al._
(2022); Gartner  _et al._ (2020). To reach these large scales, neural network
models often operate under the assumption of locality, in which interatomic
interactions are determined by atomic arrangements within a spherical region
typically on the scale of 1 nm or less Behler (2016); Behler and Parrinello
(2007); Behler (2021). These methods have created profound insights into
chemical systems and revolutionized molecular simulations Behler (2016);
Manzhos and Carrington (2021); Kapil  _et al._ (2022); Gartner III  _et al._
(2022). However, this assumption of locality results in short range models
that lack a description of long range electrostatics Yue  _et al._ (2021);
Niblett  _et al._ (2021); Behler (2021); Gao and Remsing (2022). As a result,
several approaches to modeling long range interactions in neural network
models have begun to appear Artrith  _et al._ (2011); Grisafi and Ceriotti
(2019); Ko  _et al._ (2021); Zhang  _et al._ (2022); Yao  _et al._ (2018); Gao
and Remsing (2022); Pagotto  _et al._ (2022).

The recently-developed self-consistent field neural network (SCFNN) separately
learns short and long range interactions through two coupled neural networks
Gao and Remsing (2022). This enables the SCFNN to learn long range response,
while including the impact of long range effects on short range structure and
interactions through a rapidly converging self-consistent loop. The resulting
SCFNN model can accurately describe dielectric screening and the response of a
system to external electrostatic fields. Importantly, the SCFNN uses a
physically-meaningful separation of interactions into short and long range
interactions Rodgers and Weeks (2008); Remsing  _et al._ (2016); Gao  _et al._
(2020) and focuses on learning the physics underlying the response to long
range fields. Consequently, the SCFNN is partially transferable to
environments for which it was not trained.

Training the SCFNN involves computing energies and forces for configurations
of a system, as well as the energies and forces of those same nuclear
configurations in the presence of electric fields. Despite this training being
performed on configurations and fields in the linear response regime, the
SCFNN model makes no assumptions of linear response. As a result, the SCFNN
should be able to describe non-linear response. Here, we illustrate the
transferability of the SCFNN model to the non-linear response regime by
modeling dielectric saturation in liquid water.

When an external electric field is applied to a polar liquid, the liquid will
respond to screen the field through changes in electronic and nuclear
structure Onsager (1936); Kirkwood (1939); Fröhlich (1949); Böttcher (1973);
Zangwill (2013); Chandler (1977); De Leeuw  _et al._ (1986); Yeh and Berkowitz
(1999); Zhang  _et al._ (2020); Seyedi  _et al._ (2019). This response
consists mainly of changes in the orientational nuclear structure of the
liquid, such that molecules reorient their dipoles in the direction of the
field. For small fields, the liquid’s response is linear and determined by its
static dielectric constant Onsager (1936); Kirkwood (1939); Fröhlich (1949);
Böttcher (1973); Purcell and Morin (2013); Zangwill (2013); Chandler (1977);
Hansen and McDonald (2006); Cox and Sprik (2019); Cox (2020). However, as the
field is increased in magnitude, the liquid does not continue to respond
linearly; there is a limit to the amount of dipolar reorientation that can
occur. As the maximum value of dipolar reorientation is approached, the liquid
responds non-linearly to the external electric field Debye (1929); Booth
(1951, 1955); Sutmann (1998); Alper and Levy (1990); Sandberg and Edholm
(2002); Zhang and Sprik (2016); Yeh and Berkowitz (1999); Matyushov (2015,
2018); Richert (2017); Kotodziej  _et al._ (1975); Jeanmairet  _et al._
(2019); Willard  _et al._ (2009). Consequently, dipolar fluctuations are
damped, and there is a reduction in the dielectric constant — dielectric
saturation.

The non-linear response to external electric fields underlying dielectric
saturation can be used as a test on the transferability of the SCFNN model.
Here, we show that the SCFNN model can describe dielectric saturation in water
without including configurations consistent with this effect in the training
data. After demonstrating that the SCFNN is transferable to the non-linear
response regime, we examine classic theories and quantify the structural
underpinnings of dielectric saturation, including electronic polarization
effects. We conclude with a discussion of our results in the context of neural
network models.

##  II Simulation Details

We modeled the dielectric response of water to the applied homogeneous
electrostatic fields using molecular dynamics simulations. Eight different
electric fields, 0.005, 0.01, 0.015, 0.02, 0.05, 0.10, 0.20, and 0.28 V/Å,
were applied to a cubic box of 1000 water molecules with dimension 31.2 Å. A
total of 56 independent trajectories were sampled for the three fields less
than 0.02 V/Å, and a total of 28 independent trajectories were sampled for the
rest. Each system was equilibrated for 50 ps at 300 K in the NVT ensemble,
with the temperature maintained using the Berendsen thermostat Berendsen  _et
al._ (1984). The equations of motion were integrated using a time step of 0.5
fs. The systems were periodic in all three directions. The last 25 ps of all
trajectories were used for analysis.

The SCFNN was trained following previous work Gao and Remsing (2022). Training
the SCFNN includes computing energies and forces for configurations of water
with fixed nuclear positions in the absence and presence of electric fields.
The original SCFNN model of water used electric field of magnitude 0.1 and 0.2
V/Å to learn the long range electronic response of water to applied fields Gao
and Remsing (2022). We refer to this model as SCFNN(HF) to indicate that it
was trained at high field strengths. In this study, we trained another SCFNN
model using smaller electric fields of magnitude 0.005 and 0.01 VÅ, which
corresponds to the linear polarization regime for both electronic and nuclear
response, based on previous studies Zhang and Sprik (2016); Cox and Sprik
(2019). We refer to this as the SCFNN model throughout this work. The network
architecture of both models is the same. We use a Behler-Parrinello style
network for the short range part of the SCFNN Behler and Parrinello (2007);
Morawietz  _et al._ (2016); Gao and Remsing (2022), though we expect that many
established approaches could readily be used for the short system.

The DFT calculations for training followed previous work Cheng  _et al._
(2019); Marsalek and Markland (2017); Gao and Remsing (2022) and used the CP2K
program VandeVondele  _et al._ (2005); Kühne  _et al._ (2020) with the revised
Perdew-Burke-Ernzerhof hybrid exchange-correlation functional with 25% exact
exchange (revPBE0) Adamo and Barone (1999); Zhang and Yang (1998); Goerigk and
Grimme (2011) and k-point sampling at the $\Gamma$-point only. Goedecker-
Teter-Hutter (GTH) pseudopotentials Goedecker  _et al._ (1996) were used with
TZV2P basis sets VandeVondele and Hutter (2007). DFT-D3 dispersion corrections
were used to include long ranged van der Waals interactions Grimme  _et al._
(2010). The calculations used an energy cutoff of 400 Ry and 60 Ry for the
reference grid (keyword REL_CUTOFF). Maximally localized Wannier function
centers were obtained by minimizing the spread of MLWFs within CP2K Marzari
_et al._ (2012); Berghold  _et al._ (2000). Configurations for the training
data were taken from previous work Cheng  _et al._ (2019).

Figure 1:  (a) Polarization as a function of electric field strength for the
SCFNN and SPC/E water models; SPC/E data is from previous work Zhang and Sprik
(2016). Also shown is the polarization computed using the structure predicted
by the SCFNN model but with the SPC/E dipole moment, SCFNN($\mu_{\rm SPC/E}$),
and with the value of the SCFNN dipole moment at zero field, SCFNN($\mu_{0}$).
(b) Comparison of the polarization as a function of electric field strength
for the SCFNN and SCFNN(HF) models. Error bars indicate one standard
deviation.

##  III Results and Discussion

To examine the nonlinear response of water to applied electric fields, we
apply uniform fields of magnitude E to water and compute the polarization
induced in the system according to Cox and Sprik (2019)

$\centering
P(E)=\left<P(\overline{\textbf{R}})\right>_{E}=\left<\frac{1}{V}\sum_{i=1}^{N}\mu_{i}(\overline{\textbf{R}})\right>_{E}.\@add@centering$
(1)

Here, $\left<\cdots\right>_{E}$ indicates an ensemble average over
configurations $\overline{\textbf{R}}$ in the presence of a field of strength
E, such that $P(\overline{\textbf{R}})$ is the instantaneous polarization of a
single configuration $\overline{\textbf{R}}$, $V$ is the volume of the
simulation cell containing $N$ water molecules, and
$\mu_{i}(\overline{\textbf{R}})$ is the instantaneous dipole moment of
molecule $i$. The polarization is linear at low field strength, but becomes
non-linear beyond approximately 0.02 V/Å, Fig. 1. For large enough fields,
$P(E)$ begins to plateau, indicative of dielectric saturation Alper and Levy
(1990); Yeh and Berkowitz (1999); Sutmann (1998); Sandberg and Edholm (2002);
Cox (2020); Zhang and Sprik (2016). This non-linear behavior is well-described
by the SCFNN model despite it not being trained in this regime.

We also compare the response of the SCFNN to that of the SPC/E water model
obtained in previous work Zhang and Sprik (2016). The SPC/E model displays
similar behavior although its polarization in the nonlinear regime is
significantly smaller than that of the SCFNN model. Classic treatments of
dielectric response treat water as a collection of independent dipoles, such
that the orientation and magnitude of the molecular dipole moment determines
the polarization. From this perspective, the SPC/E and SCFNN models should
have different polarization values because their average dipole moments
differ: 2.9 D for SCFNN in zero field and 2.35 D for SPC/E. To further this
comparison, we computed the polarization of SCFNN configurations assuming that
the magnitude of every dipole is the same as that of SPC/E. The resulting
polarization, indicated by SCFNN($\mu_{\rm SPC/E}$) in Fig. 1a, is very
similar to that of the SPC/E model. This agreement suggests that the change in
the orientation of molecular dipoles induced by the field is similar in the
SCFNN and SPC/E models.

Training the SCFNN model used here used lower field strengths than the
original SCFNN model Gao and Remsing (2022). To ensure that the model is
robust, we performed the same simulations using the original model, termed
SCFNN(HF) to indicate that the fields used in training this model, 0.1 and 0.2
V/Å, were higher than those used to train what we refer to as the SCFNN model,
0.005 and 0.01 V/Å. The polarization predicted by both models is essentially
identical, Fig. 1b. This excellent agreement further suggests that the SCFNN
model is learning the underlying physics responsible for screening and, as a
result, this long range part of the neural network model is transferable. The
remaining results are shown only for SCFNN unless indicated otherwise; the
SCFNN(HF) yields essentially identical results.

We can use our predictions of the polarization to estimate the field-dependent
dielectric constant, $\varepsilon(E)$, from the derivative of the polarization
with respect to the field, according to Fröhlich (1949); Böttcher (1973);
Purcell and Morin (2013); Zangwill (2013); Kusalik (1994); Zhang and Sprik
(2016)

$\centering\varepsilon(E)=1+\frac{4\pi P(E)}{E}.\@add@centering$ (2)

As field strength is increased, $\varepsilon(E)$ monotonically decreases and
begins to plateau at high field strengths (Fig. 2). This is a clear indication
of dielectric saturation at high fields. The dielectric constant of the SCFNN
model is larger than that of SPC/E but follows the same general trend.

Figure 2:  The electric field-dependence of water’s dielectric constant,
$\varepsilon(E)$. Dashed lines correspond to predictions from the Kirkwood-
Booth Equation 3, using $\varepsilon_{\infty}=1.65$ for SCFNN and
$\varepsilon_{\infty}=1$ for SPC/E. Error bars indicate one standard
deviation.

At high field strengths, the dielectric constant can be described by the
Kirkwood-Booth equation Kirkwood (1939); Fröhlich (1948); Booth (1951)

$\centering\varepsilon(E)=\varepsilon_{\infty}+\frac{7N(\varepsilon_{\infty}+2)}{3V\epsilon_{0}\sqrt{73}}\frac{\mu_{0}}{E}\mathcal{L}\left(\frac{\sqrt{73}(\varepsilon_{\infty}+2)}{6k_{\rm
B}T}\frac{\mu_{0}}{E}\right),\@add@centering$ (3)

where $\varepsilon_{\infty}$ is the high frequency dielectric constant,

$\mu_{E}=\left<\frac{1}{N}\sum_{i=1}^{N}\left|\mu_{i}(\overline{\textbf{R}})\right|\right>_{E}$
(4)

is the average magnitude of the water dipole in the presence of an electric
field of magnitude $E$, such that $\mu_{0}=\mu_{E=0}$, and
$\mathcal{L}(x)=\coth(x)-1/x$ is the Langevin function Langevin (1905). The
Kirkwood-Booth equation neglects solvent-solvent correlations beyond the first
coordination shell and assumes that the first coordination shell is
independent of field strength and described by the Bernal-Fowler model of
water Bernal and Fowler (1933). Despite these simplifications, Eq. 3 was shown
to provide a reasonable description of dielectric saturation. In part, the
accuracy of the model was achieved by adjusting the dipole moment of water
from that in the gas phase to $\mu_{0}=2.1$ D, and fitting the value of the
dipole moment in this way compensates for errors arising from the above
approximations Booth (1951). We find that the Kirkwood-Booth model with
$\mu_{0}$ determined from simulations qualitatively produces dielectric
saturation, is quantitatively accurate in the high field regime for the SCFNN
model, and describes the dielectric constant nearly everywhere for the SPC/E
model, Fig. 2.

Figure 3:  Total dipole moment of the system relative to the maximum value,
equivalent to $\left<\cos\theta_{z}\right>$, where $\theta_{z}$ is the angle
between the molecular dipole moment vector and the direction of the field.
Dashed and dotted lines correspond to predictions from the Langevin function,
Eq. 5, with $\alpha_{0}$ and $\alpha_{E}$, respectively. Error bars indicate
one standard deviation.

In a similar manner, the polarization can be estimated using the Langevin
function Debye (1929); Booth (1951); Langevin (1905),

$\centering\frac{\left<M_{z}\right>_{E}}{N\mu_{E}}=\left<\cos\theta_{z}\right>_{E}\approx\mathcal{L}\left(\alpha_{0}\right),\@add@centering$
(5)

where $\left<M_{z}\right>_{E}$ is the total dipole moment of the system in the
presence of the field E, $\theta_{z}$ is the angle made by the water dipole
moment vector and the direction of the field, $\alpha_{E}=3\mu_{E}/2k_{\rm
B}T$, and $k_{\rm B}T$ is the product of Boltzmann’s constant and the
temperature. Equation 5 arises from an independent dipole approximation,
wherein each dipolar molecule is embedded in a dielectric medium and
intermolecular correlations between the dipoles are ignored Booth (1951);
Langevin (1905). While the Langevin function produces the general trend, it
underestimates the value of $\left<\cos\theta_{z}\right>_{E}$ predicted by the
simulations for all but the smallest and largest fields. Moreover, the
agreement with the SCFNN results is deceptive because the Langevin function
assumes that the magnitude of the dipole moment is independent of field
strength, but the dipole moment can change in the SCFNN model. The good
agreement between SCFNN and the Langevin prediction results from scaling
$\left<M_{z}\right>_{E}$ by $\mu_{E}$ and not $\mu_{0}$; scaling by $\mu_{0}$
can result in values larger than one Piekara (1962). As a result, one may
anticipate that polarization in the SCFNN model includes a non-zero electronic
contribution.

Figure 4:  (a) Average molecular dipole moment of water as a function of
electric field strength. The dashed line shows a linear fit to the initial
increase between $0.02$ and $0.1$ V/Å. (b) The probability distribution of the
water dipole moment for all field strengths studied here. The darkness of the
lines are proportional to the field strength. Figure 5:  Simulation snapshots
illustrating the dipoles (blue lines) of water molecules (red oxygen and white
hydrogen) under constant electric fields of strength (a) 0.005 and (b) 0.28
V/Å, respectively.

Figure 6:  The (a-c) radial distribution functions, $g(r)$, and (d) total
X-ray scattering structure factor, $S(q)$, quantifying interatomic
correlations in water. The darkness of the lines are proportional to the field
strength.

We examine the average molecular dipole moment of water as a function of
field, $\mu_{E}$ (Fig. 4). At low fields, in the linear regime, the average
molecular dipole moment is nearly constant. At high field strengths, the
average molecular dipole moment increases with field strength, as anticipated
from the discussion above. Unlike the polarization, the increase in the dipole
moment does not saturate at high E, although it does increase nonlinearly at
the highest field values studied here. This suggests that the electronic
polarization of water is beginning to saturate as well.

The nature of the fluctuations in water’s molecular dipole moment are
characterized by its probability distribution for varying electric field
strengths (Fig. 4b). We find that the dipole moment fluctuations remain
Gaussian for all fields studied, evidenced by the constant width of the dipole
moment probability distributions. At low fields, in the linear response
regime, the dipole moment distributions are constant with increasing field,
consistent with the constant dipole moment at low fields (Fig. 4a). At high
fields, the only change in the distributions is a shift of their mean,
consistent with the non-linear increase in the dipole moment. Interestingly,
the dipole moment continues to increase at the highest fields studied, despite
the structural changes saturating at high fields. This suggests that water can
continue to respond to an applied field through electronic polarization even
when the nuclear structure can no longer produce an increase in polarization.

To quantify the contribution of electronic polarization to the total
polarization of water at high fields, we computed the polarization from SCFNN
with the magnitude of every dipole moment replaced by the average value at
zero field, $\mu_{E}\approx\mu_{0}$. At low fields, the polarization produced
by this approximation is the same as that computed in the simulations (Fig.
1). However, the polarization obtained in this approximation begins to
underestimate the actual polarization at high fields. The difference between
these two curves reflects the contribution of electronic polarization to the
total, which arises from changes in the magnitude of the molecular dipole
moment of water.

One can attempt to modify traditional theories like the Langevin function to
account for the field dependence of $\mu_{E}$. Inserting $\mu_{E}$ into the
Langevin function of Eq. 5,
$\left<\cos\theta_{z}\right>_{E}\approx\mathcal{L}(\alpha_{E})$, results in
better agreement with the SCFNN model, as shown in Fig. 3, where we estimated
$\mu_{E}$ for all fields by fitting the simulation results to a third order
polynomial. The increase in the dipole moment at higher fields increases the
estimate of $\left<\cos\theta_{z}\right>_{E}$, so much so that the SCFNN
results are now slightly overestimated. This overestimation likely arises from
the neglect of correlations in the Langevin function (the independent dipole
approximation). In water, the hydrogen bond network places constraints on the
orientations that molecules can adopt, and so water molecules cannot align
with the applied field to the same extent as independent dipoles, resulting in
a lower $\left<\cos\theta_{z}\right>_{E}$ Piekara (1962). While using
$\mu_{E}$ instead of $\mu_{0}$ may improve the predictions of Eq. 5, making
the same substitution in Eq. 3 results in minimal changes in its estimate of
the dielectric constant.

The saturation of the polarization at high field strengths occurs when water
molecules maximally orient their dipoles along the field direction, as shown
in Fig. 5. At low field strengths, the dipolar structure of the liquid is
disordered (Fig. 5a). This is consistent with water exhibiting a linear
response at low fields, because the liquid structure is essentially the same
as that at zero field (see below for more details). In contrast, at high field
strength, water dipoles preferentially align in one direction (Fig. 5b), and
the structure of water differs from that at zero field. Once this large
polarization value is reached, additional polarization can be achieved through
electronic polarization. Note that this type of polarization is absent in
rigid, fixed point charge models like SPC/E.

We first quantify the changes in water structure induced by uniform electric
fields through the site-site pair radial distribution functions (RDFs),
$g(r)$, Fig. 6. In the linear regime, the RDFs remain unchanged as the field
is increased. At high fields, the most significant changes are found in the
O-O RDF, in which the peaks sharpen and the minima deepen. There are small
changes in the third and fourth peaks of the H-H RDF, and the O-H RDF is
essentially unchanged. Importantly, the first, intramolecular peaks in the H-H
and O-H RDFs remain unchanged, suggesting that the bond lengths and angles are
not altered by the field strengths studied here and the increase in the dipole
moment with field strength arises from electronic polarization. This
electronic response is captured in the SCFNN by its ability to describe the
long range interactions between the applied field and the molecular charge
distribution.

Because the changes in the pair structure at high field occur mainly at large
distances, they may be better quantified through X-ray scattering structure
functions, $S(q)$, shown in Fig. 6d. In the linear regime, $S(q)$ is unchanged
as the field strength is increased, as found for the RDFs. In the nonlinear
regime, however, the first two peaks in $S(q)$ increase and the first minimum
decreases as the field strength is increased. The first peak consists mainly
of O-H and O-O correlations, while the second peak is dominated by O-O
correlations Skinner  _et al._ (2014), such that the changes in $S(q)$ are
consistent with O-O correlations changing most significantly at high fields.
At all field strengths, the intramolecular peaks in $S(q)$ at high $q$ are
unchanged, further suggesting that the molecular geometry is unchanged as high
fields are applied. These changes in $S(q)$ clearly indicate an increase in
intermolecular ordering at high field strengths.

Figure 7:  (a) Cosine of the angle between the dipole moment vectors of two
water molecules separated by a distance $r$ at different applied electrostatic
fields. The darkness of the lines are proportional to the field strength. The
vertical dashed line indicates the location of the first peak in $g_{\rm
OO}(r)$, which is used to define the nearest neighbor distance $r_{\rm NN}$.
(b) Electric field dependence of the cosine of the angle between the dipole
moment vectors of two water molecules separated by a nearest neighbor
distance, $r_{\rm NN}$, defined by the location of the first peak in the
oxygen-oxygen RDF, and two water molecules separated by large distances,
$r=\infty$.

Although the RDFs and $S(q)$ show some increase in translational ordering of
water, the snapshots in Fig. 5 suggests that the changes in water structure
are mainly in orientational ordering. To quantify the changes in orientational
structure, we compute $\left<\cos\theta_{\rm w}(r)\right>_{E}$, where
$\theta_{\rm w}(r)$ is the angle between the dipole moment vectors of two
water molecules separated by a distance $r$, shown in Fig. 7a. The
orientational structure of water, as quantified by $\left<\cos\theta_{\rm
w}(r)\right>_{E}$ changes significantly with the strength of the applied
electric field. At small fields, $\left<\cos\theta_{\rm w}(r)\right>_{E}$
displays a large peak near 2.5 Å that corresponds to dipolar ordering of
neighboring water molecules. After this peak, $\left<\cos\theta_{\rm
w}(r)\right>_{E}$ exhibits small oscillations and tends to small values at
large $r$; $\left<\cos\theta_{\rm w}(r)\right>_{E}$ goes to zero at large
distances in the absence of an applied field.

As E increases, the first peak in $\left<\cos\theta_{\rm w}(r)\right>_{E}$
grows until the polarization begins to saturate, at which point the peak
height changes very little with field. In addition, $\left<\cos\theta_{\rm
w}(r)\right>_{E}$ for large $r$ plateaus at higher values with large fields.
This indicates long range ordering of water molecules in the presence of the
electric field. Moreover, at the highest fields studied,
$\left<\cos\theta_{\rm w}(r)\right>_{E}$ is comparable in magnitude at all
distances, suggesting that water is approaching maximal orientational order.
This highly ordered state is responsible for the plateau in the polarization
and the saturation of the dielectric constant at high fields.

The behavior of $\left<\cos\theta_{\rm w}(r)\right>_{E}$ sheds light on the
assumptions made in the Kirkwood-Booth Equation 3. At low field strength,
$\left<\cos\theta_{\rm w}(r)\right>_{E}$ exhibits a large peak in the first
coordination shell, but is small after that, lending support to Kirkwood’s
model for dipolar order in water Kirkwood (1939); Fröhlich (1948); Booth
(1951). Furthermore, the value of $\left<\cos\theta_{\rm w}(r)\right>_{E}$ at
a typical nearest neighbor distance, $r_{\rm NN}$, is close to the value of
1/3 used in Kirkwood’s model Kirkwood (1939); Bernal and Fowler (1933), where
we have defined $r_{\rm NN}$ by the location of the first peak in the O-O RDF.
However, as the field strength increases into the non-linear regime, the
correlations between nearest neighbors increase and those beyond the first
shell become significant with increasing field strength, Fig. 7b. These field-
dependent changes in pair structure are ignored in the simplistic Kirkwood-
Booth model, and inclusion of these effects may improve its accuracy.

##  IV Conclusion

We have applied the SCFNN model of water to dielectric saturation in water.
For large uniform electric fields with magnitudes of approximately 0.05 V/Å or
larger, water responds non-linearly to the applied field and the induced
polarization starts to plateau. This plateau arises from maximal reorientation
of water molecules. The increased alignment of water dipoles restricts their
fluctuations and consequently lowers the dielectric constant. This dielectric
saturation can be described reasonably well with the Kirkwood-Booth theory
Booth (1951); Kirkwood (1939); Fröhlich (1948), and we have examined some of
the key assumptions in this model. Despite the saturation of the orientational
structure of water, the electronic structure continues to respond at the
highest field strengths studied here, evidenced by the dipole moment
increasing in magnitude, albeit non-linearly, at high fields. Despite the non-
linear increase of molecular dipole moments due to electronic polarization,
fluctuations of the molecular dipoles remain Gaussian at all fields studied
here.

The SCFNN model was not trained on configurations or fields typical of
dielectric saturation. Instead, non-linear response emerges naturally within
the SCFNN framework because the model learns the long range response
responsible for dielectric screening and how this response impacts the short
range structure and interactions in water. Our results highlight that the
SCFNN can be transferable, something that is beyond the reach of many machine
learning-based models. We expect that the transferability of the current SCFNN
model is largely limited to the long range response, while the short range
interactions will need to be retrained when additional local interactions are
introduced, such as those between water and ions. However, we anticipate that
the SCFNN idea of focusing on physical origins and length scales of molecular
interactions, and appropriately adapting the resulting neural network
structure, could enable the development of more transferable neural network
models.

###### Acknowledgements.

We acknowledge the Office of Advanced Research Computing (OARC) at Rutgers,
The State University of New Jersey for providing access to the Amarel cluster
and associated research computing resources. This work used Anvil at Purdue
through allocation CHE210081 from the Advanced Cyberinfrastructure
Coordination Ecosystem: Services & Support (ACCESS) program, which is
supported by National Science Foundation grant number 2005632.

## References

  * Car and Parrinello (1985) R. Car and M. Parrinello, “Unified approach for molecular dynamics and density-functional theory,” Phys. Rev. Lett. 55, 2471–2474 (1985).
  * Sun  _et al._ (2016) Jianwei Sun, Richard C. Remsing, Yubo Zhang, Zhaoru Sun, Adrienn Ruzsinszky, Haowei Peng, Zenghui Yang, Arpita Paul, Umesh Waghmare, Xifan Wu, Michael L. Klein, and John P. Perdew, “Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional,” Nat. Chem. 8, 831–836 (2016).
  * Chen  _et al._ (2017) Mohan Chen, Hsin-Yu Ko, Richard C. Remsing, Marcos F. Calegari Andrade, Biswajit Santra, Zhaoru Sun, Annabella Selloni, Roberto Car, Michael L. Klein, John P. Perdew, and Xifan Wu, “Ab initio theory and modeling of water,” Proceedings of the National Academy of Sciences 114, 10846–10851 (2017).
  * Tuckerman  _et al._ (1996) Mark E. Tuckerman, P. Jeffrey Ungar, Tycho von Rosenvinge, and Michael L. Klein, “Ab initio molecular dynamics simulations,” The Journal of Physical Chemistry 100, 12878–12887 (1996).
  * Marx and Hutter (2009) Dominik Marx and Jürg Hutter,  _Ab initio molecular dynamics: basic theory and advanced methods_ (Cambridge University Press, 2009).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, “Ab initio molecular dynamics for liquid metals,” Phys. Rev. B 47, 558–561 (1993).
  * Szabo and Ostlund (2012) Attila Szabo and Neil S Ostlund,  _Modern quantum chemistry: introduction to advanced electronic structure theory_ (Courier Corporation, 2012).
  * Dhattarwal  _et al._ (2022) Harender S. Dhattarwal, Jer-Lai Kuo, and Hemant K. Kashyap, “Mechanistic insight on the stability of ether and fluorinated ether solvent-based lithium bis(fluoromethanesulfonyl) electrolytes near li metal surface,” The Journal of Physical Chemistry C 126, 8953–8963 (2022).
  * Groß and Sakong (2022) Axel Groß and Sung Sakong, “Ab initio simulations of water/metal interfaces,” Chemical Reviews 122, 10746–10776 (2022).
  * Karplus and Petsko (1990) Martin Karplus and Gregory A. Petsko, “Molecular dynamics simulations in biology,” Nature 347, 631–639 (1990).
  * Hansson  _et al._ (2002) Tomas Hansson, Chris Oostenbrink, and WilfredF van Gunsteren, “Molecular dynamics simulations,” Current Opinion in Structural Biology 12, 190–196 (2002).
  * Riniker (2018) Sereina Riniker, “Fixed-charge atomistic force fields for molecular dynamics simulations in the condensed phase: An overview,” Journal of Chemical Information and Modeling 58, 565–578 (2018).
  * Venable  _et al._ (2019) Richard M. Venable, Andreas Krämer, and Richard W. Pastor, “Molecular dynamics simulations of membrane permeability,” Chemical Reviews 119, 5954–5997 (2019).
  * Brooks  _et al._ (2021) Charles L. Brooks, David A. Case, Steve Plimpton, Benoît Roux, David van der Spoel, and Emad Tajkhorshid, “Classical molecular dynamics,” The Journal of Chemical Physics 154, 100401 (2021).
  * Leven  _et al._ (2021) Itai Leven, Hongxia Hao, Songchen Tan, Xingyi Guan, Katheryn A. Penrod, Dooman Akbarian, Benjamin Evangelisti, Md Jamil Hossain, Md Mahbubul Islam, Jason P. Koski, Stan Moore, Hasan Metin Aktulga, Adri C. T. van Duin, and Teresa Head-Gordon, “Recent advances for improving the accuracy, transferability, and efficiency of reactive force fields,” Journal of Chemical Theory and Computation 17, 3237–3251 (2021).
  * Shi  _et al._ (2021) Zhihao Shi, Jian Zhou, and Runjie Li, “Application of reaction force field molecular dynamics in lithium batteries,” Frontiers in Chemistry 8 (2021).
  * Jing  _et al._ (2019) Zhifeng Jing, Chengwen Liu, Sara Y. Cheng, Rui Qi, Brandon D. Walker, Jean-Philip Piquemal, and Pengyu Ren, “Polarizable force fields for biomolecular simulations: Recent advances and applications,” Annual Review of Biophysics 48, 371–394 (2019).
  * Bedrov  _et al._ (2019) Dmitry Bedrov, Jean-Philip Piquemal, Oleg Borodin, Alexander D. Jr. MacKerell, Benoît Roux, and Christian Schröder, “Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields,” Chemical Reviews 119, 7940–7995 (2019).
  * Behler and Parrinello (2007) Jörg Behler and Michele Parrinello, “Generalized neural-network representation of high-dimensional potential-energy surfaces,” Phys. Rev. Lett. 98, 146401 (2007).
  * Zhang  _et al._ (2018) Linfeng Zhang, Jiequn Han, Han Wang, Roberto Car, and Weinan E, “Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics,” Phys. Rev. Lett. 120, 143001 (2018).
  * Hansen  _et al._ (2015) Katja Hansen, Franziska Biegler, Raghunathan Ramakrishnan, Wiktor Pronobis, O. Anatole von Lilienfeld, Klaus-Robert Müller, and Alexandre Tkatchenko, “Machine learning predictions of molecular properties: Accurate many-body potentials and nonlocality in chemical space,” The Journal of Physical Chemistry Letters 6, 2326–2331 (2015).
  * Artrith  _et al._ (2011) Nongnuch Artrith, Tobias Morawietz, and Jörg Behler, “High-dimensional neural-network potentials for multicomponent systems: Applications to zinc oxide,” Phys. Rev. B 83, 153101 (2011).
  * Blank  _et al._ (1995) Thomas B. Blank, Steven D. Brown, August W. Calhoun, and Douglas J. Doren, “Neural network models of potential energy surfaces,” The Journal of Chemical Physics 103, 4129–4137 (1995).
  * Rowe  _et al._ (2020) Patrick Rowe, Volker L. Deringer, Piero Gasparotto, Gábor Csányi, and Angelos Michaelides, “An accurate and transferable machine learning potential for carbon,” The Journal of Chemical Physics 153, 034702 (2020).
  * Noé  _et al._ (2019) Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu, “Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning,” Science 365, eaaw1147 (2019).
  * Behler (2016) Jörg Behler, “Perspective: Machine learning potentials for atomistic simulations,” The Journal of Chemical Physics 145, 170901 (2016).
  * Yao  _et al._ (2022) Nan Yao, Xiang Chen, Zhong-Heng Fu, and Qiang Zhang, “Applying classical, ab initio, and machine-learning molecular dynamics simulations to the liquid electrolyte for rechargeable batteries,” Chemical Reviews 122, 10970–11021 (2022).
  * Behler (2017) Jörg Behler, “First principles neural network potentials for reactive simulations of large molecular and condensed systems,” Angewandte Chemie International Edition 56, 12828–12840 (2017).
  * Xu  _et al._ (2021) Jiayan Xu, Xiao-Ming Cao, and P. Hu, “Perspective on computational reaction prediction using machine learning methods in heterogeneous catalysis,” Phys. Chem. Chem. Phys. 23, 11155–11179 (2021).
  * van Gerwen  _et al._ (2022) Puck van Gerwen, Alberto Fabrizio, Matthew D Wodrich, and Clemence Corminboeuf, “Physics-based representations for machine learning properties of chemical reactions,” Machine Learning: Science and Technology 3, 045005 (2022).
  * Vandermause  _et al._ (2022) Jonathan Vandermause, Yu Xie, Jin Soo Lim, Cameron J. Owen, and Boris Kozinsky, “Active learning of reactive bayesian force fields applied to heterogeneous catalysis dynamics of h/pt,” Nature Communications 13, 5183 (2022).
  * Cools-Ceuppens  _et al._ (2022) Maarten Cools-Ceuppens, Joni Dambre, and Toon Verstraelen, “Modeling electronic response properties with an explicit-electron machine learning potential,” Journal of Chemical Theory and Computation 18, 1672–1691 (2022).
  * Manzhos and Carrington (2021) Sergei Manzhos and Tucker Jr. Carrington, “Neural network potential energy surfaces for small molecules and reactions,” Chemical Reviews 121, 10187–10217 (2021).
  * Dhabal  _et al._ (2022) Debdas Dhabal, Subramanian K. R. S. Sankaranarayanan, and Valeria Molinero, “Stability and metastability of liquid water in a machine-learned coarse-grained model with short-range interactions,” The Journal of Physical Chemistry B 126, 9881–9892 (2022).
  * Gartner  _et al._ (2020) Thomas E. Gartner, Linfeng Zhang, Pablo M. Piaggi, Roberto Car, Athanassios Z. Panagiotopoulos, and Pablo G. Debenedetti, “Signatures of a liquid–liquid transition in an ab initio deep neural network model for water,” Proceedings of the National Academy of Sciences 117, 26040–26046 (2020).
  * Behler (2021) Jörg Behler, “Four generations of high-dimensional neural network potentials,” Chemical Reviews 121, 10037–10072 (2021).
  * Kapil  _et al._ (2022) Venkat Kapil, Christoph Schran, Andrea Zen, Ji Chen, Chris J Pickard, and Angelos Michaelides, “The first-principles phase diagram of monolayer nanoconfined water,” Nature 609, 512–516 (2022).
  * Gartner III  _et al._ (2022) Thomas E Gartner III, Pablo M Piaggi, Roberto Car, Athanassios Z Panagiotopoulos, and Pablo G Debenedetti, “Liquid-liquid transition in water from first principles,” Physical Review Letters 129, 255702 (2022).
  * Yue  _et al._ (2021) Shuwen Yue, Maria Carolina Muniz, Marcos F. Calegari Andrade, Linfeng Zhang, Roberto Car, and Athanassios Z. Panagiotopoulos, “When do short-range atomistic machine-learning models fall short?” The Journal of Chemical Physics 154, 034111 (2021).
  * Niblett  _et al._ (2021) Samuel P. Niblett, Mirza Galib, and David T. Limmer, “Learning intermolecular forces at liquid–vapor interfaces,” The Journal of Chemical Physics 155, 164101 (2021).
  * Gao and Remsing (2022) Ang Gao and Richard C. Remsing, “Self-consistent determination of long-range electrostatics in neural network potentials,” Nature Communications 13, 1572 (2022).
  * Grisafi and Ceriotti (2019) Andrea Grisafi and Michele Ceriotti, “Incorporating long-range physics in atomic-scale machine learning,” The Journal of Chemical Physics 151, 204105 (2019).
  * Ko  _et al._ (2021) Tsz Wai Ko, Jonas A. Finkler, Stefan Goedecker, and Jörg Behler, “A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer,” Nature Communications 12, 398 (2021).
  * Zhang  _et al._ (2022) Linfeng Zhang, Han Wang, Maria Carolina Muniz, Athanassios Z. Panagiotopoulos, Roberto Car, and Weinan E, “A deep potential model with long-range electrostatic interactions,” The Journal of Chemical Physics 156, 124107 (2022).
  * Yao  _et al._ (2018) Kun Yao, John E. Herr, David W. Toth, Ryker Mckintyre, and John Parkhill, “The tensormol-0.1 model chemistry: a neural network augmented with long-range physics,” Chem. Sci. 9, 2261–2269 (2018).
  * Pagotto  _et al._ (2022) Joshua Pagotto, Junji Zhang, and Timothy Duignan, “Predicting the properties of salt water using neural network potentials and continuum solvent theory,” (2022).
  * Rodgers and Weeks (2008) Jocelyn M Rodgers and John D Weeks, “Local molecular field theory for the treatment of electrostatics,” Journal of Physics: Condensed Matter 20, 494206 (2008).
  * Remsing  _et al._ (2016) Richard C. Remsing, Shule Liu, and John D. Weeks, “Long-ranged contributions to solvation free energies from theory and short-ranged models,” Proceedings of the National Academy of Sciences 113, 2819–2826 (2016).
  * Gao  _et al._ (2020) Ang Gao, Richard C. Remsing, and John D. Weeks, “Short solvent model for ion correlations and hydrophobic association,” Proc. Natl. Acad. Sci. U.S.A. 117, 1293–1302 (2020).
  * Onsager (1936) Lars Onsager, “Electric moments of molecules in liquids,” Journal of the American Chemical Society 58, 1486–1493 (1936).
  * Kirkwood (1939) John G Kirkwood, “The dielectric polarization of polar liquids,” The Journal of Chemical Physics 7, 911–919 (1939).
  * Fröhlich (1949) Herbert Fröhlich,  _Theory of dielectrics: dielectric constant and dielectric loss_ (Clarendon Press, 1949).
  * Böttcher (1973) C.J.F. Böttcher,  _Theory of Electric Polarization: Dielectrics in Static Fields_ , 2nd ed., Vol. 1 (Elsevier, Amsterdam, 1973).
  * Zangwill (2013) A. Zangwill,  _Modern Electrodynamics_ (Cambridge University Press, 2013).
  * Chandler (1977) David Chandler, “The dielectric constant and related equilibrium properties of molecular fluids: Interaction site cluster theory analysis,” The Journal of Chemical Physics 67, 1113–1124 (1977).
  * De Leeuw  _et al._ (1986) SW De Leeuw, John W Perram, and ER Smith, “Computer simulation of the static dielectric constant of systems with permanent electric dipoles,” Annual review of physical chemistry 37, 245–270 (1986).
  * Yeh and Berkowitz (1999) In-Chul Yeh and Max L. Berkowitz, “Dielectric constant of water at high electric fields: Molecular dynamics study,” The Journal of Chemical Physics 110, 7935–7942 (1999).
  * Zhang  _et al._ (2020) Chao Zhang, Thomas Sayer, Jürg Hutter, and Michiel Sprik, “Modelling electrochemical systems with finite field molecular dynamics,” Journal of Physics: Energy 2, 032005 (2020).
  * Seyedi  _et al._ (2019) Salman Seyedi, Daniel R Martin, and Dmitry V Matyushov, “Screening of coulomb interactions in liquid dielectrics,” Journal of Physics: Condensed Matter 31, 325101 (2019).
  * Purcell and Morin (2013) Edward M Purcell and David J Morin,  _Electricity and magnetism_ (Cambridge University Press, 2013).
  * Hansen and McDonald (2006) J. P. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ (Elsevier Ltd., 2006).
  * Cox and Sprik (2019) Stephen J. Cox and Michiel Sprik, “Finite field formalism for bulk electrolyte solutions,” The Journal of Chemical Physics 151, 064506 (2019).
  * Cox (2020) Stephen J. Cox, “Dielectric response with short-ranged electrostatics,” Proceedings of the National Academy of Sciences 117, 19746–19752 (2020).
  * Debye (1929) P Debye, “Polar molecules, the chemical catalog company,” Inc., New York , 77–108 (1929).
  * Booth (1951) F. Booth, “The dielectric constant of water and the saturation effect,” The Journal of Chemical Physics 19, 391–394 (1951).
  * Booth (1955) F Booth, “Dielectric constant of polar liquids at high field strengths,” The Journal of Chemical Physics 23, 453–457 (1955).
  * Sutmann (1998) Godehard Sutmann, “Structure formation and dynamics of water in strong external electric fields,” Journal of Electroanalytical Chemistry 450, 289–302 (1998).
  * Alper and Levy (1990) Howard E. Alper and Ronald M. Levy, “Field strength dependence of dielectric saturation in liquid water,” The Journal of Physical Chemistry 94, 8401–8403 (1990).
  * Sandberg and Edholm (2002) Lars Sandberg and Olle Edholm, “Nonlinear response effects in continuum models of the hydration of ions,” The Journal of Chemical Physics 116, 2936–2944 (2002).
  * Zhang and Sprik (2016) Chao Zhang and Michiel Sprik, “Computing the dielectric constant of liquid water at constant dielectric displacement,” Phys. Rev. B 93, 144201 (2016).
  * Matyushov (2015) Dmitry V. Matyushov, “Nonlinear dielectric response of polar liquids,” The Journal of Chemical Physics 142, 244502 (2015).
  * Matyushov (2018) Dmitry V Matyushov, “Nonlinear dielectric response of polar liquids,” in  _Nonlinear Dielectric Spectroscopy_ (Springer, 2018) pp. 1–34.
  * Richert (2017) Ranko Richert, “Nonlinear dielectric effects in liquids: a guided tour,” Journal of Physics: Condensed Matter 29, 363001 (2017).
  * Kotodziej  _et al._ (1975) HA Kotodziej, G Parry Jones, and Mansel Davies, “High field dielectric measurements in water,” Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics 71, 269–274 (1975).
  * Jeanmairet  _et al._ (2019) Guillaume Jeanmairet, Benjamin Rotenberg, Daniel Borgis, and Mathieu Salanne, “Study of a water-graphene capacitor with molecular density functional theory,” The Journal of Chemical Physics 151, 124111 (2019).
  * Willard  _et al._ (2009) Adam P Willard, Stewart K Reed, Paul A Madden, and David Chandler, “Water at an electrochemical interface—a simulation study,” Faraday discussions 141, 423–441 (2009).
  * Berendsen  _et al._ (1984) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, “Molecular dynamics with coupling to an external bath,” The Journal of Chemical Physics 81, 3684–3690 (1984).
  * Morawietz  _et al._ (2016) Tobias Morawietz, Andreas Singraber, Christoph Dellago, and Jörg Behler, “How van der waals interactions determine the unique properties of water,” Proceedings of the National Academy of Sciences 113, 8368–8373 (2016).
  * Cheng  _et al._ (2019) Bingqing Cheng, Edgar A. Engel, Jörg Behler, Christoph Dellago, and Michele Ceriotti, “Ab initio thermodynamics of liquid and solid water,” Proceedings of the National Academy of Sciences 116, 1110–1115 (2019).
  * Marsalek and Markland (2017) Ondrej Marsalek and Thomas E Markland, “Quantum dynamics and spectroscopy of ab initio liquid water: The interplay of nuclear and electronic quantum effects,” The journal of physical chemistry letters 8, 1545–1551 (2017).
  * VandeVondele  _et al._ (2005) Joost VandeVondele, Matthias Krack, Fawzi Mohamed, Michele Parrinello, Thomas Chassaing, and Jürg Hutter, “Quickstep: Fast and accurate density functional calculations using a mixed gaussian and plane waves approach,” Computer Physics Communications 167, 103–128 (2005).
  * Kühne  _et al._ (2020) Thomas D. Kühne, Marcella Iannuzzi, Mauro Del Ben, Vladimir V. Rybkin, Patrick Seewald, Frederick Stein, Teodoro Laino, Rustam Z. Khaliullin, Ole Schütt, Florian Schiffmann, Dorothea Golze, Jan Wilhelm, Sergey Chulkov, Mohammad Hossein Bani-Hashemian, Valéry Weber, Urban Borštnik, Mathieu Taillefumier, Alice Shoshana Jakobovits, Alfio Lazzaro, Hans Pabst, Tiziano Müller, Robert Schade, Manuel Guidon, Samuel Andermatt, Nico Holmberg, Gregory K. Schenter, Anna Hehn, Augustin Bussy, Fabian Belleflamme, Gloria Tabacchi, Andreas Glöß, Michael Lass, Iain Bethune, Christopher J. Mundy, Christian Plessl, Matt Watkins, Joost VandeVondele, Matthias Krack, and Jürg Hutter, “Cp2k: An electronic structure and molecular dynamics software package - quickstep: Efficient and accurate electronic structure calculations,” The Journal of Chemical Physics 152, 194103 (2020).
  * Adamo and Barone (1999) Carlo Adamo and Vincenzo Barone, “Toward reliable density functional methods without adjustable parameters: The pbe0 model,” The Journal of Chemical Physics 110, 6158–6170 (1999).
  * Zhang and Yang (1998) Yingkai Zhang and Weitao Yang, “Comment on “generalized gradient approximation made simple”,” Phys. Rev. Lett. 80, 890–890 (1998).
  * Goerigk and Grimme (2011) Lars Goerigk and Stefan Grimme, “A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions,” Phys. Chem. Chem. Phys. 13, 6670–6688 (2011).
  * Goedecker  _et al._ (1996) S. Goedecker, M. Teter, and J. Hutter, “Separable dual-space gaussian pseudopotentials,” Phys. Rev. B 54, 1703–1710 (1996).
  * VandeVondele and Hutter (2007) Joost VandeVondele and Jürg Hutter, “Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases,” The Journal of Chemical Physics 127, 114105 (2007).
  * Grimme  _et al._ (2010) Stefan Grimme, Jens Antony, Stephan Ehrlich, and Helge 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).
  * Marzari  _et al._ (2012) Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, Ivo Souza, and David Vanderbilt, “Maximally localized wannier functions: Theory and applications,” Rev. Mod. Phys. 84, 1419–1475 (2012).
  * Berghold  _et al._ (2000) Gerd Berghold, Christopher J. Mundy, Aldo H. Romero, Jürg Hutter, and Michele Parrinello, “General and efficient algorithms for obtaining maximally localized wannier functions,” Phys. Rev. B 61, 10040–10048 (2000).
  * Kusalik (1994) P.G. Kusalik, “Computer simulation study of a highly polar fluid under the influence of static electric fields,” Molecular Physics 81, 199–216 (1994).
  * Fröhlich (1948) H Fröhlich, “General theory of the static dielectric constant,” Transactions of the Faraday Society 44, 238–243 (1948).
  * Langevin (1905) Paul Langevin, “Magnétisme et théorie des électrons,” Ann. chim. et phys. , 203 (1905).
  * Bernal and Fowler (1933) John D Bernal and Ralph H Fowler, “A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions,” The Journal of Chemical Physics 1, 515–548 (1933).
  * Piekara (1962) A Piekara, “Dielectric saturation and hydrogen bonding,” The Journal of Chemical Physics 36, 2145–2150 (1962).
  * Skinner  _et al._ (2014) L. B. Skinner, C. J. Benmore, J. C. Neuefeind, and J. B. Parise, “The structure of water around the compressibility minimum,” The Journal of Chemical Physics 141, 214507 (2014).