Kiran Mathew

V S Chaitanya Kolluru

Srinidhi Mula

Stephan N Steinmann

Richard G Hennig

Implicit self-consistent electrolyte model in plane-wave density-functional theory

The ab-initio computational treatment of electrochemical systems requires an appropriate treatment of the solid/liquid interfaces. A fully quantum mechanical treatment of the interface is computationally demanding due to the large number of degrees of freedom involved. In this work, we describe a computationally efficient model where the electrode part of the interface is described at the density-functional theory (DFT) level, and the electrolyte part is represented through an implicit solvation model based on the Poisson-Boltzmann equation. We describe the implementation of the linearized Poisson-Boltzmann equation into the Vienna Ab-initio Simulation Package (VASP), a widely used DFT code, followed by validation and benchmarking of the method. To demonstrate the utility of the implicit electrolyte model, we apply it to study the surface energy of Cu crystal facets in an aqueous electrolyte as a function of applied electric potential. We show that the applied potential enables the control of the shape of nanocrystals from an octahedral to a truncated octahedral morphology with increasing potential.

I. INTRODUCTION

Interfaces between dissimilar systems are ubiquitous in nature and frequently encountered and employed in scientific and engineering applications. Of particular interest are solid/liquid interfaces that form the crux of many technologically essential systems such as electrochemical interfaces, 1-3 corrosion, lubrication and nanoparticle synthesis. [4][5][6] A comprehensive understanding of the behavior and properties of such systems requires a detailed microscopic description of the solid/liquid interfaces.

Given the heterogeneous nature and the complexities of solid/liquid interfaces, 7 experimental studies often need to be supplemented by theoretical undertakings. A complete theoretical investigation of such interface based on an explicit description of the electrolyte and solute involves solving systems of nonlinear equations with a large number of degrees of freedom, which quickly becomes prohibitively expensive even when state of the art computational resources are employed. 8 This calls for the development of efficient and accurate implicit computational frameworks for the study of solid/liquid interfaces. [8][9][10][11][12][13][14] We have previously developed a self-consistent implicit solvation model 15,[START_REF] Mathew | VASPsol -Solvation model for the plane wave DFT code VASP[END_REF] that describes the dielectric screening of a solute embedded in an implicit solvent, where the solute is described by density-functional theory (DFT). We implemented this solvation model into the widely used DFT code VASP. [START_REF] Kresse | [END_REF] sol, is freely available under the open-source Apache License, Version 2.0 and has been used by us and others to study the effect of the presence of solvent on a number of materials and processes, such as the catalysis of ketone hydrogenation, 18 ionization of sodium su-peroxide in sodium batteries, 19 the chemistry at electrochemical interfaces, 20,21 hydrogen adsorption on platinum surfaces, 22 surface stability and phase diagram of Mica in solvent, 23 etc. 24,25 In this work, we describe the extension of our solvation model VASPsol 15,[START_REF] Mathew | VASPsol -Solvation model for the plane wave DFT code VASP[END_REF] to include the effects of mobile ions in the electrolyte through the solution of the linearized Poisson-Boltzmann equation and then apply the method to electrochemical systems to demonstrate the power of this approach. Figure 1 illustrates how the solid/electrolyte system is divided into spatial regions that are described by the explicit DFT and the implicit solvent. This solvation model allows for the description of charged systems and the study of the electrochemical interfacial systems under applied external voltage. Section II outlines the formalism and presents the derivation of the energy expression. Sec. III describes the details of the implementation and Sec. IV validates the method against previous computational studies of electrochemical systems. Finally, in Sec. V we apply the approach to study how an applied electric potential can control the equilibrium shape of Cu nanocrystals.

II. THEORETICAL FRAMEWORK

Our solvation model couples an implicit description of the electrolyte to an explicit quantum-mechanical description of the solute. The implicit solvation model incorporates the dielectric screening due to the permittivity of the solvent and the electrostatic shielding due to the mobile ions in the electrolyte. The solute is described explicitly using DFT.

To obtain the energy of the combined solute/electrolyte system, we spatially divide the materials system into three regions: (i) the solute that we describe explicitly using DFT, (ii) the electrolyte that we describe using the linearized Poisson-Boltzmann equation and (iii) the interface region, that electrostatically couples the DFT and Poisson-Boltzmann equation. Fig. 1 illustrates the regions for a Pt(111) surface embedded in a 1 M aqueous electrolyte. The interface between the electrolyte and the solute system is formed self-consistently through the electronic density, and the interaction between the electrolyte and solute is described by the electrostatic potential as explained in the following sections.

A. Total energy

We define the interface region between the solute and the electrolyte using the electronic charge density of the solute, n( r). The electrolyte occupies the volume of the simulation cell where the solute's electronic charge density is essentially zero. In the interface region where the electronic charge density increases rapidly from zero towards the bulk values inside of the solute, the properties of the implicit solvation model, e.g., the permittivity, , and the Debye screening length, λ D , change from their bulk values of the electrolyte to the vacuum values in the explicit solute region. The scaling of the electrolyte properties is given by a shape function 10,15,[26][START_REF] Sundararaman | JDFTx[END_REF][START_REF] Gunceler | [END_REF]

ζ(n( r)) = 1 2 erfc log(n/n c ) σ √ 2 . (1)
Following our previous work 15 and the work by Arias et al., 10,[26][START_REF] Sundararaman | JDFTx[END_REF][START_REF] Gunceler | [END_REF] the total free energy, A, of the system consisting of the solute and electrolyte is expressed as

A[n( r), φ( r)] = A TXC [n( r)] + φ( r)ρ s ( r) d 3 r - ( r) |∇φ| 2 8π d 3 r + 1 2 φ( r)ρ ion ( r) d 3 r + A cav + A ion , (2)
where A TXC is the kinetic and exchange-correlation contribution from DFT, φ is the net electrostatic potential of the system, and ρ s and ρ ion are the total charge density of the solute and the ion charge density of the electrolyte, respectively. The total solute charge density, ρ s , is the sum of the solute electronic and nuclear charge densities, n( r) and N ( r), respectively,

ρ s ( r) = n( r) + N ( r). (3)
The ion charge density of the electrolyte, ρ ion , is given by

ρ ion ( r) = i qz i c i ( r), (4)
where c i ( r) is the concentration of ionic species i, z i , denotes the charge state, and q is the elementary charge. The concentration of ionic species, c i , is given by a Boltzmann factor of the electrostatic energy 29 and modulated in the interface region by the shape function, ζ(n( r)),

c i ( r) = ζ[n( r)]c 0 i exp -z i qφ( r) k B T , (5)
where c 0 i is the bulk concentration of ionic species i, k B the Boltzmann constant, and T the temperature.

The relative permittivity of the electrolyte, ( r), is assumed to be a local functional of the electronic charge density of the solute and modulated by the shape function, 10,15,[26][START_REF] Sundararaman | JDFTx[END_REF][START_REF] Gunceler | [END_REF] (n( r)

) = 1 + ( b -1)ζ(n( r)), ( 6
)
where b is the bulk relative permittivity of the solvent. The cavitation energy, A cav , describes the energy required to form the solute cavity inside the electrolyte and is given by

A cav = τ |∇ζ(n( r))|d 3 r, (7)
where the effective surface tension parameter τ describes the cavitation, dispersion, and the repulsion interaction between the solute and the solvent that are not captured by the electrostatic terms. 13 In comparison to our previous solvation model, 15 the main changes to the free energy expression are the inclusion of the electrostatic interaction of the ion charge density in the electrolyte with the system's electrostatic potential, φ, and the non-electrostatic contribution to the free energy from the mobile ions in the electrolyte, A ion . In our work we assume this non-electrostatic contribution from the ions to consist of just the entropy term,

A ion = k B T S ion , (8)
where S ion is the entropy of mixing of the ions in the electrolyte, which for small electrostatic potentials can be approximated to first order as 29

S ion = i c i ln c i c 0 i d 3 r ≈ - i c i z i qφ k B T d 3 r. (9)

B. Minimization

To obtain the stationary point of the free energy, A, given by Eq. ( 2), we set the first order variations of the free energy with respect to the system potential, φ, and the electronic charge density, n(r), to zero. 10,15,[26][START_REF] Sundararaman | JDFTx[END_REF][START_REF] Gunceler | [END_REF] Taking the variation of A[n( r), φ( r)] with respect to the electronic charge density, n( r), yields the typical Kohn-Sham Hamiltonian 30 with the following additional term in the local part of the potential

