Supporting Information for “Stabilization of AgI’s polar surfaces by the aqueous
environment, and its implications for ice formation”
Thomas Sayer1 and Stephen J. Cox1, a)
Department of Chemistry, University of Cambridge, Lensfield Road,
Cambridge CB2 1EW, United Kingdom
(Dated: 12 June 2019)
In this supporting information, we show that while using a ‘mirrored slab’ geometry
may remove the net dipole of the simulation cell, it does not provide polarity com-
pensation. A derivation of the CNC conditions is presented in the case that R1 6= R2,
which is the relevant case for AgI. We also present results not included in the main
article for various combinations of fields and slab mobility, as well as results using
different force fields and cell sizes. Finally, we give values of the parameters used in
nonelectrostatic interactions.
a)Electronic mail: sjc236@cam.ac.uk
1
I. MIRRORED SLAB GEOMETRY: REMOVING THE NET DIPOLE VS.
POLARITY COMPENSATION
In the main article, we compared the behavior of AgI in contact with pure H2O at
Dz = 0 against simulations performed at DCNC and ECNC. Recall that the Hamiltonian
for Dz = 0 is equivalent to the commonly used Yeh-Berkowitz correction,
1 which adds a
force Fz,i = −4piqiPz to each particle i in the simulation cell. In those simulations, the
system comprised a single AgI slab immersed in H2O. We found a stark contrast in the
behavior of the system under CNC conditions versus at Dz = 0, with a large potential drop
of ∆xtlφ = −46.2 V across the crystal slab in the case of the latter. As the force applied by
the YB correction is proportional to the net dipole of the simulation cell, another approach
that has been taken is to use a ‘mirrored slab’ geometry, in which two slabs are used in the
simulation, positioned in such a way that their dipoles cancel. If the net polarization of the
simulation cell vanishes on average, and standard tin-foil boundary conditions are used, then
we argue that this is equivalent on average to the Dz = 0 ensemble. This is because tin foil
boundary conditions correspond to the Ez = 0 ensemble, giving 〈Dz〉Ez = Ez+4pi〈P 〉Ez = 0.
We therefore do not expect the mirrored slab geometry to enforce CNC conditions.
To test the above hypothesis, we have performed a simulation with Lz = 21.75 nm. An
immobile AgI crystal is placed with its lowermost Ag-(0 0 0 1) face at z = −7.42 nm and its
uppermost I-(0 0 0 1) face at z = −4.33 nm. A second crystal is then placed such that its
dipole points in the opposite direction, with its lowermost I-(0 0 0 1) face at z = +4.33 nm
and its uppermost Ag-(0 0 0 1) face at z = +7.42 nm. In the region −4.33 < z < +4.33 nm
we placed 750 TIP4P/2005 water molecules. The remaining volume was occupied by vacuum.
All other simulation settings are the same as outlined in the main article. The resulting
electrostatic potential profile is shown in Fig. S1 (a). We can see that the potential drop
across each slab closely resembles that seen for the simulations presented in the main article
at Dz = 0. Indeed, we find |∆xtlφ| ≈ 46.2 V/nm in excellent agreement with the Dz = 0
result presented in Fig. 3. We have also performed the simulations with the polarity of
the AgI crystals reversed, such that water is in contact with the Ag-(0 0 0 1) faces instead,
and again find |∆xtlφ| ≈ 46.2 V/nm, see Fig. S1 (b). These results clearly demonstrate
that removing the net dipole of the simulation cell is not the same as providing polarity
compensation.2 As in the Dz = 0 ensemble, the mobile slab is unstable in this mirrored slab
2
(a)
(b)
FIG. S1. Electrostatic potential profile φ(z) obtained from the mirrored slab geometry at Ez = 0.
The dashed silver and pink lines indicate the positions of the Ag-(0 0 0 1) and I-(0 0 0 1) surfaces
respectively. The region −4.33 < z < +4.33 nm is occupied by water, and the remaining volume
(|z| > 7.42 nm) is occupied by vacuum. As in the single slab geometry at Dz = 0, we find
|∆xtlφ| ≈ 46.2 V/nm. In (a), the AgI crystals are arranged such that the I-(0 0 0 1) faces are in
contact with water, whereas in (b), the Ag-(0 0 0 1) faces are in contact with water.
geometry [see Fig. 4 (a)]. Moreover, we also find that the structure of water is similar to
that at Dz = 0, as shown by the bond orientation statistics in Fig. S2 (a). In Fig. S2 (b)
we also show the charge density profile ρq(z) = (LxLy)
−1∑
i qiδ(z − zi) close to one of the
Ag-(0 0 0 1) interfaces. Results are also shown for a single slab geometry with different fields
and slab mobility, for comparison.
II. DERIVATION OF THE CNC CONDITIONS
Consider the system shown schematically in Fig. S3, which comprises a polar crystal
surrounded on either side by an aqueous electrolyte solution. The system is placed under
three dimensional periodic boundary conditions, however, for the ensuing discussion, we are
only concerned with the periodicity in the direction normal to the crystal’s surfaces. The
3
(a)
(b)
FIG. S2. (a) Bond orientation statistics in the mirrored slab geometry with water in contact with
Ag-(0 0 0 1) [see Fig. S1 (b)]. The results are consistent with the Dz = 0 ensemble [see Fig. 5 (c)].
A value of cos θ = −1 indicates a bond directed away from the surface, and cos θ = +1 indicates a
bond directed toward the surface. (b) Charge density profile ρq(z) from the Ag-(0 0 0 1) interface.
The z-coordinate has been shifted such that the Ag-(0 0 0 1) interface is located at z = 0, and
the liquid extends to z < 0. Results are shown for different simulation geometries and fields, as
indicated by the legend.
length of the periodic supercell in this direction—which we take to be the z direction—is Lz.
We consider the ‘insulator centered cell’ (ICS) where the crystal is centered at z = 0, and
the cell spans −Lz/2 < z < Lz/2. The crystal is made up of n+ 1 planes of ions, each with
a surface charge density ±σ0. (n is an odd integer.) As the crystal is a dielectric material, it
can sustain finite electric fields in its interior. Here, we assume the field alternates between
values E1 and E2 between the different ionic planes, as shown in Fig. S3. We have also
assumed that the dielectric constant xtl is constant throughout the material. On the other
hand, the surrounding electrolyte is a conductor, and the electric field in its interior strictly
vanishes. In other words, the electrolyte screens the crystal’s surface charge. This screening
induces planes of surface charge density ±σ near the surface of the crystal, and the regions
between these planes and the crystal are referred to as the Helmholtz layers. The electric
field in the Helmholtz layer is denoted EH, and its dielectric constant by H. As the electric
displacement field must be continuous across the boundary separating the Helmholtz layer
and the electrolyte, EH 6= 0 implies a finite polarization in the electrolyte. [This can be seen
most readily by considering 4piP = lim→∞ −1( − 1)D.] An important implication of this
4
FIG. S3. Schematic of the continuum model used in the derivation of the CNC conditions. As
in the main article (Fig. 1), the crystal comprises alternating layers with charge density ±σ0. E1
and E2 are the electric fields between layers in the crystal (assumed to have constant dielectric
constant xtl), with respective polarizations P1 and P2. EH is the electric field in the Helmholtz
layer (situated a distance `H from the crystal), which also has a polarization PH. The leftmost and
rightmost orange dashed lines indicate the cell boundaries. Due to the uniform polarization in the
electrolyte, there are formally charge densities ±σ at these boundaries.
is that there are formally surface charge densities of ±σ at the boundaries of the simulation
cell (see Fig. S3).
This configuration of charges produces dielectric polarization in each continuum region
that alters the net charge at the boundaries. It is these net charges that determine the
dipole moment. Explicitly they are,
σ1 = σ0 − PH − P1,
σ2 = σ0 − P1 − P2,
σH = σ − PH.
The dipole moment of the simulation cell is then,
Pz,cell =
1
L
∑
i
σizi. (S1)
5
Performing this sum leads to,
LPz,cell = −σ1
[(
n+ 1
2
)
R1 +
(
n− 1
2
)
R2
]
+σH
[
2`H +
(
n+ 1
2
)
R1 +
(
n− 1
2
)
R2
]
+σ2
(
n− 1
2
)
R2. (S2)
In the case that R1 = R2, Eq. S2 reduces to the result previously derived in Ref. 3.
One of the central themes of the main article was the notion of polarity compensation.
In macroscopically sized crystals, this means that the average field Extl = −∆xtlφ/w across
the crystal must vanish in order to avoid a divergent surface energy. Clearly, the potential
drop across the slab is the sum of electric fields inside the slab, appropriately weighted by
the width of each region,
wExtl =
(
n+ 1
2
)
R1E1 −
(
n− 1
2
)
R2E2, (S3)
where the width of the slab is given as
w =
(
n+ 1
2
)
R1 +
(
n− 1
2
)
R2. (S4)
By construction, at CNC the potential drop across the crystal is zero. This allows us to
write,
(n+ 1)R1E1 = (n− 1)R2E2. (S5)
We now express the surface charge densities by applying the Maxwell displacement equations
across the boundaries, resulting in the following set of simultaneous equations,
HEH = 4piσ (S6)
HEH + xtlE1 = 4piσ0 (S7)
xtlE2 + xtlE1 = 4piσ0. (S8)
We can solve Eq. S5 for E1 and then substitute it into Eq. S8 to find E2 in terms of R1, R2,
σ0 and xtl. Solving the simultaneous equations S6–S8 ultimately yields,
σCNC =
(n+ 1)R1
(n+ 1)R1 + (n− 1)R2σ0, (S9)
6
where we have added the subscript ‘CNC’ to acknowledge that this is the value of σ that
ensures ∆xtlφ = 0. We again note that for R1 = R2 we recover the result obtained in Ref. 3.
We will now show how σCNC is related to DCNC, the value of the electric displacement
field that enforces CNC conditions. This was derived in Ref. 4 for the case R1 = R2. Here
we generalize the result to R1 6= R2. Ultimately we will find that the difference between
the two scenarios manifests itself entirely in the value of σCNC. In other words, for the ICS
geometry we will find,
DCNC = −4piσCNC, (S10)
with the value of σCNC determined by Eq. S9. To demonstrate this, we note the Maxwell
field equations across the boundaries,
EH = 4piσH (S11)
EH + E1 = 4piσ1
E2 + E1 = 4piσ2.
We can use these to write Pz,cell given by Eq. S2 in terms of electric fields rather than surface
charge densities,
4piLPz,cell = −E1
(
n+ 1
2
)
R1 + 2EH`H + E2
(
n− 1
2
)
R2. (S12)
Like the potential drop across the crystal, the total potential drop across the cell ∆cellφ =
−LEz can itself be written as the sum of all electric fields, appropriately weighted by the
widths of the regions they occupy. This leads to,
Ez =
1
L
[
E1
(
n+ 1
2
)
R1 − 2`HEH − E2
(
n− 1
2
)
R2
]
, (S13)
allowing us to write,
4piPz,cell = −Ez. (S14)
At this point it is important to realize that Ez is the Maxwell field and includes all possible
contributions. Thus Ez comprises contributions from both the cell polarization and applied
fields. In the current context, this manifests itself in the form of the induced charges ±σ
at the cell boundaries. The full polarization must also include the charge induced at the
boundaries of the supercell,
4piPz = −4piσ − Ez. (S15)
From the Maxwell relation Dz = Ez + 4piPz, Eq. S10 follows immediately when enforcing
CNC conditions.
7
(a)
(b)
FIG. S4. Trial-and-error approach for locating ECNC. (a) φ(z) for Ez = 0 V/nm (dark blue) and
Ez = ECNC (cyan). (See Fig. 2.) At Ez = 0 V/nm there exists a potential drop ∆xtlφ across the
crystal, as indicated. (b) ∆xtlφ vs Ez exhibits an approximately linear relationship, as indicated
by the dashed line. The cyan and dark blue data points correspond to ECNC and Ez = 0 V/nm,
respectively.
III. ESTABLISHING CNC CONDITIONS UNDER CONSTANT Ez
In the main paper, we discussed the trial-and-error procedure for finding ECNC. This is
displayed pictorially in Fig. S4. To help orient the reader, Fig. S4 (a) is a replica of Fig. 2.
To find ECNC we simply perform simulations at different values of Ez and measure the value
of ∆xtlφ that results. The value of ECNC that causes ∆xtlφ to vanish (or in practice, to
be sufficiently small e.g. |∆xtlφ| < 0.2 V), is the value of ECNC. In general this procedure
is greatly facilitated by the approximately linear relation between ∆xtlφ and ECNC, which
allows us to extrapolate to ECNC with just a few data points. Results from a typical trial-
and-error procedure are shown in Fig. S4 (b).
8
IV. VALUES OF Ez AND Dz USED IN SIMULATIONS OF CNC
CONDITIONS
We have performed a number of simulations with different values of Ez and Dz in order
to approximate CNC conditions. The different values used are summarized in Table S1,
where we also give the corresponding values of ∆xtlφ. Due to relaxation of the Ag
+ and
I– ions, in the case of the mobile slab there is no reason to expect the values of ECNC and
DCNC to remain unchanged from those for the immobile crystal. For the AgI + pure H2O
system, we find that Ez = −0.56 V/nm is a good approximation for ECNC in both cases,
but that DCNC changes from Dz = −14.99 V/nm to Dz = −15.74 V/nm. For the AgI +
electrolyte system, on the other hand, using the same value of Ez = −0.31 V/nm for the
immobile and mobile slab leads to a potential difference of ∆xtlφ ≈ +0.39 V across the
crystal in the latter case. This was the value used in the simulations of ice formation at
ECNC with a mobile slab. The average value of the electric displacement field at 298 K was
found to be 〈Dz〉Ez=−0.31 ≈ −15.88 V/nm. Performing a simulation with a mobile slab at
Ez = −0.26 V/nm (which gives ∆xtlφ ≈ −0.04 V) gives 〈Dz〉Ez=−0.26 ≈ −15.51 V/nm, and
this was the value used in the simulations of ice formation at DCNC with a mobile slab.
Nevertheless, we observe similar results between the two ensembles, suggesting that the AgI
crystal—and ice formation in its presence—is sufficiently robust to such relatively small
potential drops across the crystal, at least for the simulation cell sizes investigated in the
main article. See e.g. Figs. S5-S10.
V. RESULTS NOT PRESENTED IN THE MAIN ARTICLE
In this section we give results not presented in the main article. Figs. S5–S8 show orienta-
tion statistics at Ag-(0 0 0 1) for different combinations of fields, temperatures, and whether
or not the slab is mobile. The reader is referred to the figure captions for details. Good
agreement with those results presented in the main manuscript is seen in all cases. Fig-
ure S9 shows the number of ice-like molecules, Nice, vs. time, t, in the absence of dissolved
ions. To determine Nice, we have used the CHILL+ algorithm.
5 As the crystal structure
of AgI resembles that of ice, we have included the Ag+ and I– ions in the analysis, and
we have also shifted the y coordinates of the Ag+ and I– ions by −0.2651 nm to enhance
9
TABLE S1. Summary of Ez and Dz values used for the CNC simulations presented in the main
article (reported in V/nm), along with corresponding ∆xtlφ (reported in V). For the case of a
mobile AgI slab in contact with NaCl electrolyte, these were performed at a value close to ECNC,
and have therefore been marked with an asterisk (see text).
System ECNC DCNC ∆xtlφ at ECNC ∆xtlφ at DCNC
immobile AgI, pure H2O −0.56 −14.99 −0.06 +0.15
immobile AgI, NaCl electrolyte −0.31 −14.99 +0.02 +0.15
mobile AgI, pure H2O −0.56 −15.74 −0.08 +0.10
mobile AgI, NaCl electrolyte −0.31∗ −15.51 +0.39∗ −0.10
the performance of the algorithm at Ag-(0 0 0 1). (This shifting is only performed during
the CHILL+ analysis.) Figure S10 shows the corresponding plots for the electrolyte system.
Both of these results are discussed in Sec. IIIB in the main article. In Fig. S11 we show
snapshots from simulations (DCNC and ECNC for pure H2O in contact with AgI) in which ice
formation is seen to occur at I-(0 0 0 1) as well as Ag-(0001). We stress that ice formation
is still seen to occur preferentially at Ag-(0001), even in the simulations shown in Fig. S11.
Due to this prefrence for ice formation at Ag-(0001), all the results for P (cos θ) shown in
the main article and in this supporting information have focused on the contact layer at this
interface. In Fig. S12 we present P (cos θ) for the contact layer at I-(0 0 0 1) from Dz = 0, and
under CNC conditions. We see analogous behavior to that reported at the Ag-(0 0 0 1), i.e.,
at Dz = 0 no water molecules direct O–H bonds toward the negatively charged interface,
whereas under CNC conditions, a significant proportion do. In Fig. S13, we present orienta-
tion statistics of water dipoles at Ag-(0 0 0 1) for Dz = 0 and at DCNC. This is an alternative
way of elucidating water’s structure at the interface, rather than looking at statistics of O–H
bonds.
All results presented so far have been obtained with the TIP4P/2005 water model,6 and the
Madrid model for the NaCl ions.7 While we may expect the thermodynamics and kinetics of
ice formation to depend sensitively on force field parameters, the CNC conditions have been
derived from a continuum model (see Sec. II). We therefore expect the CNC conditions to
be robust to the force field used for the water and dissolved ions. To demonstrate this, we
10
have performed simulations at 298 K using 864 SPC/E water8 (qO = −0.8476e) and three
Joung-Cheatham NaCl ion pairs9 (qNaCl = 1e). The resulting φ(z) for an immobile crystal
are shown in Fig. S14, both for at Dz = 0 and DCNC = −14.99 V/nm. At Dz = 0, we find
|∆xtlφ| = 46.2 V, while at DCNC we find |∆xtlφ| = 0.15 V. These are in excellent agreement
with the results obtained with the TIP4P/2005 and Madrid models. In Fig. S15 we show
the corresponding orientation statistics of O–H bonds at Ag-(0 0 0 1). Similar behavior is
observed compared to results obtained with the TIP4P/2005 and Madrid models.
In addition to being insensitive to the water and NaCl force field parameters, the CNC
conditions also appear to be independent of Lx and Ly. Indeed, doubling the extent of
the simulation cell in these dimensions while keeping Lz fixed (such that the system is four
times as large) results in the φ(z) profiles shown in Fig. S16. These results for the immobile
crystal agree excellently with the those found in the main article (see Fig. 3). Importantly,
the field in the solvent is the same for both simulation setups, Ez,solv ≈ −0.39 V/nm. This
leads to similar structures as reported in the main article, as seen in Fig. S17. With such
a large field in the solvent, we expect a similar proton ordered ice structure to form in this
system too. For the mobile crystal, however, we have found that ECNC changes slightly from
−0.56 V/nm to −0.55 V/nm, whereas for the system studied in the main article, we found
no discernible difference upon moving to a mobile crystal in the presence of pure water (see
Table S1).
VI. FURTHER SIMULATION DETAILS
As discussed in Sec. V of the main article, we used the reparametrized10 PRV potential11
to model AgI, as detailed in Ref. 12. To this end, we found it convenient to used a look-up
table to compute non-electrostatic interactions for the PRV potential. These interactions
take the form,
uPRV(r) =
H
rη
− C
r6
− P
r4
(S16)
where r is the interatomic separation, and H, η, C and P are parameters that take specific
values for the different types of interactions (i.e. Ag+–Ag+, Ag+–I– and I––I– ). Like all
other non-electrostatic interactions, these were truncated and shifted at r = 10 A˚. uPRV
was then tabulated with a spacing ∆r = 0.02 A˚. In Fig. S18 (a) we show the output from
LAMMPS of the interpolated potential, at a resolution δr ≈ 0.015 A˚. It is clear that the
11
(a) (b)
FIG. S5. Orientation statistics for pure H2O in contact with an immobile AgI crystal, at ECNC.
(a) Simulation performed at 298 K with water in its liquid state. (b) After ice formation at 242 K.
Good agreement with the simulations performed at DCNC is observed, see Figs. 5 (d) and 6 (d).
(a) (b) (c)
FIG. S6. Orientation statistics of water (NaCl electrolyte) in the contact layer and away from the
crystal, in the presence of an immobile AgI crystal. Simulations at 298 K at (a) DCNC and (b)
ECNC. In both cases, no orientation order is observed away from the surface. No water molecules
in the contact layer direct O–H bonds toward the surface. (c) After ice formation at 242 K at
ECNC. These results are consistent with those seen at DCNC, see Fig. 7 (b).
interpolation agrees well with the analytic form given by Eq. S16. In Fig. S18 (b) we show
radial distribution functions g(r) between the different species obtained for molten AgI at
923 K (250 AgI ion pairs, in a cubic simulation box with dimension L = 26.1069 A˚). These
appear to agree well with those reported by Bitria´n and Trulla`s.12
In Table S3 we report the parameters for the Lennard-Jones (LJ) interactions between
12
(a) (b)
FIG. S7. Orientation statistics after ice formation for pure H2O in the contact layer and away from
the crystal, in the presence of a mobile AgI crystal. Simulations have been performed both at (a)
ECNC and (b) DCNC. The results agree well with those presented in Fig. 6 (d), where an immobile
AgI crystal was used.
TABLE S2. Parameters used in the PRV potential, see Eq. S16. η is dimensionless. H, C and P
are given in eV A˚η, eV A˚4 and eV A˚6, respectively.
interaction pair η H P C
Ag+–Ag+ 11 0.16 0 0
Ag+–I– 9 1310.0 14.9 0
I––I– 7 5328.0 29.8 84.5
different species, where the LJ potential is,
uLJ(r) = 4LJ
[(
σLJ
r
)12
−
(
σLJ
r
)6]
. (S17)
Interactions between H2O, Na
+ and Cl– use the Madrid model.7 Interactions between H2O,
Ag+/I– use the parameters of Hale and Kiefer,13 as reported by Fraux and Doye.14 Remaining
interactions between Na+/Cl– and Ag+/I– were computed using Lorentz-Berthelot mixing
rules. The hydrogen atoms of the water molecules do not interact via the LJ potential, and
have been omitted from Table S3.
13
(a) (b)
FIG. S8. Orientation statistics of water (NaCl electrolyte) in contact with a mobile AgI crystal,
at ECNC. (a) Simulation performed at 298 K with water in its liquid state. Good agreement with
the simulations performed with an immobile slab is observed, see Fig. S6. (b) After ice formation
at 242 K. Again, good agreement is seen with the immobile slab simulations. Moreover, the results
are also consistent with those performed at DCNC, see Fig. 7 (b).
TABLE S3. LJ interaction parameters (LJ, σLJ) reported to four decimal places. LJ is reported
in kcal/mol, and σLJ is reported in A˚ngstrom. ‘−’ indicates no LJ interaction is specified between
this pair of species.
O Na+ Cl– Ag+ I–
O (0.1852, 3.1589) (0.1896, 2.5134) (0.0148, 4.2687) (0.5467, 3.1710) (0.6220, 3.3420)
Na+ (0.3519, 2.2174) (0.3439, 2.9051) (0.7536, 2.7002) (0.8574, 2.8712)
Cl– (0.0184, 4.8491) (0.1723, 4.0161) (0.1960, 4.1871)
Ag+ − −
I– −
REFERENCES
1I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 (1999).
2J. Goniakowski, F. Finocchi, and C. Noguera, Rep. Prog. Phys. 71, 016501 (2008).
3T. Sayer, C. Zhang, and M. Sprik, J. Chem. Phys. 147, 104702 (2017).
14
(a)
(b)
(c)
FIG. S9. Nice vs t for pure H2O in contact with AgI. Blue, green and orange lines indicate
the results obtained with an immobile AgI crystal, while purple, brown and pink indicate results
obtained with a mobile AgI crystal. (a) Simulations performed at Dz = 0 (immobile slab only).
(b) Simulations performed at DCNC. (c) Simulations performed at ECNC. The black dashed line
indicates the number of AgI ions, which have also been included in the CHILL+ analysis.
15
(a)
(b)
FIG. S10. Nice vs t for NaCl electrolyte in contact with AgI. Color scheme as in Fig. S9. (a)
Simulations performed at DCNC. (b) Simulations performed at ECNC. In contrast to the pure
water system, ice formation in the presence of ions appears to be preceded by an induction time
followed by rapid growth. The black dashed line indicates the number of AgI ions, which have also
been included in the CHILL+ analysis.
4T. Sayer, M. Sprik, and C. Zhang, J. Chem. Phys. 150, 041716 (2019).
5A. H. Nguyen and V. Molinero, J. Phys. Chem. B 119, 9369 (2014).
6J. L. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).
7A. Benavides, M. Portillo, V. Chamorro, J. Espinosa, J. Abascal, and C. Vega, J. Chem.
Phys. 147, 104501 (2017).
8H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
16
(a) (b) (c)
FIG. S11. Examples of ice formation at I-(0 0 0 1) for pure H2O in contact with AgI. We stress
that ice formation occurs preferentially at Ag-(0 0 0 1), even in the simulations shown in this figure.
(a) and (b) show snapshots from simulations performed at DCNC and ECNC, respectively, with an
immobile crystal. (c) Snapshot from simulation performed at ECNC with a mobile crystal.
9I. S. Joung and T. E. Cheatham III, J. Phys. Chem. B 112, 9020 (2008).
10F. Shimojo and M. Kobayashi, J. Phys. Soc. Jpn. 60, 3725 (1991).
11M. Parrinello, A. Rahman, and P. Vashishta, Phys. Rev. Lett. 50, 1073 (1983).
12V. Bitria´n and J. Trulla`s, J. Phys. Chem. B 112, 1718 (2008).
13B. N. Hale and J. Kiefer, J. Chem. Phys. 73, 923 (1980).
14G. Fraux and J. P. Doye, J. Chem. Phys. 141, 216101 (2014).
17
−1.0 −0.5 0.0 0.5 1.0
cos θ
0.0
0.5
1.0
1.5
2.0
P
(c
os
θ)
Dz = 0, pure water
DCNC, pure water
ECNC, electrolyte
FIG. S12. Orientation statistics of water at I-(0 0 0 1), at Dz = 0 and under CNC conditions
(see legend). The results are analogous to those at Ag-(0 0 0 1): At Dz = 0, essentially no water
molecules direct O–H bonds toward the surface, whereas under CNC conditions, many water
molecules do direct O–H bonds toward the negatively charged surface. Note that in contrast to
our analysis at Ag-(0 0 0 1), cos θ = −1 indicates a bond directed toward the surface, and cos θ = +1
indicates a bond directed away from the surface.
18
(a) (b)
(c) (d)
FIG. S13. Orientation statistics of water dipoles at Ag-(0 0 0 1). θµ is the angle between a water
dipole and the z-axis of the simulation cell. (a) and (b) show results for liquid water at 298 K
for Dz = 0 and at DCNC, respectively. (c) and (d) show results after ice formation at 242 K for
Dz = 0 and at DCNC, respectively. A value of cos θ = −1 indicates a dipole directed away from
the surface, and cos θ = +1 indicates a dipole directed toward the surface. All results used an
immobile crystal.
19
−4 −2 0 2 4
z [nm]
−20
−10
0
10
20
φ
(z
)
[V
]
DCNC
Dz = 0
FIG. S14. φ(z) using SPC/E water and Joung-Cheatham NaCl (three ions pairs) at 298 K. The
striking similarity to the results obtained with TIP4P/2005 and the Madrid model demonstrate the
robustness of the CNC conditions to the underlying atomic model of the electrolyte. Results have
been obtained with an immobile crystal.
20
(a)
(b)
FIG. S15. Orientation statistics of water (NaCl electrolyte) in contact with an immobile AgI
crystal, at 298 K. Results have been obtained with the SPC/E and Joung-Cheatham force fields.
21
−4 −2 0 2 4
z [nm]
−20
−10
0
10
20
φ
(z
)
[V
]
ECNC
DCNC
Dz = 0
FIG. S16. Comparing φ(z) from different ensembles for AgI (n = 17) in contact with pure water
at 298 K, with the system doubled in the x and y dimensions. At Dz = 0 we find ∆xtlφ ≈ −46.2 V.
At ECNC we find ∆xtlφ ≈ −0.04 V, and |Ez,solv| ≈ 0.39 V/nm. The result obtained at DCNC =
−14.99 V/nm (dotted line) agrees well with the ECNC result.
(a) (b)
FIG. S17. Bond orientation statistics at 298 K for a system at whose extent in x and y has been
doubled. (a) and (b) show P (cos θ) at Dz = 0 and DCNC, respectively, both for water molecules
in the contact layer (blue circles) and in bulk solvent (orange squares). cos θ = +1 and cos θ = −1
indicate O–H bonds directed immediately toward and away from Ag-(0 0 0 1), respectively. All
results have been obtained with an immobile crystal.
22
(a) (b)
FIG. S18. (a) Non-electrostatic portion of uPRV vs interatomic separation. The symbols have been
output from LAMMPS at a resolution of approximately 0.0015 nm, while the solid lines show the
analytic form of the PRV potential (truncated and shifted at r = 1 nm). (b) Radial distribution
functions from simulations of molten AgI at 923 K. These appear to agree well with those reported
in Ref 12.
23