Combining Simulations and Spectroscopy to Design Protective Coatings on Sulfide Solid Electrolytes

Introduction

All-solid-state batteries (ASSBs) have attracted significant interest owing to their potential to deliver higher energy density and improved safety compared with state-of-the-art Li-ion batteries. [1, 2, 3, 4, 5] In contrast to conventional Li-ion batteries, which rely on flammable liquid electrolytes, ASSBs employ solid electrolytes (SEs), thereby reducing safety risks associated with electrolyte leakage and flammability. Beyond safety, SEs can enable operation over wider temperature ranges, [1, 6] support faster charging by mitigating concentration polarization in the electrolyte, [4] and have the potential to enable higher energy and power densities through the use of Li-metal anodes and bipolar cell architectures. [7] Over the past decades, several classes of SEs with high ionic conductivities have been developed, including oxide, polymer, halide, and sulfide solid electrolytes. [8, 9, 10] Among these, sulfide solid electrolytes (SSEs) are considered particularly promising because of their high room-temperature ionic conductivity and their good processability. Within this family, argyrodite-based SSEs with the general formula $\mathrm{Li_6PS_5X}$ (X = Cl, Br, I) have received considerable attention due to their ionic conductivities of up to 5 mS $\mathrm{cm^{-1}}$ at room temperature. Despite these advantages, their adoption in industrially relevant applications remains limited. [11] Besides the low intrinsic electrochemical stability of SSEs, commercialization is also hindered by the pronounced moisture sensitivity of SSEs, which leads to the release of toxic $H_{2}S$ gas and degradation of the active material, even in dry rooms with low levels of water. This reactivity stems from the P–S framework, where the hard acid P readily reacts with the hard base O from atmospheric $H_{2}O$, replacing the softer sulfur atom. Various strategies have been explored to improve the moisture stability of SSEs, including the incorporation of softer acids to replace phosphorus and reduce susceptibility to attack by hard bases, as well as partial substitution of sulfur by oxygen. Recently, several groups have focused on coating SSEs with moisture-repellent layers, [12, 13, 14, 15, 16, 17, 18] in particular, octadecylphosphonic acid coatings proved to significantly improve the retention of ionic conductivity during dry-room exposure, while simultaneously improving electrochemical stability with both the anode and the cathode. [19]

Despite these efforts, the molecular-level principles governing coating effectiveness remain poorly understood. While a few studies include comparisons of different head groups, [13, 14] systematic investigations of head group chemistry are still lacking. Moreover, to the best of our knowledge, no studies have examined the influence of alkyl chain length. As a result, it is unclear which molecular features primarily control interfacial adhesion and which control moisture exclusion. [20]

Addressing these challenges requires an atomistic-level understanding of the organic–inorganic interface, including surface reconstruction, head-group anchoring, the conformation and organization of the adsorbed molecular layer, and the structure and dynamics of interfacial water. Such detailed insight, however, is difficult to obtain experimentally owing to the limited number of characterization tools and remains challenging to achieve using conventional simulation methods. Recent experimental studies have demonstrated that molecular coatings can substantially improve the stability of sulfide electrolytes. In these studies, [18, 21, 22, 23] static DFT calculations provide valuable atomistic insight into coating adsorption and interfacial interactions. However, the computational cost of DFT limits simulations to small systems and short timescales, preventing the study of collective finite-temperature phenomena such as surface reconstruction, molecular self-assembly at realistic coverages, and water penetration. These processes require large-scale simulations extending over nanoseconds and involving thousands of atoms.

The last decade saw the rise of machine-learning interatomic potentials (MLIPs), [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34] which hold the promise of overcoming this cost-accuracy trade-off. MLIPs, once trained on accurate quantum-mechanical reference data, reproduce

DFT accuracy at a cost only marginally higher than empirical force fields. [35, 36, 37, 38, 39, 40, 41, 42, 43] MLIPs have already provided detailed mechanistic insight into the bulk and surface properties of lithium thiophosphates and related solid electrolytes, including their ionic and thermal transport properties. [44, 45, 46, 47, 48, 49, 50]

More recently, foundation MLIPs [51, 52, 54] [53] appeared, such as PET-MAD, [51] which is built on the Point Edge Transformer architecture (PET) [55] and trained on the chemically diverse Massive Atomic Diversity (MAD) dataset. [56] PET-MAD proved to be sufficiently accurate to allow direct application to a variety of different systems, including SSE. [19, 51, 57, 58] Here we fine-tune PET-MAD to LPSCI surfaces and organic coating adsorbates using an $r^{2}$SCAN reference dataset (Sec. II A). Starting from a foundation potential pretrained across chemically diverse atomic environments provides a transferable representation that can subsequently be adapted to the specific LPSCI–organic interfaces considered here. This strategy enables the large-scale, finite-temperature simulations of coated surfaces and water interfaces that underpin this study. In this work, we present an integrated design study of molecular coatings for SSEs, combining targeted experiments with large-scale molecular dynamics simulations to reveal how coating molecules dynamically behave and interact with water at realistic interfaces. Figure 1 shows an overview of the studied systems. Based on our previous work showing promising performance of phosphonic acid (PA) head groups with long alkyl chains, [19] we expanded the head group selection to include also carboxylic acid (CA) and sulfonic acid (SA) to study how head group chemistry affects surface binding. Additionally, we selected alkyl chains with different lengths to investigate the effect on surface morphology. The molecules are denoted according to the alkyl chain length and head group, where Cn indicates the total number of carbon atoms in the alkyl chain and the suffix specifies the head group. Machine-learning-driven atomistic simulations establish the design rules step by step: we first identify the relevant LPSCl surface, investigate how different head groups interact with the surface, analyze how the adsorbed chains organize as a function of alkyl-chain length and surface coverage, and finally quantify how the resulting molecular layer protects the surface from water penetration and affects ionic conductivity. At every step, theoretical predictions are complemented with relevant experimental trends, such as conductivity measurements of the coated SSE, interfacial adhesion (FTIR and TGA), and the chain-length dependence of moisture protection (H₂S evolution), for carboxylic, sulfonic, and phosphonic acids of varying chain lengths.

Methods

Potential architecture and dataset generation