V solv = δ (n) δn |∇φ| 2 8π + φ δρ ion δn + τ δ|∇ζ| δn + k B T δS ion δn . (10)
Taking the variation of A[n( r), φ( r)] with respect to φ( r) yields the generalized Poisson-Boltzmann equation

∇ • ∇φ = -ρ s -ρ ion , (11)
where (n(r)) is the relative permittivity of the solvent as a local functional of the electronic charge density.

We further simplify the system by considering the case of electrolytes with only two types of ions present, whose charges are equal and opposite, i.e. c 0 1 = c 0 2 = c 0 and z 1 = -z 2 = z. Then the ion charge density of the electrolyte becomes 29

ρ ion = ζ[n( r)]qzc 0 exp -zqφ k B T -exp zqφ k B T = -2ζ[n( r)]qzc 0 sinh zqφ kT (12)
and the Poisson-Boltzmann equation becomes

∇ • ∇φ = -ρ s + 2ζ[n( r)]qzc 0 sinh zqφ k B T . (13)
For small arguments x = zqφ k B T 1, sinh(x) → x and we obtain the linearized Poisson-Boltzmann equation

∇ • ∇φ -κ 2 φ = -ρ s (14)
with

κ 2 = ζ[n( r)] 2c 0 z 2 q 2 k B T = ζ[n( r)] 1 λ 2 D , (15)
where λ D is the Debye length that characterizes the dimension of the electrochemical double layer.

The implicit electrolyte model is given by the additional local potential, V solv , of Eq. ( 10) and the linearized Poisson-Boltzmann equation given by Eqs. ( 14) and ( 15). This model has four key approximations: (i) The ionic entropy term of Eq. ( 9) treats the ionic solution as an ideal system and assumes that any interactions between the ions are small compared to k B T . (ii) The cations and anions of the electrolyte have equal and opposite charge, simplifying the Poisson-Boltzmann equation. (iii) The electrostatic potential in the electrolyte region is small, such that zqΦ k B T , which leads to the linear approximation of the Poisson-Boltzmann equation (14). To be quantitative, arguments to the sinh function in Eq. ( 13) of less than 0.25 lead to errors below 1%. Hence, in the region where the shape function is unity (and thus ∇• ≈ 0), charge excesses of up to 0.5 c 0 are still faithfully reproduced. (iv) The ions are treated as point charges. The finite size of the ions in solution would limit their maximum concentration in the electrolyte, and the model could lead to unphysically high concentrations near the electrodes. These four approximations could be overcome and will be considered in future extensions of the VASPsol model.

III. IMPLEMENTATION

We implement the implicit electrolyte model described above into the widely used DFT software Vienna Abintio Software Package (VASP). [START_REF] Kresse | [END_REF] VASP is a parallel plane-wave DFT code that supports both ultra-soft pseudopotentials 31,32 as well as the projector-augmented wave (PAW) 33 formulation of pseudopotentials. Our software module, VASPsol, is freely distributed as an open-source package and hosted on GitHub at https://github.com/ henniggroup/VASPsol. [START_REF] Mathew | VASPsol -Solvation model for the plane wave DFT code VASP[END_REF] The main modifications in the code are the evaluation of the additional contributions to the total energy and the local potential, given by Eqs. ( 2) and (10), respectively. Corrections to the local potential require the solution of the linearized Poisson-Boltzmann equation given by Eq. ( 13) in each self-consistent iteration. Our implementation solves the equation in reciprocal space and makes efficient use of fast Fourier transformations (FFT). 15 This enables our implementation to be compatible with the Message Passing Interface (MPI) and to take advantage of the memory layout of VASP. We use a pre-conditioned conjugate gradient algorithm to solve the linearized Poisson-Boltzmann equation with the preconditioner G 2 + κ 2 -1 , where G is the magnitude of the reciprocal lattice vector and κ 2 is the inverse of the square of the Debye length, λ D . The solution of the linearized Poisson-Boltzmann equation (13) provides a natural reference electrostatic potential by setting the potential to zero in the bulk of the electrolyte. However, plane-wave DFT codes, such as VASP, implicitly set the average electrostatic potential in the simulation cell to zero, not the potential in the electrolyte region. For simplicity, we do not modify the VASP reference potential but provide the shift in reference potential that needs to be added to the Kohn Sham eigenvalues and the Fermi level. Furthermore, the shift of the electrostatic potential to align the potential in the electrolyte region to zero, or to any other value, modifies the energy of the system. This energy change is given by ∆E ref = n electrode ∆U ref , where n electrode is the net charge of the electrode slab and ∆U ref is the change in reference potential, i.e., the shift to align the potential in the electrolyte region to zero.

Electrode Potential U (V) vs. PZC Charge (e) Grand canonical energy F-F 0 (meV) 0 -1 -2 -3 -4 -5 -6 -7 -0.3 -0.2 -0
To validate the change in reference potential, Fig. 2 shows the grand canonical electronic energy, F (U ), as a function of the applied potential, U , of the Pt (111) electrode slab. The grand-canonical electronic energy, F , is the Legendre transformation of the free energy, A, of the system, F (U ) = A(n)-n electrode U , and is always lowered when charge is added or removed from the neutral slab. The grand canonical energy, F , exhibits the expected quadratic behavior 34 and the maximum at the neutral slab when the energy change ∆E ref due the change in reference potentials is included.

IV. VALIDATION

We further validate the model by comparing with existing experimental and computational data. First, we compute the potential of zero charge (PZC) of various metallic slabs and compare them against the experimental values. Second, we compute the surface charge density of a Pt(111) slab as function of the applied external potential and compare the resulting value of the double layer capacitance with previous computations and experiments.

The DFT calculations for these two benchmarks and the following application are performed with VASP and the VASPsol module using the PAW formalism describing the electron-ion interactions and the PBE approximation for the exchange-correlation functional. [START_REF] Kresse | [END_REF]33,[START_REF] Perdew | [END_REF] The Brillouin-zone integration employs an automatic mesh with 50 k-points per inverse Ångstrom with only one kpoint in the direction perpendicular to the slab.

We consider an electrolyte that consists of an aqueous solution of monovalent anions and cations of 1 M concentration. At room temperature this electrolyte To properly resolve the interfacial region between the solute and the implicit electrolyte region and to obtain accurate values for the cavitation energy requires a high plane-wave basis set cutoff energy of 1000 eV. However, in our validation calculations and applications, we observe that the contribution of the cavitation energy to the solvation energies are negligibly small. We, therefore, set the effective surface tension parameter τ = 0 for all following calculations. This also removes the requirement for such a high cutoff energy, and we find that a cutoff energy of 600 eV results in converged surface energies and PZC.

To change the applied potential, we adjust the electron count and then determine the corresponding potential from the shift in electrostatic potential as reflected in the shift in Fermi level. This takes advantage of the fact that the electrostatic potential goes to zero in the electrolyte region for the solution of the Poisson-Boltzmann equation, providing a reference for the electrostatic potential.

Figure 3 compares the PZC for Cu, Ag, and Au surface facets calculated with the implicit electrolyte model as implemented in VASPsol with experimental data. The electrode PZC is defined as the electrostatic potential of a neutral metal electrode and is given by the Fermi energy with respect to a reference potential. A natural choice of reference for the implicit electrolyte model is the electrostatic potential in the bulk of the electrolyte, which is zero in any solution of the Poisson-Boltzmann equation. In experiments, the reference is frequently chosen as the standard-hydrogen electrode (SHE). The dashed line in Fig. 3 presents the fit of U pred = U exp +∆U pred SHE to obtain the shift, ∆U pred SHE , between the experimental and the computed PZC, i.e. the potential of the SHE for our electrolyte model. The resulting ∆U pred SHE = 4.6 V compares well with the previously reported computed value of 4.7 V for a similar solvation model. [START_REF] Gunceler | [END_REF] Figure 4 shows the calculated surface charge density, σ, of a Pt(111) surface as a function of applied electrostatic potential, U . We find that the PZC for the Pt(111) electrode is 0.85 V, in good agreement with previous computational studies. 26,[START_REF] Gunceler | [END_REF] However, this value is at the high end of the experimentally reported results. [START_REF] Trasatti | Modern Aspects of Electrochemistry[END_REF]37 This discrepancy is likely due to adsorption on Pt(111). The work by Sekong and Groß shows that Pt(111) is covered at potentials below 0.5 V by hydrogen, followed by the socalled double layer region between 0.5 and 0.75 V with no adsorbates present and OH adsorption at more positive potentials. 25 The presence of adsorbates in the experiments alters the PZC and precludes a direct comparison with our calculations.

