Active Learning and Neural Network Potentials Accelerate Molecular Screening of Ether-based Solvate Ionic Liquids

Solvate Ionic Liquids (SIL) have promising applications as electrolyte materials. Despite the broad design space of oligoether ligands, most reported SILs are based on simple tri- and tetraglyme. Here, we describe a computational search for complex ethers that can better stabilize SILs. Through active learning, a neural network interatomic potential is trained from density functional theory data. The learned potential fulﬁlls two key requirements: transferability across composition space, and high speed and accuracy to ﬁnd low-energy ligand-ion poses across conﬁgurational space. Candidate ether ligands for Li + , Mg 2+ and Na + SILs with higher binding aﬃnity and electrochemical stability than the reference compounds are identiﬁed. Lastly, their properties are related to the geometry of the coordination sphere

form a stable coordination complex where oxygen donors in the ligand coordinate the cation. This results in SILs with better thermal and oxidative stability than the two components alone and high ionic conductivity like a conventional IL. The properties of these SILs can be tuned by adjusting the ligand chemistry. [1][2][3][4] SILs synthesized with tri-and tetraethylene glycol dimethyl ethers (triglyme (G3) and tetraglyme (G4), respectively) and lithium bis(trifluoromethylsulfonyl)imide (LiTFSI), have been extensively studied in experiments 1,3,5,6 and molecular dynamics simulations, [7][8][9][10][11][12] and demonstrate the potential of SILs as a promising alternative to conventional ILs for electrolyte applications. Hence, finding SILs with desired melting temperature, oxidative stability, transport properties and cationligand association can advance battery technologies. We screen a combinatorial library of 10 3 ether molecules in search seeking to maximize ligand-ion binding energy, the electrochemical stability of the chelate, and the size of the complex. Neural-network-accelerated atomistic simulations allow identifying a number of synthetically accessible ligands with stronger affinities than the baselines and mapping the chemical space of oligoether-ion interactions for SILs. 3D  The exhaustive configurational sampling needed to evaluate every ion and ligand pair is computationally demanding, because the ion-ligand potential energy surface (PES) has many shallow local minima. Algorithms such as simulated annealing or seeding gradientbased optimization from different random guesses can perform such global optimization but are typically too costly to be combined with high-accuracy quantum chemical methods. Classical force fields are possible alternatives to the expensive ab initio methods for sampling.
However, they describe electrostatic interactions poorly and struggle to capture the interactions between ions and lone pairs. This results in inaccurate structure guesses which, in the best case, require extensive refinement with, and in the worse case, do not correspond with the global minimum of the quantum chemical PES. In order to combine higher accuracy and fast sampling, we trained neural network (NN) interatomic potentials based on graph convolutions (GC) 16,17 using density functional theory (DFT) energies and atomic forces for a set of reference chemistries. The GCNN potential was then deployed to quantify ligand-ion binding. An active learning (AL) strategy was used to acquire training data efficiently by carrying out DFT only on areas with high model uncertainty. 18,19 AL augments the data space and automatically improves the model by iteratively ac-quiring training data and deploying molecular dynamics (MD) simulations with the learned potential. Additional DFT simulations are carried out on points in the trajectory for which the model is uncertain. CGNN models trained with this approach showed high accuracy both in terms of energies and forces at the end of the AL loop, and also in predicting converged ion-binding structures that are local minima and agree with the DFT refined ones (Fig. 1). all the configurations used in the training data and colored the data by different generations ( Fig. 1 B), see SI for details). The configurations in the 0 th generation do not cover enough configurational space and hence showed poor performance when the learned potential is used to run MD (Figure 1 A). As the evolution progresses, the diversity in configurational space was systematically improved (higher dispersity of points in the projected 2D t-sne plot) in keeping with the validation accuracy of the model.
We then applied the trained NN models to sample the conformations of complexes between 1086 oligoether and each of the three ions. The low-energy geometries were refined using a long-range-corrected hybrid DFT functional, ωB97X-D3/def2-TZVP. 26 We first applied our workflow on well-studied crown ether binding with Li + and Na + . It was found that In the X:X:X:X short notation, : represents the oxygen and X is the number of carbons spacers or caps. Dashed vertical lines mark the energy of baseline glyme molecules (G3 and G4).
the calculated binding energies correlate well with experimental data measured in gas-phase and correlate reasonably with solvent data (see Table S2.  Figure 2, are predicted to be more stable than Li(G3) + and Li(G4) + baseline molecules by 0.2-0.3 eV while maintaining a similar E H level.
Li + binding can be improved by increasing the number of oxygen atoms and by increasing spacer length (Fig. 3). The longer carbon spacers facilitate the arrangement of oxygen atoms around the ion and therefore to the formation of favorable oxygen coordination structures.
In particular, having 3-C spacers is beneficial because it contributes to the formation of a six membered-ring with low strain energy. 30 Linear energy decomposition analysis suggests that larger rings contribute more to the binding energy ( Fig. S5 in SI). This hypothesis is further supported by geometrical analysis using τ , an index ranging from 0 to 1 that captures the deviation from an ideal coordination structure, a value of 1 corresponding to the lowest-energy coordination structure (tetrahedral and bi-pyrmidal for 4-and 5-coordinated structures respectively). 31 Strongly-binding ethers have high τ 4 and τ 5 . In comparison, the baseline G3 and G4 complexes are closer to square planar (4-coordination) and square pyramidal (5-coordination) geometries, resulting form the constraints imposed by the short spacers.
Because the hydrodynamic radius of SIL complexes is related to their self-diffusion, and hence to their conductivity, we have estimated the solvent-excluded volume of the converged ion-ligand complexes using a spherical probe with radius 1.9 A. Fig. S6 in SI reports how the larger ethers with longer spacers and more oxygen atoms and thus higher affinity are typically bulkier and are expected to be more viscous.
The HOMO level of the complex was used as a surrogate for oxidative stability (SI Figure. S7 reports the correlation between ionization potential and HOMO levels). As shown in Batteries based on sodium and magnesium are promising and earth-abundant alternatives for lithium, but they lack effective and stable electrolytes, so a similar analysis was performed for Mg 2+ and Na + complexes. Fig S4 reports   Experiments show relatively weaker binding in Na[Glyme] than the smaller lithium ion, which results in higher presence of uncoordinated glymes which compromise the stability of the electrolytes. 34 Our simulations are in keeping, the larger radius of sodium ions results in weaker ion-dipole interactions and lower affinity than Li + by about 1 eV. We thus focused on potential candidates from the G4 and G5 families (Fig. 2), and observed relatively minor improvement over the baseline.
Design of optimal ligands for a desired ion requires exploring a complex compositional and configurational design space. By combining fast and transferable machine learning with accurate DFT simulations through active learning, over 1,000 candidate oligoethers have been evaluated as ligands for Li + , Mg 2+ and Na + SILs. The tradeoff between affinity, electrochemical stability and size have been charted, and candidate ethers with better predicted properties than the baseline compounds have been identified.

Conflicts of interest
There are no conflicts to declare.

Data availability
The training data and optimized geometries are available for download at https://github.com/learningmatter-mit/.