The MLIP architecture used in this work is based on the PET architecture [55] as implemented in metatrain. [59] Unlike many commonly used MLIP frameworks, [24, 25, 34, 60] the base PET model operates directly on the Cartesian coordinates of neighboring atoms, avoiding the need to construct symmetric invariant or equivariant descriptors. This design provides a favorable compromise between model expressivity and computational cost. Rotational invariance is instead learned during training through data augmentation. Recent work [61] showed that non-symmetrized PET can learn rotational symmetry and predict observables with a high degree of rotational equivariance. For further details of the architecture see the supplementary information (SI) or the relevant publications. [51, 55]

To obtain an accurate description of LPSCl surfaces and their interactions with coating molecules, we fine-tuned a PET-MAD model on a material-specific dataset containing 8058 structures. Fine-tuning was performed by low-rank adaptation (LoRA), [62] in which the pretrained weights are kept fixed and low-rank corrections to them are optimized.

Reference energies and forces were obtained from single-point DFT calculations using the Vienna Ab initio Simulation Package (VASP) code [63] with the renormalized Strongly Constrained and Appropriately Normed (r{sup}2SCAN) functional. {citep}[64, 65, 66] Gaussian smearing with a width of 0.08 eV was applied, together with D3-BJ dispersion corrections [67, 68]` and dipole corrections along the surface-normal direction for slab calculations.

The dataset contains bulk LPSCl and LPS, clean LPSCl and LPS surfaces, isolated molecular species, and combined surface–adsorbate structures. The initial bulk Li₆PS₅Cl structure was obtained from the Materials Project database. [69] Surface slabs were generated using pymatgen [70] by cleaving low-index crystallographic orientations with Miller indices up to 3. For each surface orientation, multiple terminations were considered by shifting the cleavage position along the surface normal. During slab construction, terminations that directly break P–S bonds within the PS₄³⁻ tetrahedra were avoided. A 5 nm vacuum region was added normal to the slab to suppress interactions between periodic images.

The molecular part of the dataset includes organic and inorganic adsorbates relevant to surface functionalization and degradation. Surface-adsorbate configurations were generated by initially placing the molecules above the LPSCl surface, followed by relaxation and molecular-dynamics sampling with the PET-MAD universal potential. [51] To increase the diversity of thiophosphate surface environments, structures from previous work on lithium thiophosphate surfaces were also included and recomputed using the same DFT reference protocol. [46]

The final dataset was split into training, validation, and test subsets using an 8:1:1 ratio. Model accuracy was assessed on the test set by comparing ML-predicted energies and forces with the corresponding DFT reference values. The fine-tuned PET-MAD model achieves an energy RMSE of 5.53 meV atom⁻¹ and a force RMSE of 78.17 meV Å⁻¹. For this system, fine-tuning a foundation model is clearly advantageous with respect to training from scratch, leading to faster training convergence and lower errors (compare § 0.8.1 for a detailed analysis).

Simulation protocol

All simulations used the fine-tuned PET-MAD potential. Molecular dynamics simulations were performed in LAMMPS [71] with a time step of 1 fs and a Nosé–Hoover thermostat. The bulk LPSCl cell was first relaxed at 300 K in the NPT ensemble to obtain the equilibrium lattice parameters. Surface slabs were then constructed using these bulk-relaxed cell parameters and equilibrated at 300 K by 200 ps of NVT dynamics, followed by energy minimization. This procedure allows Li redistribution and surface reconstruction while keeping the slab cell consistent with the relaxed bulk reference used for surface-energy evaluation (definition in SI Sec. S1.6).

Adsorption models were built by placing protonated and lithiated carboxylate, sulfonate, and phosphonate molecules of controlled chain length approximately 5 Å above the relaxed LPSCl $\{111\}$ surface (468-atom slab, lateral cell of $20.3 \times17.6$ Å$^{2}$), in the vacuum region, with the head group oriented toward the surface. The surface–adsorbate systems were relaxed by 500 ps of NVT dynamics at 300 K followed by minimization, after which adsorption energies were evaluated (SI Sec. S1.6). With 1, 9, or 18 adsorbed molecules, these systems contain 490–1310 atoms.

Water penetration was simulated by placing a pre-relaxed liquid water layer of 92 molecules approximately 5 Å above the coated surface, followed by 1 ns of MD simulation at 300 K (750–1590 atoms in total). Water molecules whose oxygen reached an O–Li distance below 2.5 Å were counted as penetrating the coating; the resulting areal density, normalized by that of the uncoated LPSCl surface, gives the relative water accessibility, which provides a dynamic accessibility metric rather than a reaction energy.

Surface-stability descriptor

To quantify how strongly the local surface coordination environment deviates from the bulk reference, we define an area-normalized coordination-deviation descriptor, $D$. This construction is motivated by previous works showing that surface stability can often be rationalized in terms of broken bonds, under-coordination, or local coordination descriptors, [72, 73, 74, 75] and recent work on lithium thiophosphate surfaces. [46] For LPSCl, however, a generic broken-bond count is not sufficient, because the relevant surface stability is controlled by whether chemically distinct anion environments preserve their bulk-like Li coordination. We therefore define $D$ as a bulk-referenced Li-coordination deviation around three anion environments:

Here, $\alpha$ denotes one of three chemically distinct anion environments: sulfur atoms belonging to $\mathrm{PS}_4^{3-}$ tetrahedra, isolated sulfide-like sulfur atoms, and chlorine atoms. $N_{\alpha}^{\mathrm{surf}}$ is the number of surface atoms of type $\alpha$, $\overline{\mathrm{CN}}_{\alpha}^{\mathrm{surf}}$ is the average Li coordination number of this surface atom type, $\overline{\mathrm{CN}}_{\alpha}^{\mathrm{bulk}}$ is the corresponding bulk reference value for the same chemical environment, and $A$ is the surface area of the supercell. Surface atoms were defined as S and Cl atoms within $d_{\mathrm{surf}} = 3.5$ Å of either vacuum-facing surface.

For reference, bulk LPSCl gives a time-averaged $D$ close to zero, with a standard deviation of $1.36 \times10^{-3}$ Å$^{-2}$. Surfaces with positive values of $D$ indicate a larger loss of bulk-like Li coordination, while negative values indicate local Li over-coordination.

Sulfur atoms were classified according to their local P coordination. A sulfur atom was assigned to $\mathrm{S}_{\mathrm{PS}_4}$ if it belonged to a $\mathrm{PS}_4^{3-}$ tetrahedron, identified by the presence of a neighboring P atom within the P–S bonding cutoff. Sulfur atoms without such P coordination were assigned to $\mathrm{S}_{\mathrm{iso}}$. This distinction avoids mixing the Li coordination environments of framework sulfur and isolated sulfide-like sulfur, which have different bulk coordination references. For each atom type, Li coordination numbers were calculated using a smooth distance-based cutoff, whose detailed definition is given in SI S1.2. The descriptor $D$ should therefore be interpreted as a relative, bulk-referenced coordination-deviation density rather than a direct count of broken bonds. A value of $D = 0$ corresponds to bulk-like average Li coordination in the analyzed surface region. Positive values indicate reduced Li coordination relative to the corresponding bulk environment, while small negative values can occur when surface relaxation brings mobile Li atoms closer to S or Cl sites and increases their smooth coordination weight. By separating $\mathrm{S}_{\mathrm{PS}_4}$, $\mathrm{S}_{\mathrm{iso}}$, and Cl contributions, $D$ directly measures whether the cleaved and relaxed surface preserves the chemically distinct local coordination motifs of bulk LPSCl.

Materials and sample preparation

All chemicals and solvents were used without any further purification and characterization unless otherwise specified. The argyrodite Li$_6$PS$_5$Cl was obtained from NEI Corporation and all handling of electrolyte powder was performed in an argon glovebox (O$_2$ & H$_2$O < 0.5 ppm). Propionic acid, hexanoic acid, octanoic acid, pentadecanoic acid, sodium octylsulfonate, and octylphosphonic acid were obtained from Sigma-Aldrich or Fluorochem. Octyl sulfonic acid was obtained by using an ion-exchange resin (Dowex 50W-X8) pre-treated with HCl, followed by low-temperature vacuum drying.

Prelithiated acids were obtained by dissolving the acid in double-distilled (DD) water followed by stoichiometric addition of 1 N LiOH to neutralize each acidic proton. LPSCl coated with different acids was prepared by wet mixing the corresponding acids at specific molar percentages relative to LPSCl, in order to keep the ratio between the acid head group and the electrolyte surface area constant. For 500 mg of LPSCl, the respective (Li-)acid was dissolved by stirring overnight in 2 mL of dry tetrahydrofuran (THF) (H$_2$O < 10 ppm, determined by Karl Fischer titration). LPSCl was then added, and the resulting suspension was stirred for 1 h to ensure uniform coating. For standard samples used in chemical characterization, the THF was subsequently evaporated overnight at 80 °C. To evaluate adhesion, the LPSCl–acid suspension was transferred to a centrifuge tube and centrifuged for 10 min at 4200 rpm, followed by the addition of 2 mL of fresh THF, which was repeated 2 more times to completely remove unreacted acid molecules.

Experimental characterization

Raman imaging was carried out using a WITec Alpha300R confocal Raman microscope (WITec, Ulm, Germany) equipped with Zeiss objectives (10× or 20×). The system was coupled via optical fiber to a 532 nm diode laser with a maximum output power of 125 mW. Samples were placed in airtight holders sealed with polyethylene film to prevent exposure during measurement. Fourier transform infrared - attenuated total reflection (FTIR-ATR) spectra were recorded on a PerkinElmer Frontier spectrometer or on an ALPHA II Bruker spectrometer inside the glovebox (O$_2$ & H$_2$O < 0.5 ppm) in the range of 4000-250 cm$^{-1}$. Thermal gravimetric analysis (TGA) was performed using Mettler Toledo TGA/DSC 3 using 40 µl aluminum crucibles between 25 °C and 600 °C at a heating ramp of 10 °C min$^{-1}$ under nitrogen gas. The ionic conductivity of LPSCl, both with and without additives, was evaluated by electrochemical impedance spectroscopy at 25

°C using a BioLogic VSP-300 potentiostat. Measurements were performed inside a glovebox over a frequency range of 1 Hz to 7.0 MHz with an applied AC perturbation of 10 mV. Approximately 30 mg of LPSCl powder was uniaxially pressed into an 8 mm diameter pellet using an ASC-T cell (Sphere Energy) under a load of 2 tons. Ionic conductivity was calculated as $\sigma= L/(RA)$, where $L$ is the pellet thickness (cm), $R$ is the measured resistance (Ω), and $A$ is the electrode area (cm$^2$). The resulting conductivity in S cm$^{-1}$ was converted to mS cm$^{-1}$ by multiplying by 1000.

H$_2$S evolution was measured using a home-built setup consisting of a plastic desiccator, a modular humidity generator (MHG32-ProUmid) and a H$_2$S sensor (MODBUS RS485). After stabilization of the relative humidity (RH) at 20% the H$_2$S evolution was measured.

Results

LPSCl surface stability

We first analyzed the relative stability of clean LPSCl surfaces generated from different crystallographic orientations using our fine-tuned PET model. For each Miller family, several terminations were constructed by shifting the cleavage position along the corresponding surface normal; Figure 2 shows the lowest-energy termination obtained for each orientation. The calculated surface formation energies range from approximately 0.57 to 0.68 J m$^{-2}$. Vertical and horizontal error bars denote the standard deviations of $\gamma$ and $D$ within the selected lowest-energy termination class, with both quantities averaged over 20 minimization frames.

Among the screened surfaces, the {111} termination exhibits the lowest surface formation energy. The corresponding Wulff construction predicts that the equilibrium morphology is dominated by the {011} and {111} facet families, which together account for approximately two thirds of the total surface area. This is qualitatively consistent with previous first-principles Wulff constructions of LPSCl, which also identified the symmetry-equivalent {110}/{011} family as a dominant exposed surface. [76] These low-energy surfaces show small values of the signed coordination descriptor $D$, indicating that the S–Li and Cl–Li coordination environments remain close to the bulk reference after relaxation (slightly negative values of $D$, as for {111}, reflect local Li compensation allowed by the smooth coordination-number definition; see SI S1.2). The structural snapshots in Figure 2 provide a microscopic interpretation of this trend: surfaces with low $\gamma$ preserve LPSCl-like surface clusters with approximately correct $\mathrm{Li}_6\mathrm{PS}_5\mathrm{Cl}$ stoichiometry (red dashed boxes), i.e., local structural units in which neighboring $\mathrm{PS}_4^{3-}$ tetrahedra form a framework supporting the Li, S, and Cl atoms. When the surface cleavage preserves these clusters without strongly breaking or twisting them, the local S/Cl coordination remains close to the bulk reference and the surface formation energy is low. In contrast, high-energy surfaces correspond to disrupted local cluster environments (gray boxes in Figure 2), through two distinct mechanisms: the $\{113\}$ termination cuts directly through an LPSCl-like cluster, whereas the $\{001\}$ surface fails to preserve its surface cluster during finite-temperature relaxation, as mobile Li atoms redistribute near the vacuum interface. Both mechanisms lead to a larger coordination deviation and a higher surface formation energy; a detailed discussion, including the role of finite-temperature Li mobility compared to static 0 K calculations, [46] is provided in SI S1.3. Because the $\{111\}$ surface has the lowest surface formation energy, retains a bulk-like local LPSCl cluster environment, and belongs to one of the two dominant facet families in the predicted equilibrium morphology, it is selected as the representative stable LPSCl surface for the subsequent adsorption-energy calculations and water-penetration simulations.

Head group chemistry and interfacial binding

We examine the binding of the coating molecules to the LPSCl surface, focusing on the combined effects of alkyl-chain length, head group chemistry, and lithiation state. Figure 3a shows the computed single-molecule adsorption energy as a function of alkyl-chain length for the carboxylic, sulfonic, and phosphonic acids (CA, SA, PA) and their lithiated forms (CALi, SALi, PALi). For all head groups, adsorption is favorable and increases with increasing chain length. This indicates that the alkyl chain also provides additional stabilization of the molecules through contact with the LPSCl surface.

Experimentally, coating was achieved using a wet approach similar to our earlier work (see § 0.2). [19] The FTIR spectra (Figure S6) show that molecular species remain associated with LPSCl after processing. Samples prepared from protonated and prelithiated additives display similar head-group regions after coating, which is compatible with a change in protonation state during processing, but does not by itself establish complete in-situ lithiation. Raman spectra (Figure S7) indicate that the main LPSCl vibrational features are retained after THF treatment and coating.

For the carboxylate series, Figure 3b shows a qualitative association between the calculated adsorption energy and the change in decomposition temperature relative to the corresponding additive alone. More strongly adsorbing carboxylates tend to show smaller decreases in decomposition temperature after association with LPSCl. However, the number of experimental points is limited, and the decomposition temperature also depends on the intrinsic thermal properties of the free and surface-associated species. The TGA comparison should therefore be regarded as qualitative support rather than as a direct experimental measurement of adsorption energy. The head-group-dependent TGA data are provided in Figure S8, and qualitatively corroborate the computed adsorption-energy ordering among CA, SA, and PA.

Within the tested surface models, the lithiated species generally show more favorable adsorption energies than their protonated counterparts. A possible explanation is that the Li atoms associated with the lithiated head groups can participate in coordination with both the head-group O atoms and nearby surface S/Cl sites, whereas the acidic H atoms mainly form more localized contacts with surface anions (S1 Section S2.5). Differences among the head groups may similarly reflect differences in their electronic structure, bond polarity, acidity, charge distribution, and oxygen-donor ability, in addition to the number and arrangement of available O coordination sites and the local adsorption geometries sampled at the surface. [77, 78] In particular, PA/PALi tends to show the most negative adsorption energies, while the differences between CA and SA are smaller. These trends should nevertheless be interpreted cautiously because the head groups can adopt different configurations and protonation states, and the available experiments do not isolate intrinsic adsorption sufficiently to establish a definitive head-group ranking. The calculations are therefore used to identify broad trends rather than a universal ordering of binding strength. Additional head-H/head-Li RDFs and representative adsorption configurations are provided in the Supporting Information (Figure S14 and Figure S13) as qualitative structural context.

The finite-temperature simulations also directly show spontaneous deprotonation events for the sulfonic acid head group, with the detached H relocating toward sulfide-like surface S sites. This observation is consistent with the comparatively high proton lability of the SA head group and provides a possible microscopic pathway for proton transfer during surface binding. Together with the similarity of the coated FTIR spectra obtained from protonated and prelithiated additives, it is compatible with deprotonation during coating or solvent processing. [19, 79] The fact that deprotonation occurs spontaneously over the course of short simulations provides strong evidence for a strong thermodynamic driving force, even though the statistics are not sufficient to determine the actual deprotonation rate. The analysis in SI Section S2.6 is presented as direct evidence of a deprotonation event within the simulated SA trajectories.

The ionic conductivity after coating with different head groups, protonated and lithiated, was measured experimentally and is shown in Figure 3c. The decrease in conductivity after THF treatment indicates that the solvent-processing step itself affects LPSCl, potentially through mild chemical interactions and/or solvent-induced changes in particle morphology that increase grain-boundary resistance. [80, 81, 82] Prelithiated additives generally retain higher conductivity than the corresponding protonated additives, while the conductivity trend does not follow the computed adsorption energies in a one-to-one manner. This lack of direct correspondence may reflect several coupled effects, including interfacial chemistry, amount of surface-bound organic molecules, particle contacts, and local Li mobility. The simulations provide a nanoscale picture of how head-group chemistry and lithiation state modify the interfacial environment underlying the experimentally observed differences in conductivity retention. For further analysis of the role of chain length and surface coverage, this study concentrates on CALi, as it provides the most favorable practical balance between calculated adsorption and conductivity retention among the tested head groups, rather than because it is the strongest-binding species.

Chain-length-dependent coating morphology

Building on the adsorption-energy analysis and the head group comparison above, we next examine how the adsorbed molecular layer organizes as a function of alkyl-chain length and surface coverage. We use lithiated carboxylates on the LPSCl {111} surface as a representative model system.

Figure 4 examines how the conformation of carboxylate chains changes with alkyl-chain length and coverage for lithiated carboxylate. The vertical molecular extension is quantified by $\langle z_{\mathrm{C,max}}\rangle-\langle z_{\mathrm{O}}\rangle$, where $z_{\mathrm{C,max}}$ is the height of the highest carbon atom in the alkyl chain, and $z_{\mathrm{O}}$ is the height of the anchoring oxygen atoms. For a single adsorbed $\mathrm{CH}_{3}(\mathrm{CH}_{2})_{n-1}\mathrm{COOLi}$ molecule, this height remains small and changes only weakly with chain length. This trend is illustrated by the bottom snapshots in Figure 4: both the short-chain (C2, bottom left) and long-chain (C14, bottom right) molecules lie nearly flat on the surface after 500 ps relaxation. Thus, increasing the chain length does not make an isolated molecule stand up; instead, the longer chain remains in contact with the substrate, consistent with the chain-length dependence of the adsorption energy discussed in Section III B.

The behavior changes as surface coverage increases. For nine adsorbed molecules, the molecular height increases with chain length, indicating that intermolecular packing begins to lift the alkyl chains away from the surface. The middle snapshots in Figure 4 show this transition directly: the short-chain molecules (C2, middle left) form a compact layer close to the surface, whereas the long-chain molecules (C14, middle right) form a thicker and more disordered adsorbed layer. For eighteen adsorbed molecules, the vertical extension increases much more strongly with chain length. This is consistent with the top snapshots, where the short-chain system (C2, top left) remains relatively thin, while the long-chain system (C14, top right) forms an upright, dense molecular layer. Therefore, the molecular orientation is controlled by both chain length and surface coverage. At low coverage, chain-surface contact favors flat-lying configurations, whereas at high coverage, chain-chain packing drives the chains into more upright conformations.

The effects of this coverage-dependent structure can be seen in the lower panel of Figure 4, which shows the experimental ionic conductivity of $\mathrm{C}_{n}$-CALi-coated LPSCl as a function of alkyl-chain length. A denser or thicker organic layer may improve surface coverage but can also introduce additional resistance to Li-ion transport across grain boundaries. The dependence on alkyl chain length is relatively weak but systematic: short-chain additives lead to stronger conductivity losses, likely due to stronger interaction with, or penetration into, the sulfide framework, as supported by FTIR and conductivity analysis (Figure S17 and Figure S18). In contrast, very long chains (e.g., C14) also reduce conductivity, which is attributed to hindered $\mathrm{Li}^{+}$ transport through the non-polar alkyl region and the higher effective mass loading at constant molar ratios. Among the four chain lengths tested (C2, C5, C7, and C14), ionic transport is best preserved for C7.

Moisture protection: water penetration and H2S evolution

Finally, we assess whether the coated layers protect the surface from moisture by combining the atomistic water-penetration simulations with the experimental degradation rate in a single comparison (Figure 5).

Computationally (Figure 5a), the coatings are compared by the relative water accessibility defined in § 0.2 B. Within the lithiated-carboxylate model, water protection depends strongly on chain length and surface coverage: once the head group is anchored, the morphology of the hydrocarbon layer becomes decisive for suppressing water access. At intermediate coverage, short and medium chains still allow water to reach the surface, whereas longer chains progressively reduce water accessibility by filling the interfacial gaps. At high coverage, the same suppression is achieved at shorter chain lengths. Notably, irrespective of surface coverage, the onset of the water protection in both cases (C14 and C7 at coverages of 9 and 18 molecules, respectively, see Figure 5) is at a coating thickness of approximately 8 Å (compare Figure 4 and Figure S23), provided that the chain length is sufficient to reach this thickness. Since coating thickness is jointly determined by chain length and surface coverage, these two parameters act as coupled descriptors of the protective layer. The simulations, therefore, identify chain length and surface coverage as the two variables controlling the continuity of the hydrocarbon barrier and the resulting accessibility of water to reactive surface sites.

Experimentally, the consequences of water penetration are reflected in the amount of $\mathrm{H}_{2}\mathrm{S}$ released into a closed chamber. Figure 5b shows the initial $\mathrm{H}_{2}\mathrm{S}$ evolution rates for 10 mol% coated LPSCl at 20% RH, together with a THF-treated uncoated reference, extracted from the time traces in Figure S24. Short-chain coatings remain close to this level, whereas longer chains progressively lower the initial $\mathrm{H}_{2}\mathrm{S}$ slope, in qualitative agreement with the simulated suppression of water penetration. The experimental reduction is less pronounced than the idealized simulations would suggest, as expected for powder-processed coatings that are likely to be inhomogeneous and may not reach the high-density, defect-free coverage represented by the highest-coverage simulation models. Moisture protection thus improves with increasing chain length, with C14 giving the lowest initial $\mathrm{H}_{2}\mathrm{S}$ evolution rate, whereas ionic conductivity is best preserved for C7 (Figure 4). Since protection sets in only once the coating reaches the thickness identified above and thicker layers increasingly hinder Li-ion transport, the most favorable balance among the tested chain lengths is expected for the shortest chains that reach this thickness, i.e. in the intermediate regime between C5 and C14.

with shorter chains as the coating density increases. The experimental H$_2$S measurements follow the same chain-length trend, supporting the water-exclusion mechanism revealed by the simulations.

Conclusions

In this work, we presented a framework that combined large-scale MLIP simulation with targeted experiments to establish molecular design principles for moisture-protective coatings on SSEs. For the Li$_6$PS$_5$Cl system, we identified the low-energy {111} surface as a representative interface for analyzing adsorption, coating morphology, conductivity retention, and water accessibility. We further demonstrate the utility of foundation MLIPs as priors for modeling complex atomistic systems, showing that fine-tuned models converge more rapidly and achieve lower energy errors than models trained from scratch on the same dataset.

The calculations show that interfacial anchoring depends on head-group chemistry and lithiation state, with lithiated species generally exhibiting more favorable adsorption than their protonated counterparts. Differences among the head groups are comparatively smaller and should be interpreted as broad trends rather than a universal binding-strength ranking.

Once the head group is anchored, the organization of the alkyl chains becomes critical for moisture protection. Isolated molecules remain largely flat on the surface, whereas increasing coverage drives intermolecular packing and progressively more upright, brush-like configurations. Longer chains and higher coverage reduce direct water access to the Li-rich surface, and the same level of protection can be achieved with shorter chains as the coating density increases. The experimental H$_2$S measurements follow the same chain-length trend, supporting the water-exclusion mechanism revealed by the simulations.

Importantly, moisture protection and ionic transport do not vary in the same way with chain length. While the water barrier improves progressively with increasing chain length and becomes significant at coating thicknesses of approximately 8 Å, the measured ionic conductivity shows a non-monotonic dependence with a maximum around C7. This suggests that intermediate chain lengths offer the best balance between moisture exclusion and Li-ion transport. This optimum corresponds to the shortest chains that still reach the protective thickness, while longer chains form thicker layers that increasingly hinder Li-ion transport.

Taken together, these results show that effective molecular coatings require a balance between interfacial anchoring, hydrophobic surface coverage, and preservation of Li-ion transport. The simulations provide critical insights for experimental trends that are otherwise probed only indirectly, and establish chain length and surface coverage as coupled design variables, acting through the coating thickness, for controlling the protection–transport trade-off. This study combines simulations and experiments to identify trends and mechanistic insights into a complex passivation process, and validates the use of fine-tuned universal MLIPs to enable atomistic simulations of unprecedented scale and complexity at first-principles accuracy. The weak dependence of H$_2$S evolution on coating suggests that future studies should focus on understanding the mechanisms of both water and Li transport across the coating, as well as the impact of defects and partial coverage. The availability of a validated MLIP provides the necessary infrastructure to facilitate these challenging joint modeling–experimental investigations, which have the potential to enable rational selection of passivation layers tailored to specific electrolyte compositions and processing conditions, accelerating the path toward commercially viable moisture-tolerant SSEs.

Data Availability

The data to reproduce the results in the paper, such as input files and the model used to perform the simulations, will be available on Materials Cloud [83] upon publication.

Acknowledgments

A.C., M.C., D.T., and L.B. acknowledge support from the Swiss National Science Foundation (SNSF) (Grant No. Sinergia CRSII5_202296). H.T. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) as a Walter-Benjamin Fellow (Project No. 534232519). L.B. thanks Prof. Salentinig from the University of Fribourg for the use of the Confocal Raman Spectrophotometer. This work was supported by grants from the Swiss National Supercomputing Centre (CSCS) under Project No. lp95.

References

Computational model validation and surface analysis

Machine-learning potential training and validation

The calculations were performed using machine-learning interatomic potentials based on the Point Edge Transformer (PET) architecture. [55] PET has demonstrated high accuracy across molecular and materials benchmarks. [51, 55, 58]

The model used in this work was initialized from the pretrained universal PET-MAD v1.0.1 checkpoint [51] and subsequently fine-tuned on the LPS/LPSCl dataset described in § 0.2 A. The network has two message-passing layers, each containing a two-layer transformer with an embedding dimension of 256 and 8 attention heads, and uses a cutoff radius of $4.5\,\text{\AA}$. Fine-tuning was performed by low-rank adaptation (LoRA) [62] with rank 8: the pretrained weights were kept fixed and only the low-rank updates (16 384 parameters, 0.42% of the model) were optimized. Training used batches of 850 atoms; the model used for all simulations is the checkpoint at epoch 1100, selected on the validation set as described below. Energy and force mean-squared errors entered the loss with weights of 0.1 and 1.0, respectively.

To assess the benefit of pretrained initialization, the same PET architecture was also trained from random initialization using the same material-specific dataset and evaluated on the same held-out test set. For this comparison, the model trained from scratch and the LoRA fine-tuned model were trained with identical data and settings (learning-rate schedule, batch of 850 atoms, cutoff, loss weights and random seed); the only difference is the initialization. Both models were trained for more than 1400 epochs, with a checkpoint saved every 50 epochs. For each training approach, the reported model is not the last checkpoint but the one with the lowest validation loss among all checkpoints up to epoch 1100 (epoch 900 for the model trained from scratch and epoch 1100 for the fine-tuned model). Note also that beyond 1100 epochs the model trained from scratch starts to overfit: its training error keeps decreasing faster than the validation error. On the other hand, the fine-tuned model is converged at this point: the validation energy and force RMSEs of all later checkpoints up to epoch 1500 differ from those of the selected checkpoint by less than 0.4 meV atom$^{-1}$ and 0.6 meV \AA$^{-1}$, respectively. As shown in Figure S1, fine-tuning reduces the energy RMSE from 7.26 to 5.53 meV atom$^{-1}$ (about 25%) and the corresponding MAE from 4.33 to 2.83 meV atom$^{-1}$. The force errors of the two models are of similar magnitude, with RMSEs of 82.90 meV \AA$^{-1}$ (scratch) and 78.17 meV \AA$^{-1}$ (fine-tuned) and MAEs of 45.6 and 47.8 meV \AA$^{-1}$, respectively. The small difference in force RMSE comes mainly from a few large errors (enlarged markers in Figure S1): the largest force-component errors of the model trained from scratch reach 7.4 eV \AA$^{-1}$, mostly on O atoms in a single test structure: an off-stoichiometric Li–O–P slab containing O–O bonded species, with DFT forces of up to 15 eV \AA$^{-1}$ on these atoms. This chemistry is represented by only 15 of the 6448 training structures, none with the same composition, and the fine-tuned model predicts all atoms of this structure with errors below 1.4 eV \AA$^{-1}$. Excluding the 20 atoms with the largest errors, the force RMSEs of the two models are 74.4 (scratch) and 76.1 meV \AA$^{-1}$ (fine-tuned). The fine-tuned model also converges faster and more smoothly: after 100 epochs its validation energy RMSE is 9.6 meV atom$^{-1}$, compared with 33.5 meV atom$^{-1}$ for the model trained from scratch, which needs about 500 epochs to reach a similar value (11.0 meV atom$^{-1}$; values averaged over 20 epochs). Overall, fine-tuning lowers the energy error, gives fewer extreme force outliers and converges much faster, while typical force errors remain of similar magnitude. This suggests that the representations acquired during pretraining on the broad MAD dataset [56] provide a useful prior for local environments that are only sparsely represented in the material-specific dataset, rather than requiring all relevant local interactions to be learned from scratch.

Coordination-number definition

Coordination numbers were computed using a smooth distance-based switching function. For an atom $i$, the coordination number with respect to neighboring atoms $j$ is

where $r_{ij}$ is the interatomic distance and

In this work,

Atoms below 0.5 Å therefore contribute with weight 1, atoms between 0.5 and 3.5 Å contribute with a smoothly decreasing weight, and atoms beyond 3.5 Å are not counted. The species-resolved coordination numbers used in Eq. (1) are

and

The same switching function was used for bulk and slab structures, so that

measures the surface coordination relative to the bulk reference for species $\alpha$.

Interpretation of the surface descriptor

For some terminations, including $\{111\}$, the signed coordination descriptor $D$ becomes slightly negative. This does not imply a negative number of broken bonds. Because the coordination number is smooth, Li atoms moving closer to surface $S$ or Cl sites during relaxation increase their fractional coordination weight. A small negative $D$ therefore indicates local Li compensation or slight over-coordination, and the relative trend across terminations is more meaningful than the absolute value.

The present protocol combines 300 K bulk equilibration, 300 K NVT slab relaxation, and final minimization. It therefore captures Li redistribution and surface reconstruction that are suppressed in static 0 K calculations. This interpretation is consistent with previous work showing that molecular-dynamics-induced reconstruction, surface density, and local coordination affect the stability of lithium thiophosphate surfaces. [46]

Surface termination selection

For a given Miller family, different cleavage positions generate distinct surface terminations with different local stoichiometries and coordination environments. Figure S2 illustrates this for the $\{111\}$ family, where representative terminations span a wide range of surface formation energies, from 0.568 to 1.064 J m$^{-2}$.

Because the equilibrium morphology is determined by the most stable surface available for each orientation, we use the lowest-energy termination of each Miller family in the Wulff construction and in the comparison of surface stabilities. Higher-energy terminations are therefore treated as metastable cleavage configurations rather than representative equilibrium facets.

Robustness of the Wulff morphology

To assess the sensitivity of the predicted equilibrium morphology to the set of surface orientations included in the Wulff construction, we progressively expanded the input from the three $\{001\}$, $\{011\}$, and $\{111\}$ families to all low-index surfaces with $\max(h,k,l) \le2$, and finally to the complete set of surface families considered in this work.

As shown in Figure S3, additional surface orientations progressively subdivide the equilibrium shape and reduce the individual area fractions of the dominant facets. Nevertheless, the qualitative morphology remains unchanged: the $\{011\}$ and $\{111\}$ families are the two largest contributors in all three constructions. In the full surface set, they account for 40.6% and 26.9% of the total surface area, respectively, corresponding to approximately two thirds of the predicted equilibrium surface.

Surface and adsorption energies

The surface formation energy of a slab with $N_{\mathrm{slab}}$ atoms and surface area $A$ (two equivalent surfaces) is

where $E_{\mathrm{bulk}}/N_{\mathrm{bulk}}$ is the energy per atom of bulk LPSCl. The adsorption energy of $N$ molecules is

where $E_{\mathrm{mol}}$ is the energy of the isolated molecule; negative values indicate favorable adsorption. All energies are averages over minimized MD frames.

Head-group chemistry and interfacial binding

Adsorption trends on low-energy LPSCl surfaces

To test whether the adsorption trends obtained for LPSCl $\{111\}$ are surface-specific, adsorption energies were also evaluated on the low-energy $\{011\}$ surface. Figure S4 reproduces the broad tendencies observed on $\{111\}$: adsorption generally becomes more favorable with alkyl-chain length, and the lithiated models tend to adsorb more strongly than the corresponding protonated models. Phosphonate-based species give the most negative values within the tested configurations. Because the head groups can adopt different geometries and protonation states, the quantitative separation between head-group families should be interpreted cautiously.

For comparison across the full chain-length range, Figure S5 summarizes the average adsorption energies of the protonated and prelithiated acids on LPSCl $\{111\}$.

Spectroscopic evidence for coating and phase stability

The ATR–FTIR spectra in Figure S6 provide evidence that molecular species remain associated with LPSCl after processing. Samples prepared from protonated and prelithiated additives show similar head-group regions after coating, while several features of the pristine protonated acids are suppressed or shifted. These observations are compatible with changes in protonation state during processing, but the spectra alone do not establish complete in-situ lithiation or quantify the fraction of deprotonated molecules.

Raman spectra show that the main vibrational features of LPSCl are retained after THF treatment and molecular coating (Figure S7).

Thermal stability and washing-based adhesion assessment

The TGA decomposition temperatures provide an experimental comparison of the thermal stability of the free and surface-associated additives. Figure S8 compares the different C7 head groups and the carboxylic-acid chain-length series. Prelithiated additives, particularly PALi, show high decomposition temperatures; however, these values also reflect the intrinsic thermal stability of the corresponding salts and therefore provide only qualitative support for the calculated ordering of adsorption energies.

Adhesion was also assessed indirectly by washing coated particles repeatedly with THF and quantifying the remaining organic content by TGA. The correlation with the computed adsorption energies is weak, particularly for short-chain acids (Figure S9), because the procedure can reflect kinetic removal and surface reactions rather than equilibrium adsorption. The washing data should therefore be treated as qualitative support rather than a direct adhesion measurement.

The raw C7 TGA traces below show the washing response for the three head-group families used in the main comparison. Figure S10 presents the carboxylic-acid series.

The corresponding sulfonic-acid and sulfonate traces are shown in Figure S11.

Figure S12 provides the equivalent phosphonic-acid and phosphonate data.

Representative adsorption configurations

Figure S13 shows representative relaxed structures for protonated and lithiated head groups on LPSCl {111}. The snapshots illustrate the range of local environments sampled by the different head groups, but individual configurations should not be interpreted as unique adsorption motifs or as direct evidence for the origin of the adsorption-energy differences.

Head-H and head-Li surface-contact RDFs

To provide additional structural context for the protonated and prelithiated models, Figure S14 compares peak-normalized RDFs involving $H$ or $Li$ initially associated with the molecular head group. The protonated systems show short-range $H$–$Cl$ features and weaker or broader $H$–$S$ features, whereas the prelithiated systems show head-$Li$–$Cl$ and head-$Li$–$S$ distances that overlap approximately with the corresponding bulk $Li$–$Cl$ and $Li$–$S$ distributions.

These RDFs indicate that protonated and prelithiated heads sample different local interfacial environments, but they do not by themselves establish specific bonding motifs or explain the calculated adsorption-energy differences. The curves are peak-normalized, represent ensemble averages over several adsorption configurations, and do not include $H$ atoms after they detach from the molecular head group. The comparison is therefore used only as qualitative structural information.

Finite-time deprotonation in the SA system

$H$ detachment was observed most clearly for the sulfonic acid head group (SA) within the finite MD trajectories. Figure S15 compares intact SA, a trajectory classified as deprotonated SA, and fully prelithiated SALi on the LPSCl $\{011\}$ and $\{111\}$ surfaces. In these sampled configurations, the deprotonated case shows intermediate adsorption strength and is not equivalent to fully prelithiated SALi.

Within the inspected trajectories, detached H atoms were found mainly near isolated sulfide-like S sites rather than S atoms within $PS_4$ tetrahedra (Figure S15). Figure S15 shows one representative event for the $n = 7$ SA molecule on LPSCl $\{111\}$. Because only a limited number of proton-transfer events were sampled, this analysis illustrates a possible pathway but does not provide a statistical deprotonation probability or establish the dominant experimental protonation state.

Chain-length and coverage effects

Effect of molecular coverage on adsorption energy

Figure S16 shows the adsorption energy per molecule of lithiated carboxylate chains on LPSCl $\{111\}$ at different coverages. Although the molecular conformation changes from flat-lying at low coverage to more upright at high coverage, the chain-length dependence remains similar. Increasing coverage therefore mainly shifts weak interactions from chain–surface to chain–chain contributions without substantially changing the head-group anchoring chemistry.

The larger error bars for the single-molecule systems reflect the broader range of adsorption sites and flat-lying configurations sampled by isolated molecules. At higher coverage, more sites are occupied simultaneously and the adsorption energy is averaged over several molecules, reducing the apparent variation.

Coverage estimates and experimental additive loading

To relate the simulated coverages to experimental additive loading, the LPSCl surface area was estimated by assuming spherical particles. Because the particle size and morphology are irregular, these values should be regarded as approximate. For particles of diameter $D$ and density $\rho= 1.64\ \mathrm{g\,cm^{-3}}$,

The molar fraction required for complete coverage was estimated as

Average particle size	Particle surface area	Specific surface area	mol% for full coverage
	$\mathrm{\mu m}^{2}$	$\mathrm{m}^{2}\,\mathrm{g}^{-1}$	mol%
0.1 $\mathrm{\mu m}$	0.0314	36.6	8.2
1 $\mathrm{\mu m}$	3.14	3.7	0.8
5 $\mathrm{\mu m}$	78.5	0.73	0.2
Literature solid*		2.8	0.6
Literature liquid-phase synthesis*		20.6	4.6

where $N_A$ is Avogadro’s number and $MW_{\mathrm{LPSCl}}$ is the molecular weight of LPSCl. Eighteen molecules on a 3.56 $\mathrm{nm}^{2}$ slab correspond to 0.20 $\mathrm{nm}^{2}$ per molecule, consistent with reported carboxylic head-group areas of 0.19–0.21 $\mathrm{nm}^{2}$. [85] The surface coverage of the experimental samples was not measured, so these values only indicate the order of magnitude of the loading required for a monolayer. More accurate estimates would require direct BET measurements.

Chain-length-dependent experimental characterization

The chain-length-dependent conductivity discussion in the main text is supported by the FTIR spectra in Figure S17, which compare the pristine carboxylic acids and the corresponding coated LPSCl samples across the investigated chain lengths.

The decrease in ionic conductivity with increasing additive loading, while chain length remains constant, supports the conclusion that the coating mildly hinders $\mathrm{Li}^{+}$ transport.

The raw TGA traces for the chain-length series are provided below. Figure S19 shows the short-chain C2 carboxylic-acid system. Figure S20 gives the corresponding intermediate-chain C5 data. The long-chain C14 carboxylic-acid system is shown in Figure S21.

Moisture protection

Water–Li radial distribution analysis

$O_{\mathrm{water}}$–Li RDFs were calculated for LPSCl $\{111\}$ surfaces coated with lithiated carboxylates at coverages of 9 and 18 molecules. A first peak near $r \approx2$ Å indicates direct short-range access of water oxygen atoms to Li-rich interfacial sites, whereas suppression of this peak indicates reduced water penetration.

As shown in Figure S22, short-chain coatings, especially $n = 2$ and $n = 3$, retain pronounced first peaks. The peak decreases strongly with chain length and is nearly absent for $n = 9$, $n = 11$, and $n = 14$. Suppression is stronger at 18-molecule coverage, consistent with a denser coating that more effectively limits direct water–Li contact.

Water accessibility versus vertical chain extension

To relate the water barrier to the height of the coating, Figure S23 shows the water accessibility (Figure 5) as a function of the vertical chain extension $\langle z_{\mathrm{C,max}}\rangle- \langle z_{\mathrm{O}}\rangle$ (Figure 4) at coverages of 9 and 18 molecules. At both coverages, the water accessibility decreases with increasing extension and vanishes at approximately 8 Å, which is reached by $n = 14$ at 9 molecules (7.8 Å) and by $n = 7$ at 18 molecules (7.7 Å). Below this extension, the denser 18-molecule coating gives a lower water accessibility than the 9-molecule coating at a similar extension, indicating that the packing density also contributes to the barrier.

Time-resolved H₂S evolution during humid exposure

Figure S24 shows the time-resolved H₂S concentration during exposure of THF-treated and Cn-CALi-coated LPSCl powders to humid gas at 20% RH. The THF-treated reference shows the largest H₂S accumulation, whereas all coated samples show reduced evolution. C14-CALi gives the lowest accumulation over the measured interval, consistent with improved moisture blocking by longer hydrophobic chains. The initial rates reported in the main text were extracted from these traces.