The linear slope of σ(U ) is the double-layer capacitance. Fitting a quadratic polynomial to the data and evaluating the slope at the pzc yields a double-layer capacitance for the Pt(111) surface at 1 M concentration of 14 µF/cm 2 . This agrees perfectly with the previously reported computed result [START_REF] Gunceler | [END_REF] and is close to the experimental value of 20 µF/cm 2 . 37 Our implementation of the electrolyte model neglects the finite volume of the ions in the electrolyte. Hence, at larger applied potentials, the model may predict unphysically large ion concentrations to compensate the surface charge. To understand how this approximation limits the applied potential, we perform calculations on a Cu(111) slab in an electrolyte with 2 M concentration for a range Position (Å) 111) surface, at a potential of about -1.9 V vs. the potential of the SHE. The plane-wave cutoff of 400 eV leads to numerical noise, which disappears at 600 eV.

of applied potentials from -3 to +3 V.

Figure 5 shows the extrema of the net ion concentration in the double-layer region as a function of the applied potential (written by VASPsol into the file RHOION). The negative and positive concentrations correspond to excess cations and anions near the surface, respectively. For a large range of potentials from -2 to +3 V, the net ion concentration depends nearly linear on the applied potential. At negative potentials below -2 V, we observe that electrons leak from the slab into the electrolyte region, rendering the results of the electrolyte model unphysical. Hence, care should be taken when utilizing the model at large negative potentials.

The inset in Fig. 5 illustrates the average net ion concentration along the direction perpendicular to the Cu(111) surface for an applied potential of -1 V. The ion concentration decays exponentially away from the solid-liquid interface and approaches zero in the bulk electrolyte. We selected an electrolyte region of 30 Å to illustrate the need to converge the results with respect to the size of the electrolyte domain. Here, the size is not quite sufficient to reach a negligible ion concentration. While this does not affect the considered extrema of the ion concentration, care must be taken to ensure convergence of the results with the size of the electrolyte region. As general guidance, an electrolyte region > 10λ should suffice. 38 At a large applied potential of +3 V, we find a net ion concentration of 0.2 M, which is 10% of the electrolyte concentration, far from the 50% excess that would introduce a 1% error as discussed following Eq. 14. Here, the concentrations of anions and cations are 2.1 and 1.9 M, respectively, and the error introduced by the linear approximation to the sinh function is less than 0.1%. Thus, even for a considerable potential of +3 V, the linear electrolyte model does not result in unphysically large ion concentrations.

Finally, we validate the accuracy of the implemented analytic force evaluation. As an example, we have cho- sen the adsorption of a Na atom on the Au(111) surface. The surface charge was set to 0.4 e -for the symmetric p(2 × 2) surface, corresponding to an electrochemical potential of about -1.9 V vs. the SHE, while the neutral Na@Au(111) surface corresponds to a potential of -2.6 V. Fig. 6 shows that the agreement between the analytical gradient and the minimum for the energy is excellent, provided a sufficiently high cutoff energy of 600 eV is used.

V. CRYSTAL SHAPE CONTROL BY APPLIED POTENTIAL

The implicit electrolyte model, VASPsol, enables the study of electrode surfaces in realistic environments and under various conditions, such as of an electrode immersed in an electrolyte with an external potential applied. To illustrate the utility of the solvation model, we determine how an applied potential changes the equilibrium shape of a Cu crystal in a 1 M electrolyte. The crystal shape is controlled by the surface energies. We show that the surface energies are sensitive to the applied potential and that each facet follows different trends, leading to opportunities for the practical control of the shape of nanocrystals.

Using the MPInterfaces framework, 40 we construct slabs of minimum 10 Å thickness for the(111), (100), (110), (210), ( 221), (311) and (331) facets of Cu. The different facets for Cu are chosen based on the crystal facets included in the MaterialsProject database. 39 We calculate the surface energy of each slab in vacuum and electrolyte by relaxing the top and bottom layers and keeping the middle layers of atoms fixed to their bulk positions. We employ a vacuum spacing of 30 Å to ensure that the electrostatic interactions between periodic images of the slabs are negligibly small. From the surface energies, we determine the resulting shape using the Wulff construc- tion as implemented in the pymatgen framework. 39 Figure 7 compares the crystal shape of Cu in vacuum obtained in our calculations with MaterialsProject. The MaterialsProject calculations used a cutoff energy of 400 eV and a vacuum spacing of 10 Å which both are lower than the cutoff energy of 600 eV and the vacuum spacing of 30 Å used in this work. We find that our predicted surface energies agree with MaterialsProject within 2% for all the facets. However, we note that even such small deviations in surface energies matter for the details of the crystal shape, where the high index (310) and (221) facets do not show up in our calculation of the equilibrium crystal shape of Cu. These small but observable deviations demonstrate that the crystal shape is sensitive to small variations in the surface energies.

For the calculation of the Cu surface energies in the 1 M electrolyte, we vary the number of electrons in the slabs and determine the change in the corresponding applied electrostatic potential as described in Sec. III. We calculate the surface energies for each of the above facets in vacuum over a range of applied potentials of approximately ±1 eV around the corresponding PZC. We then perform a spline interpolation to obtain the surface energies as a function of the applied potential. This transformation enables the comparison of the surface energies for a given applied potential and is made possible by the absolute reference potential provided by the solvation model. From the results we obtain the equilibrium crystal shape as a function of applied potentials by the Wulff construction.

For copper, the (111) facets exhibit the lowest surface energy over the whole potential range. Figure 8 shows how the ratios of the surface energies of the copper (hkl) facets relative to the (111) facet vary as a function of

Copper crystals at applied potential -1 -0.4 0.4 U (V) (111) (100)
FIG. 9. Variation of the equilibrium crystal shape of Cu with applied electrode potential. The equilibrium shape transitions from an octahedron at negative potentials to a truncated octahedron upon increasing applied potentials.

applied potential. As expected, the surface energies of all facets systematically increase with potential, however, their ratios relative to the (111) facet decreases. Due to this, the crystal shape changes as a function of applied potential. Figure 9 illustrates the change in crystal shape with increasing applied potential. At negative potentials, the (111) facets dominate and lead to an octahedral shape. With increasing potential, the ratio of the (100) to (111) surface energies decreases sufficiently to lead to a change in shape to a truncated octahedron. The transition occurs around a potential of -0.56 V for Cu. With further increase of potential, the area of the (100) facets increases. Even though the ratio of the surface energies of other facets are also lowered, the change is not sufficient for other facet to show up in the equilibrium crystal shape.

The prediction for the changes in the equilibrium shape of Cu crystal as a function of applied electric potential demonstrate the utility of the VASPsol solvation model for computational electrochemistry simulations using density-functional theory and provides the opportunity to design nanocrystals of specific shape through the application of an electrostatic potential during growth at an applied potential.

VI. CONCLUSIONS

Solid/liquid interfaces between electrodes and electrolytes in electrochemical cells present a complex system that benefits from insights provided by computational studies to supplement and explain experimental observations. We developed and implemented a selfconsistent computational framework that provides an efficient implicit description of electrode/electrolyte interfaces within density-functional theory through the solution of the linearized Poisson-Boltzmann equation. The electrolyte model is implemented in the widely used DFT code VASP, and the implementation is made available as a free and open-source module, VASPsol, hosted on GitHub at https://github.com/ henniggroup/VASPsol. We showed that this model enables ab-initio electrochemistry studies at the DFT level. We validated the model by comparing the calculated electrode potentials of zero charge for various metal electrodes with experimental values and the calculated double-layer capacitance of Pt(111) in a 1 M solution to previous calculations and experiments. We tested the validity of the linear approximation at large applied potentials and the accuracy of the analytic expressions for the forces. To illustrate the usefulness of the model and the implementation, we apply the model to determine the equilibrium shape of Cu crystal as a function of applied electrostatic potential. We predict a change in shape from an octahedron to a truncated octahedron with increasing potential. These calculations illustrate the utility of the implicit VASPsol model for computational electrochemistry simulations.

2 FIG. 1 .

 21 FIG.1. Spatial decomposition of the solid/electrolyte system into the solute, interface, and electrolyte regions for a fcc Pt(111) surface embedded in an electrolyte with relative permittivity r = 78.4 and monovalent cations and anions with concentration 1 M. The solute is described by densityfunctional theory, the electrolyte by an implicit solvation model. The interface region is formed self-consistently as a functional of the electronic density of the solute. All properties are shown as a percentage of their maximum absolute value. The inset shows the relatively invisible peak in the Hartree potential at the interface.

FIG. 3 .

 3 FIG.3. Comparison between the computed and the experimental potential of zero charge (PZC) with respect to the standard hydrogen electrode (SHE). The dashed line is a fit of U pred = U exp + ∆U pred SHE to determine the theoretical potential of the SHE, ∆U pred SHE . The experimental values are taken from Ref.35.

10 FIG. 4 .

 104 FIG. 4. Surface charge density of the Pt(111) surface as a function of applied electrostatic potential computed with the implicit electrolyte module, VASPsol. The slope of the curve determines the capacitance of the dielectric double layer.

FIG. 5 .

 5 FIG.5. The extrema of the ion concentrations over the potential range of -3 to +3 V relative to the PZC in an electrolyte with ion concentration of 2 M. The negative and posconcentrations correspond to excess of cations and anions, respectively, in the double layer. The solid blue line is a linear fit between -1 and +2 V, indicating the large linear potential range. The inset shows the xy averaged ion concentration profile along the z-direction at a potential of -1 V; the central region corresponds to the solute slab with zero ion concentration.

FIG. 6 .

 6 FIG.6. Analytical gradient and relative energy as a function of the distance of a Na ion from the Au(111) surface, at a potential of about -1.9 V vs. the potential of the SHE. The plane-wave cutoff of 400 eV leads to numerical noise, which disappears at 600 eV.

FIG. 7 .

 7 FIG.7. Comparison of the predicted shape of a Cu crystal in vacuum (a) for this work obtained by the Wulff construction using the surface energies of the (100), (110), (111), (210), (221), (311), and (331) facets compared with (b) the Materi-alsProject database.39

FIG. 8 .

 8 FIG.8. The surface energies, γ hkl , of the various facets of Cu with respect to surface energies, γ111, of the (111) facet as a function of applied potential. The dashed horizontal line indicates the ratio of the (100) and (111) surface energies, γ100/γ111 equals √ 3, which corresponds to a change in shape from an octahedron to a truncated octahedron.

  The grand canonical electronic energy, F (U ) of a charged Pt(111) slab relative to the neutral slab as a function of the electrode potential U . The grand canonical energy displays the maximum correctly at the neutral slab when the reference electrostatic potential is chosen such that the potential is zero in the bulk of the electrolyte and the necessary correction due to the reference potential ∆E ref = n electrode ∆U ref is included in the energy.

	-0.04	-0.02	0.0		0.02	0.04
		.1	0.0	0.1	0.2	0.3
	FIG. 2.

ACKNOWLEDGMENTS This work was supported by the National Science Foundation under awards Nos. DMR-1542776, ACI-1440547, and OAC-1740251. This research used computational resources of the Texas Advanced Computing Center under Contract Number TG-DMR050028N and of the University of Florida Research Computing Center. Part of this research was performed while the authors were visiting the Institute for Pure and Applied Mathematics (IPAM), which is supported by the National Science Foundation (NSF).


Ruben Staub

Stephan N Steinmann
email: stephan.steinmann@ens-lyon.fr

Parameter-Free Coordination Numbers for Solutions and Interfaces

Keywords: numbers: 68, 08, De

Coordination numbers are among the central quantities to describe the local environment of atoms and are thus used in various applications such as structure analysis, fingerprints and parameters. Yet, there is no consensus regarding a practical algorithm, and many proposed methods are designed for specific systems. In this work, we propose a scale-free and parameter-free algorithm for nearest neighbor identification. This algorithm extends the powerful Solid-Angle based Nearest-Neighbor (SANN) framework to explicitly include local anisotropy. As such, our Anisotropically corrected Solid-Angle based Nearest-Neighbor (ASANN) algorithm provides with a fast, robust and adaptive method for computing coordination numbers. The ASANN algorithm is applied to flat and corrugated metallic surfaces to demonstrate that the expected coordination numbers are retrieved without the need for any system-specific adjustments. The same applies to the description of the coordination numbers of metal atoms in AuCu nano-particles and we show that ASANN based coordination numbers are well adapted for automatically counting neighbors and the establishment of cluster expansions. Analysis of classical molecular dynamics simulations of an electrified graphite electrode reveals a strong link between the coordination number of Cs + ions and their position within the double layer, a relation that is absent for Na + , which keeps its first solvation shell even close to the electrode.

I. INTRODUCTION

Detailed analysis of chemical systems allows to gain insight into their working mechanisms, to rationalize their properties and sometimes even to predict the impact of chemical modifications such as functionalization of molecules or alloying of surfaces. Roughly two big categories can be distinguished: energetic and structural characterizations. While energy decomposition analysis provides detailed energetic insight in molecular chemistry, solid interfaces and in liquid phase, [1][2][3][4] it is the coordination environment that commonly defines the building block for local structural interpretations. The main applications in chemistry involve hydration numbers of ions 5,6 , which are of particular interest in the context of biochemistry, 7 the link between structure and activity in heterogeneous catalysis 8,9 and the characterization of alloys. 10,11 Intuitively, the coordination number of a particle is defined as the number of nearest neighboring particles. However, this does not provide with a rigorous definition: what is neighboring after all? Chemists tend to link neighboring to the idea of bonding: two items are connected if there is a strong interaction between them. But again, "strong" is not defined and most of the time the strength of the interaction is not easily accessible. Therefore, it is far more common and convenient to think in terms of spatial proximity, echoing the initial terminology of "neighbors".

In chemistry and physics, many different algorithms are used to define the connectivity between atoms, using more or less additional parameters (e.g., nature of the atoms or reference bond lengths). We propose to classify the algorithms based on their input, using only two criteria: (i) Usage of parameters (e.g., particles nature, reference bond lengths). Any algorithm that depends on more than a set of point coordinates falls into this category and is denoted with a P 1 , while parameter-free algorithms are denoted by P 0 . (ii) Local adaptivity. By this we mean that the connectivity between two particles depends on the presence of other surrounding particles and denote it by A 1 , while the absence of local adaptivity corresponds to A 0 . Since these two criteria are non-exclusive, Figure 1 provides relevant examples of algorithms for the four possible combinations.

The simplest coordination number algorithm is based on tabulated "typical" bond lengths: only if the atoms are closer, they are connected. Hence, it falls into the class P 1 A 0 . This simple and rapid algorithm is most popular in visualization softwares, 12 where FIG. 1: Overview of available coordination number algorithms. We report only popular algorithms of relevance for the context of the ASANN algorithm development, with the references provided in the text. bond length cut-offs are typically constructed from predefined covalent radii 13 of considered elements.

Instead of using predetermined typical/average bond lengths that might not be of relevance for the considered sample, coordination cut-offs can be estimated from the sample itself. This lead to the development of multiple parameter-free methods 7,10,11,14-19 using a cut-off based on the radial distribution function g(r). However, such algorithms require enough data to perform meaningful radial distribution analysis. Finally, these algorithms have an inherently mean-field description of coordination since they use averaged cut-offs.

Hence, local fluctuations are not necessarily well handled as these algorithms do not take the remaining neighboring particles into account, i.e., it is a P 0 A 0 algorithm.

Constructing locally adaptive coordination numbers requires an analysis of the surroundings of each possible connection, which necessarily increases the computational cost. The best known algorithm to extract context-dependent coordination numbers uses a Voronoi tesselation, 18,20 where space is decomposed into regions (polyhedra for 3D structures), based on which particle is the closest. Two particles are connected if their corresponding regions share a face. This algorithm is parameter-free and thus falls into the P 0 A 1 category, which means it can be applied without any additional knowledge on the set of particles considered.

However, it has been shown that Voronoi-based coordination numbers have a tendency to be overestimated, [18][19][20][21] which can be alleviated by re-weighting. 11 Furthermore, 3D Voronoi tesselations are computationally expensive [START_REF] Ledoux | 4th International Symposium on Voronoi Diagrams in Science and Engineering[END_REF] . In the same P 0 A 1 category, the Solid-Angle based Nearest-Neighbor algorithm (SANN) algorithm 21 was developed as an efficient alternative to Voronoi tesselations. SANN attributes a solid angle to each possible neighbor and iteratively increases the cut-off radius to reach a sum of 4π for the solid angles, corresponding to a fully surrounded coordination environment. This elegant method is faster than the Voronoi tesselation and, furthermore, reduces the overestimation of the coordination numbers, without the introduction of any parameters 21 . However, SANN suffers from coordination number overestimation at interfaces [START_REF] Higham | [END_REF] .

Recently, the Relative Angular Distance (RAD) algorithm [START_REF] Higham | [END_REF] was developed to amend the SANN overestimation at interfaces. The reduction of the coordination number is achieved by accounting for blocking: the coordination of two particles can be obstructed by particles between the considered pair. Compared to other available blocking-based methods 24,25 RAD does not require predefined cut-off distances, but still requires additional parameters, so that the algorithm falls in the P 1 A 1 class. An earlier P 0 A 1 algorithm 26 accounted for blocking by considering only bonds that are not intersecting a triangle formed by shortest bonds. However, this idea was not further explored as its computational cost is expected to be large, and we found the method not suitable to describe surfaces.

Herein we propose to augment the SANN algorithm with a correction term that takes the local anisotropy into account, while preserving all properties of SANN, i.e., low computational cost, no parameters, being locally adaptive (P 0 A 1 class) and being scale free, i.e., an isotropic scaling of all coordinates does not change the coordination number. This scheme is called Anisotropically corrected Solid-Angle based Nearest-Neighbors (ASANN).

To set the stage, we introduce the SANN algorithm in section II A. Then, in section II B, we identify the SANN isotropic description as source of coordination number overestimation and propose a natural solution to tackle local anisotropy through ASANN. After the computational details in section III, section IV contains tests and applications of ASANN for characterizing metal surfaces, alkali metal ion coordination at an electrified interface and for the construction of model Hamiltonians for alloy nanoparticles.

II. ALGORITHMS

A. The SANN algorithm

The reader is strongly encouraged to read the original SANN article 21 for details. Nevertheless, the main points are outlined below to make the ASANN algorithm understandable. SANN defines a coordination number CN SANN (i) of particle i to be the number of the m nearest neighbors found within the individual coordination sphere of radius R (m) i (also called shell radius). SANN elegantly provides a self-consistent definition of m and the associated coordination sphere radius R (m) i as detailed below.

Let r i,j be the distance between particle i and its j-th nearest neighbor (denoted simply as j in the following). By definition of R (m) i and m, one can already write the inequality:

r i,m ≤ R (m) i < r i,m+1 (1)
Within SANN, the m nearest neighbors are defined such that the sum of their solid angles (with respect to i) is equal to 4π. This requirement is related to the intuitive idea that the nearest neighbors of i are the smallest set of particles fully covering its field of view. Indeed, each nearest neighbor j of particle i is associated with a solid angle Ω i,j , which is defined by

Ω i,j = 2π(1 -cos(θ i,j
)), where the angle θ i,j is depicted on Figure 2. surrounding particles are considered neighbors if and only if they are within the coordination radius).

The above definitions are summarized in the following equation:

4π = m j=1 Ω i,j = m j=1 2π(1 -cos(θ i,j )) = m j=1 2π 1 - r i,j R (m) i (2)
One can rewrite the two fundamental Eqs 1 and 2 into a much simpler form:

R (m) i = m j=1 r i,j m -2 < r i,m+1 (3)
The inequality of Eq. 3 allows to easily compute the CN SANN (i) = m by iteratively increasing R (m) i from m = 3 until the inequality is respected (see Figure 3). The convergence of the algorithm was mathematically proven 21 .

SANN is scale free, i.e., a uniform scaling of all coordinates will not affect the coordination number. This property can have counter-intuitive consequences when the analyzed phase features significant density fluctuations. This is particularly problematic when applied to the gas-liquid interface, where thermal fluctuations might create nearly isolated particles.

The scale-free algorithm will assign a significantly larger coordination sphere to such an isolated particle, which potentially leads to a higher coordination number compared to particles in dense regions. However, such problematic cases can be spotted by monitoring the coordination sphere radius.

Furthermore, in situations where the environment of j is more dense than the one of i, it is possible that j is considered a neighbor of i, even though i is not considered a neighbor of j. Further along this line, we note that the solid-angle based approach behind SANN is inherently linked to the assumption of locally identical coordination radii. This hypothesis is the only meaningful approach within a parameter-free approach and represents a first order expansion for the coordination sphere. As a consequence, directional bonding (e.g., in organic molecules) which leads to locally more compact geometries is can lead to an overestimation of the coordination number compared to conventional views.

B. Anisotropically corrected Solid-Angle based Nearest-Neighbors, ASANN algorithm

The SANN algorithm implicitly uses a locally isotropic description of the coordination shell, i.e., the angular distribution of neighboring particles is not taken into account. At surfaces or tips there are, however, simply no neighbors in some directions, leading to empty sections in the coordination sphere. Therefore, the SANN radius must expand to reach more neighbors in the non-empty directions in order to compensate and reach a total sum of 4π. Hence, one observes an over-estimation of coordination numbers when neighboring particles are not isotropically distributed. [START_REF] Higham | [END_REF] This issue was already reported in the literature, and led to the development of the relative angular distance (RAD) algorithm. However, RAD requires the particles sizes as input parameters, and is fundamentally different from SANN. Here, we propose an Anisotropically corrected Solid-Angle based Nearest neighbors (ASANN) algorithm, which is based on the SANN framework, but includes a parameter-free correction term to correctly describe intrinsically anisotropic coordination shells.

The basic idea of the ASANN algorithm is that some directions are more relevant for coordination than others. Therefore, instead of requiring the neighbors to fill up the whole field of view, ASANN only requires to fill up the relevant neighboring region of space (see Figure 4).

An estimator of the first order local anisotropy is obtained based on the barycentre. context of the Voronoi tessalationfor its canonicity. 11 Mathematically, this can be written:

G (m) i = m j=1 Ω i,j r i,j m j=1 Ω i,j (4)
This SANN-based barycentre is in agreement with the properties expected for retrieving the first order local anisotropy. Indeed, it identifies a global privileged direction in the distribution of coordinating neighbors and, more importantly, it provides with a quantitative estimation of the first order local anisotropy magnitude α (m) i defined by:

α (m) i = || G (m) i || R (m) i (5)
Instead of considering a coordination sphere, the ASANN algorithm uses a coordination spherical cap C which focuses on the locally relevant region of space for coordination. C i is the spherical cap centered on i having the same anisotropic properties as the neighboring particles through our first order local anisotropy estimator. This is achieved by setting the radius

R C i = R (m) i and requiring that its barycentre G C i is confounded with the SANN-based barycentre G (m) i : C i is centred on i R C i = R (m) i G C i = G (m) i (6)
The solid angle Ω C i associated with this spherical cap (with respect to i) can be written as:

Ω C i = 4π(1 -γ C i ) (7)
where γ C i ∈ [0, 1[ corresponds to the angular correction coefficient, which only depends on the first order local anisotropy magnitude α (m) i

(see section S1 in the SI for a derivation):

γ C i = α (m) i + α (m) i 2 + 3α (m) i 3 (8)
Therefore, instead of filling up the whole field of view (

m j=1
Ω i,j = 4π), ASANN neighbors only cover the spherical cap. This leads to the fundamental equation of ASANN:

m j=1 Ω i,j = m j=1 2π 1 - r i,j R (m ) i = Ω C i = 4π(1 -γ C i ) (9)
where the ASANN coordination number m and the coordination cap radius R (m ) i are defined by extending the SANN coordination radius definition to include the angular correction

coefficient γ C i : R (m ) i = m j=1 r i,j m -2 × (1 -γ C i ) < r i,m +1 (10)
Note that for isotropic distributions (γ C i = 0), Eq. 10 goes back to Eq. 3, which means that SANN and ASANN coincide for this case.

Using the inequality of Eq. 10, the ASANN coordination number m can be easily com-

puted by iteratively computing R (m ) i , starting with m = 2(1 -γ C i ) + 1,

and increasing

m until the inequality is respected.

In principle, this angular correction could be computed self-consistently, i.e., the angular correction term is recomputed at each iteration along with the ASANN coordination radius.

This self-consistency breaks the analytic, direct link between the fundamental equation (Eq. 9) and the coordination radius definition (Eq. 10). Therefore, we only explore the "per- which also provides information regarding the computational cost of ASANN compared to SANN: Essentially, there is a factor of two between them, since the iterations have to be performed twice. Knowing that ASANN coordination number is lower or equal to the one from SANN (see SI for a proof), there are, in general, less ASANN than SANN iterations.

Based on general properties of ASANN, partially inherited from SANN, we show the expected applicability domain of the two algorithms in Figure 6. The first major limitation is related to low coordination numbers (e.g., four) in combination with isotropic distributions (e.g., a tetrahedron). In this case, the coordination radius R

i of Eq. 3 is twice the distance between atom i and its nearest neighbors j. If there are atoms closer than 2r ij , the inequality is not satisfied and the coordination shell is expanded to include next nearest neighbors. This is the case for most organic molecules (see Figs. S2-S5 and the associated discussion), but also in ionic compounds where small cations are located in tetrahedral sites as shown in the SI (Tables S4 to S15), where we have analyzed a variety of salts, oxides and perovskites.

The particularity of these systems is most clearly seen for the case of a graphene sheet (Fig. S5): the barycenter of neighboring atoms is confounded with the central atom, so there is no first order anisotropy and ASANN is equivalent to SANN. There are three nearest neighbor atoms at a distance r C-C and the next nearest neighbors are already at √ 3r C-C , again due to the directional bonding, but which is well below the R

(3) i = 3r C-C (Eq. 3 for
the perfect symmetry of this system). If there would not be an anisotropic preference due to bond formation characteristic for an sp 2 carbon atom, there would also be "space" for at least two more particles one above and one below the plane within such a coordination sphere. But this is not the case and, therefore, for geometrical reasons, ASANN significantly overestimates the conventional coordination number of 3, by adding the 6 next nearest neighbors, reaching a CN ASANN =9 for carbon in graphene. The same geometrical reasons and deviations from homogeneous distributions apply to the ionic compounds referred to above. Since ASANN is built on a first order anisotropy, the application to nanopores or slits is also at the border of its applicability: if the SANN coordination sphere spans two surfaces, the first order anisotropy vanishes and ASANN does not improve upon SANN.

Higher order anisotropy corrections might correct SANN in this situation as well as for the isotropic, low-coordination environments mentioned above. Second, density fluctuations are problematic for the scale free algorithms. Therefore, cross-checks with the coordination radii are required for biphasic systems and the liquid-gas interface. As we will show by numerical results below, ASANN extends, however, the applicability of SANN to (inhomogeneous) solutions, the solid-liquid interface and rough surfaces including nano-particles.

The time complexity of (A)SANN is in the worst-case O(nk) for pre-treated input, where n is the number of particles, and k is the maximum coordination number. In fact, the most time consuming step is due to the input pre-processing, required for both SANN and ASANN algorithms. This pre-processing comprises the computation of pair-wise distances (and vectors for ASANN) in O(n 2 ), and sorting those distances in O(n 2 log(n)). It is, therefore, this step that is also responsible for the space complexity of (A)SANN in O(n 2 ). In practice, this pre-processing quickly becomes the computational bottleneck for large systems. As a consequence, particular care was taken for optimizing these steps in our Python 3 implementation, heavily relying on the Numpy library [START_REF] Oliphant | Guide to NumPy[END_REF][START_REF] Walt | [END_REF] . We thus combine low-level language performance for large systems, with convenient Python 3 syntax and modularity.

As an example, the computation of the Cs coordination numbers in a 742 atom system (vide infra), took around 0.14 seconds and the computation of all coordination numbers in a 2942

atom system took about 2 seconds (93% is spent on input pre-processing) on a personal laptop (Core i5-6300HQ CPU).

III. COMPUTATIONAL DETAILS

The presented algorithms have been implemented using Python3 [START_REF] Van Rossum | Python 3 Reference Manual[END_REF] and are available in the supporting information and also on the webpage of the authors. The implementation relies on the Numpy library 27,28 whenever possible and ASE 12 for reading common coordinate formats. A particular effort was made to provide with extensively commented code, including documentation and illustrative examples.

Figures have been generated with a custom ASE 12 fork to visualize the coordination numbers.

All DFT computations are performed using VASP [START_REF] Kresse | [END_REF]31 . The PBE exchange-correlation functional 32 was applied and a 400 eV energy cut-off was chosen for the plane-wave basis set. The PAW pseudo-potentials were used for describing the electron-ion interactions. 33,34 Some of the optimized metallic surfaces and metal clusters have been previously published by our group. [35][36][37] The graphite electrode was set up as follows: A symmetric six layer p(2×2) graphite slab has been built based on the optimized bulk geometry using the PBE-dDsC 38 level of theory. A 9 × 9 × 1 Monkhorst-Pack K-point grid was used to integrate the Brillouin zone. To assess the charge distribution at negative potentials, the number of electrons was varied and the linearized Poisson-Boltzmann equation was utilized to neutralize the unit cell as implemented in the VASPsol package. 39,40 Similarly to our previous study on electroreduction in aprotic medium, 41 the isodensity value defining the cavity has been lowered to 5 10 -4 e/ Å3 . In the case of graphite, this value ensures an absence of implicit solvent between the dispersion bound graphene layers.

The potential of zero charge of graphite is, according to our setup, -0.2 V vs SHE, which compares reasonably well with experiments (0.1 to 0.3 V depending on the nature of the carbon electrode [42][43][44] ) considering the idealized morphology of our electrode and the expected systematic error due to the simplified electrolyte model as discussed in the context of the potential dependent adsorption of pyridine on a gold single crystal surface 45 . A negative charge of -0.25 electrons in the unit cell (-0.125 electrons per surface) led to a potential drop of 2.4 V compared to the potential of zero charge and thus corresponds to a significant negative polarization. Based on the DFT computations, the unit cell was replicated four times in each in-plane direction, resulting in 128 carbon atoms per layer and a total charge of 4 e-.

Following a setup very similar to ref 46, the Lennard-Jones (LJ) parameters for graphite are taken from UFF. 47 The initial solvent distribution is obtained from the predefined TIP3P 48 box with about 35 Å of water surrounding the system, resulting in a box of 716 water molecules. To simulate a 1 M electrolyte, 15 Na + (or Cs + ) and 11 Cl -ions, described by their standard AMBER force field parameters, are added using tleap of the AmberTools.

MM computations are performed with NAMD 2.9 49 for 5 ns (time step of 2 fs), applying a Langevin thermostat (300 K) with a damping coefficient of 1 ps -1 . The Langevin barostat

(1 bar) is applied with a piston period of 100 fs and a decay of 150 fs.

IV. RESULTS AND DISCUSSION

A. Metallic surfaces

As a first test, we apply ASANN to various metallic surfaces. The selection of surfaces is intended to span a wide range of typical cases to assess the applicability domain of SANN and ASANN.

For face centered cubic (fcc) metallic surfaces, the metal-metal contact distance in the bulk is a reliable measure for identifying nearest neighbors. Reassuringly, closed packed, flat surfaces are well described by SANN and ASANN, giving identical results that are in perfect agreement with cutoff based methods (see SI). In other words, flat surfaces are not anisotropic enough to induce an expansion of the SANN coordination sphere to the next nearest neighbors.

An increased local anisotropy can be obtained for nano-rods. The large rod, constructed from the corresponding bulk without any geometry optimization (Figure 7), is still well and identically described by SANN and ASANN.

However, stepped surfaces (see Figures 8 and9) are more challenging for SANN, while ASANN consistently retrieves the cutoff based coordination number. In particular, SANN is very sensitive to small local surface induced contractions, while ASANN is more robust, avoiding coordination number overestimations seen in Figure 9 for SANN.

To simulate low-coordinated atoms, 9 add-atoms have also been analyzed. Based on visual inspection (and confirmed by cutoff based methods), the square of addatoms on a (100) surface is characterized by two nearest neighbors from the square (the third is not "touching", i.e., it is a next-nearest neighbor) and four from the (100) plane below, yielding a coordination number of six. Similarly, an addatom of the triangle on a (111) surface possesses two nearest neighbors from the triangle and three from the supporting (111) plane, resulting in a coordination number of five. Recalling that the flat (100) surface of fcc closed packed metals has a coordination number of 8 for its surface atoms, the CN SANN =8 for surface addatoms reveal a clear overestimation (see Figures 10 and11), qualitatively failing to identify the low-coordination environment. This contrasts with CN ASANN that is in full agreement with visual inspection and cuttoff based methods.

To challenge the algorithms even further, we apply them to an optimized, amorphous gold nanoparticle (see Figure 12). As can be seen, the particle features atoms that clearly form tips and, therefore, we expect low CNs for them. However, just like for the addatoms on the planar surfaces, SANN fails to correctly identify these low-coordinated atoms. For these highly aniostropic distributions, the SANN coordination radius has to expand significantly to As a rule of thumb, we empirically found that SANN rarely finds coordination numbers below 7 in these systems, illustrating the surface-related overestimation reported earlier. [START_REF] Higham | [END_REF] ASANN, however, performs very well in all tested systems, correctly identifying low-coordination environments. As shown in the supporting information, this robust, parameter-free performance is also transmitted to generalized coordination numbers 9 of various surfaces, dispensing with the need to define a characteristic distance or manual counting.

B. Coordination numbers at charged interfaces

The interaction of cations with negatively charged electrodes is of growing interest in the electro-reduction of CO and CO 2 . [50][51][52] In a different context, electrolyte solutions at charged silica interfaces are widely studied to understand the behavior of confined electrolytes. 53 Finally, the advent of carbon based supercapacitors 54,55 has spurred simulations of electrolyte solutions at charged carbon interfaces. 56,57 Such simulations are also of particular interest to understand the intercalation mechanism of alkali metal ions into graphitic surfaces. 58 Here, we consider a simple model system: aqueous salt solutions in contact with a perfect graphite electrode. On the one hand, the experimental capacitance of a graphitic electrode is not strongly dependent on the alkali cation in the electrolyte, 59 which might suggest that the double layer structure is similar for all alkali metal cations. On the other hand, Cs + ions are, in contrast to Na + ions, known to be soft and large, implying a low charge density at their surface. As a consequence, Cs + ions are expected to partially loose their solvation shell to approach the negatively charged electrode surface more closely, which increases both the VdW interactions and the Coulomb attraction.

We probe the difference between Cs + and Na + based on classical molecular dynamics simulations of a graphite carbon electrode at a potential of about -2.5 V (see Computational Details) in ∼1 M aqueous CsCl and NaCl solutions. To analyze the double layer, we investigate the existence of a link between coordination number of ions in water (determined by ASANN, SANN and a g(r) based cut-off algorithm) and their adsorption state. presence near the anode. The strong first cation excess peak is followed by a nearly charge compensated zone before a second, although lower, cation excess peak occurs which then decays to zero towards the bulk solution. Besides, the anion (Cl -) distribution is barely structured and independent on the nature of the cation. The details of the charge-excess profiles are quite different: Cs + excess density has the first peak at about 3.5 Å, while the corresponding one for Na + is at 4.8 Å. This difference indicates that Cs + can indeed approach the electrode more closely. The second peak is, accordingly, shifted from 6 Å to 7.5

Å when going from caesium to sodium. while the second layer is indistinguishable from the bulk solution in terms of the coordination numbers. This can be rationalized based on Van der Waals (VdW) radii: a fully solvated Cs ion cannot approach the graphitic electrode closer than roughly 7 Å (r C + r Cs + + 3 Å, where 3 Å is a rough estimate for the size of a water molecule).

For Na + ions, however, the coordination numbers are the same in the two excess regions as in the bulk, despite the qualitative similarity between the two cation excess profiles.

Even the nearest Na + layer (∼ 4.8 Å) is not close enough to the surface to induce a partial desolvation (see Table II). In other words, in contrast to Cs + , the solvation shell of Na + is strong and stable across the double layer and the three methods agree on an approximate average coordination number of six, indicating minor fluctuations and a rather isotropic distribution, in agreement with spectroscopic data. 6 To quantitatively determine if ASANN provides with "good" coordination numbers, we assess their ability to discriminate between "close" and "bulk" Cs ions. Specifically, we want to quantify if CN ASANN is statistically lower (i.e., below CN thresh ) at the interface than in the bulk. For the sake of comparison, we apply the same analysis also to SANN and to the cutoff based algorithm.

To have a quantifiable measure whether or not CN ASANN is lower at the interface than in the bulk, we explore the relation between closeness and loss of coordination via a Receiver Operating Characteristic (ROC) curve analysis. The ROC curve is one of the most rigorous and well-established tools to evaluate the discriminating power for a binary classifier (here Here, the TPR is the ratio of CN i <CN thresh for the i Cs ions "close" to the electrode to the number of Cs ions "close" to the electrode overall, while the False Positive Rate is the ratio of CN i <CN thresh for Cs ions not "close" to the electrode to the number of Cs ions not "close" to the electrode overall. The threshold value can be varied between the minimum and maximum observed with a given algorithm. Here, we define "close" as the Cs -surface distance ≤ 4.7 Å.

Such an analysis is reported on Figure 14. Let us first look at the ASANN results to explain how to read a ROC curve. Since the CNs are reasonable descriptors, for low CN thresh (e.g., 7), only "close" Cs ions are identified (FPR=0%), but not all of them (TPR≈20%), as shown in Table S16a. Similarly, at high CN (e.g., 10), the TPR reaches ≈100%, i.e., all "close" Cs ions have CNs below this value, but about 50% (FPR) of the ions not close to the surface have CNs below this value. Hence, increasing the threshold further does not change the TPR, but the FPR reaches 100%.

A major ROC curve analysis is performed through the Area Under the Curve (AUC), as it provides an aggregate measure of classification performance 63 . In this particular example, the AUC is associated with the probability that a random interfacial Cs ion has a lower coordination number than a random bulk Cs ion.

FIG. 14: ROC curve for classifying Cs ions for being "close" to the carbon anode based on the coordination number being below a given threshold. When moving from TPR=20% to TPR=55% with ASANN (blue curve), CN thresh is increased from 7 to 8 (see Table S16a in the SI for the raw data). This exploits the geometrical constraint responsible for a partial desolvation when close enough to the surface. The distance criterium for being close enough to the surface is chosen to be 4.7 Å, since it corresponds to the empirical frontier between the partially desolvated adsorption zone and the fully solvated adsorption zone for Cs ions near the graphite anode (see Fig. 13).

With an AUC of 0.90, we find a strong correlation between CN ASANN and the proximity of Cs + and the surface. This finding is coherent with ASANN providing reasonable coordination numbers in this system. The same analysis performed with SANN only yields an AUC of 0.81. The cutoff-based algorithm yields an AUC of 0.86, outperforming SANN, but not reaching the performance of ASANN. This is rather remarkable, given that the cutoff algorithm is tuned to the system by analyzing the radial distribution function, while ASANN does not know anything about the specificities of the system. On a physical level, the major difference between the cutoff algorithm and ASANN is that the ASANN coordination sphere adapts itself to the local environment, while the cutoff does not. This means that the cutoff based CNs are more sensitive to thermal fluctuations. Indeed, we found cutoff-based coordination numbers to be slightly noisier (≈ +8%) than their ASANN/SANN counterparts within each region of the studied system (see Table I), rationalizing the poorer AUC performance.

From this comparison it can be concluded that fixed-cutoff based algorithms are somewhat challenged when dealing with liquids, while SANN is not suitable for interfaces. In contrast, ASANN combines the best features of these algorithms and performs well. Therefore, ASANN is most relevant for describing solid/liquid interfaces.

C. Application to Au NP energy fitting

Exploring the energy landscape of alloy nano-particles from first principles is computationally very expensive. Therefore, cheap energy evaluations based on simplified Hamiltonians are often applied. [START_REF] Chen | [END_REF][65][66] . Furthermore, if the model Hamiltonian captures the essential physics, it can also be used to gain insight into the driving force of the alloy formation.

Inspired from cluster expansions, we here test the use of a simple model Hamiltonian H m for bimetallic materials. H m depends only on the number of atoms of each element and the number of bonds between each possible elements pairs:

H m = χ Γ χ N χ + χ,χ ≥χ Γ χ,χ N χ,χ (11)
In other words, it is a two-body model Hamiltonian. For alloys of elements with quite different atomic radii, the alloy nano-particles are not necessarily very regular and selecting characteristic distances is not straight forward. Therefore, an automatic, parameter free algorithm to identify bonds is particularly beneficial for these systems. In the following, we apply this simplified model to bimetallic Au-Cu materials, where H m becomes

H m = Γ Au N Au + Γ Cu N Cu + Γ Au,Au N Au,Au + Γ Au,Cu N Au,Cu + Γ Cu,Cu N Cu,Cu = Γ • N (12)
with Γ the set of parameters to fit. N is constructed by enumeration of atoms and bonds of each type. The automatic bond enumeration was performed using ASANN and SANN for comparison purposes. In order to have consistent bond description, a bond between atoms A and B is only counted if and only if B is found as nearest neighbor of A and reciprocally. S1 andS2 in the SI), but the overall model performance is better with ASANN-based bonds (MAD near 1.5 eV for ASANN, versus 1.8 eV for SANN). This is especially true when using only the nanoparticles as test set. The predicted energy for the nanoparticles with respect to the DFT reference energies yields a coefficient of determination R 2 = 0.72 and R 2 = 0.59 with ASANN and SANN based 2-body patterns, respectively (see Figure 15). Similarly, the pure bulks are better described by ASANN (MAD=0.37 eV), than SANN (MAD=0.60 eV). Hence, the better performance of SANN for exclusively fitting the nanoparticles set can be traced back to error cancellation between coordination number overestimation and missing higher-order many-body terms at the surface. In summary, we conclude that ASANN provides a balanced and physical sound description between bulk and surface-rich nanoparticles. Therefore, ASANN can also be applied to alloys of elements with quite different atomic radii. Similarly, ASANN could be used to provide a simple estimation of the surface of nanoparticles, by considering the union (or surface sum) of the planar edge of each coordination cap. A more rigorous approach using ASANN to define the surface of nanoparticles would be through the definition of a locally adaptative coordination radius, defining the local coarseness of the surface.

V. CONCLUSION

In this work we have extended the solid-angle based nearest-neighbor (SANN) algorithm to account for the local anisotropy, defining the anisotropically corrected solid-angle based nearest-neighbor (ASANN) method. ASANN allows to efficiently compute parameter-free coordination numbers and significantly improves over SANN at interfaces, where the coordination environment is highly anisotropic, such as nano-particles and partially desolvated ions at solid/liquid interfaces, where SANN leads to overestimations of the coordination number. In all situations tested, the ASANN coordination numbers are either identical or physically more sound than the ones provided by SANN and on par with system-specific, cutoff based coordination numbers. Therefore, ASANN provides a robust method to compute coordination numbers for describing alloy nano-particles, extended surfaces and (electrified) solid/liquid interfaces.

FIG. 2 :

 2 FIG. 2: Schematic representation of the solid angle definition Ω i,j = 2π(1 -cos(θ i,j )) associated with neighbor j at distance r i,j of central particle i, and using a coordination radius of R (m) i .

FIG. 3 :

 3 FIG.3: Illustration of the SANN iterative algorithm in 2D. Blue dots represents neighbors of the considered particle (green dot), while other particles are displayed in purple. For each increasing number of neighbors m (starting from 3), a coordination sphere radius is computed such that the SANN equality is respected (sum of neighboring solid angle requirement). The algorithm continues until the coordination radius is well-defined (i.e.

FIG. 4 :

 4 FIG. 4: Illustration of the SANN coordination number overestimation, and founding idea of the ASANN algorithm: tackle local anisotropy by focusing on relevant regions for coordination. This leads to the consideration of a coordination spherical cap instead of a sphere. The chosen spherical cap reproduces the local anisotropy, through our first order local anisotropy estimator (using a SANN-based barycentre with solid-angle weighting).

  FIG. 5: Schematic representation of the ASANN workflow: Based on the SANN coordination sphere, the barycenter and the corresponding angular correction is computed, leading to the ASANN coordination spherical cap.

FIG. 6 :

 6 FIG. 6: Venn diagram depicting the applicability domains of both SANN and ASANN. By design, ASANN extends the applicability domain of SANN to systems presenting local anisotropy.

FIG. 7 :

 7 FIG. 7: Transversal cut of a periodic nano-rod of an fcc metal, with superimposed coordination numbers. SANN and ASANN coordination numbers are identical.

FIG. 9 :

 9 FIG. 9: Side view of a stepped gold surface (644), with superimposed SANN (a) and ASANN (b) coordination numbers.

FIG. 10 :

 10 FIG. 10: Representation of a p(3×3) Cu(100) surface with a square of 4 Cu addatoms. Coordination numbers from SANN (a) and ASANN (b) are superimposed.

  FIG. 11: Representation of a p(3×3) Cu(111) surface with a triangle of 3 Cu addatoms. Coordination numbers from SANN (a) and ASANN (b) are superimposed.

FIG. 12 :

 12 FIG. 12: Low energy, amorphous Au 38 nanoparticle. Coordination numbers from SANN (a) and ASANN (b) are superimposed. SANN overestimation is particularly visible at the tips.

Figure 13

 13 Figure13displays the structuring of the electrolyte at the electrode surface. As commonly observed60 , for both electrolyte solutions the distribution curve of the excess cation density presents two well-defined peaks, corresponding to two distinct layers of preferred

FIG. 15 :

 15 FIG. 15: Parity plot for the total energy predicted by a 2-body model Hamiltonian fitted using the SANN vs. ASANN coordination algorithm. The whole set of 27 nanoparticles + 6 Bulk structures was used as training set, and the test set was composed of the 27 Au-Cu bimetallic nanoparticles only for graphical considerations.

Table I and

 I II report the ASANN and SANN coordination numbers of Cs + and Na + ,

		Average coordination numbers of Cs + ions
	Method	Desolvated trapping layer	Solvated trapping layer	Bulk-like layer
		h < 4.7	Å	4.7 < h < 7	Å	h > 7	Å
	ASANN	7.4 (1.1)	9.6 (1.2)		9.4 (1.2)
	SANN	9.1 (1.1)	10.6 (1.2)		10.5 (1.2)
	fixed cutoff	7.1 (1.1)	9.2 (1.3)		9.0 (1.3)

respectively. For reference purposes, the radial distribution function of the O-Cs + was analysed and a first coordination sphere radius cutoff of 4.0 Å was identified. From this radius, cutoff-based coordination numbers have been computed. All three algorithms agree that on average at least one water molecule is lost when Cs + is located in the first layer,

TABLE I :

 I Comparison of average Cs + ions coordination numbers determined by ASANN, SANN and a fixed cutoff (4.0 Å) for each defined region. The data is displayed in the form:

mean (standard deviation).

TABLE II :

 II Comparison of average Na + ions coordination numbers determined by ASANN, SANN and a fixed cutoff (3.2 Å) for each defined region. The data is displayed in the form: mean (standard deviation).close/not close) depending on a threshold value (here CN ASANN ).[START_REF] Spackman | Proceedings of the Sixth International Workshop on Machine Learning[END_REF][START_REF] Fawcett | Pattern Recognition Letters ROC Analysis in Pattern Recognition[END_REF] The ROC curve of a class of binary classifiers is constructed by plotting the True Positive Rate (TPR) against the False Positive Rate (FPR) computed on a data set, for different threshold values.

B. Cordero, V. Gmez, A. E. Platero-Prats, M. Revs, J. Echeverra, E. Cremades, F. Barragn, and S. Alvarez, Dalton Trans. , 2832 (2008).

ACKNOWLEDGMENTS

We thank Paul Clabaut, Mathilde Iachella, David Loffreda and Kamila Kazmierczak for providing us the geometries of nanorods and nanoparticles, bulk alloys and cobalt surfaces, respectively. The authors thank the SYSPROD project and AXELERA Pôle de Compétitivité for financial support (PSMN Data Center).

VI. SUPPLEMENTARY MATERIAL

Additional Theorems, Figures and Tables are available in the supplementary material.