Zachary A. H. Goodwin
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Malia B. Wenny
Department of Chemistry and Chemical Biology, Harvard University, Cambridge, MA 02138, USA

Julia H. Yang
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Andrea Cepellotti
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Kyle Bystrom
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Anders Johansson
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Lixin Sun
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Simon Batzner
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Albert Musaelian
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Jarad A. Mason
Department of Chemistry and Chemical Biology, Harvard University, Cambridge, MA 02138, USA

Boris Kozinsky
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA
Robert Bosch LLC, Research and Technology Center, Cambridge, MA 02142, USA

Nicola Molinari
John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA
Robert Bosch LLC, Research and Technology Center, Cambridge, MA 02142, USA

Chemical Transferability and Accuracy of Ionic Liquid Simulations with Machine Learning Interatomic Potentials

Ionic liquids (ILs) are an exciting class of electrolytes finding applications in many areas, from energy storage, to solvents for performing reactions in, to fluids for storing molecular gasses. Moreover, ILs have been touted as “designer” solvents, as they can be mixed to form complex fluids with tailored physiochemical properties. To make predictions of these properties, classical molecular dynamics (MD) have often been employed instead to simulate the dynamics of ILs, as it is computationally inexpensive. However, its accuracy is often under question, which motivates some to use density functional theory (DFT), but this is often prohibitively expensive, which imposes severe limitations on system sizes and time scales. Machine learning interatomic potentials (MLIPs) offer a way to bridge this simulation gap, offering a cheaper method than DFT but without significant loss in accuracy. As using MLIPs for ILs is still a relatively unexplored, several unanswered questions remain to see if MLIPs can be transformative in the space of ILs. As ILs are often not pure, but mixed together or with additives, we first demonstrate that a MLIP can be trained to be chemically transferable whilst only being trained on a few compositions, using a salt-in-IL as a test case. Having demonstrated that a transferable model can be trained, we investigate the accuracy of a ML IP for a novel IL trained on less than 200 DFT frames. To test the predictions of this model, we compare it against DFT and experiments, and find reasonable agreement. Overall, we hope these results can be used as guidance for future development of MLIPs for ILs, and complex mixtures of them, to make accurate predictions of these novel electrolytes.

§ INTRODUCTION

Ionic liquids (ILs) are a unique and highly promising class of electrolytes that contain no solvent in their neat form <cit.>. This lack of solvent, such as organic carbonates or water, makes ILs advantageous for numerous applications from supercapacitors/batteries <cit.>, to a medium to store gasses <cit.>, to themselves acting as “green” solvents for chemical reactions <cit.>. This is because ILs have an extremely low vapour pressure, are non-flammable and can withstand large voltages without decomposing <cit.>. These properties, in part, come from the ionic nature of the molecular species which comprise the IL <cit.>. These molecular ions are typically large and highly asymmetric to ensure the electrolyte remains liquid at room temperature <cit.>. Moreover, ILs can be mixed together to tune the desired physiochemical properties, creating a class of “designer” solvents with a huge chemical space <cit.>.

As the ions of ILs are often large and complicated molecular species, this has meant atomistic simulations are necessary to accurately and reliable make predictions of these complex electrolytes <cit.>. Typically, classical molecular dynamics (MD) has been used to simulate ILs and calculate physiochemical properties <cit.>, largely because of the computationally cheap and transferable force fields <cit.>. The accuracy of these predictions from classical force fields is often under question, however <cit.>. On the other hand, ab initio MD (AIMD) has also been used to simulate ILs <cit.>. However, the computational expense of running density functional theory (DFT) is prohibitively expensive in comparison, and limits the simulations to short time scales and system sizes <cit.>.

Machine learning interatomic potentials (MLIPs) have promised a way to bridge this efficiency-accuracy gap in atomistic simulations <cit.>. In the context of ILs, there have been several works which develop MLIPs for specific electrolytes to demonstrate this advancement, such as the work of Montes-Campos et al. <cit.>, Dajnowics et al. <cit.> and Ling et al. <cit.> for room temperature ILs, and the work of Tovey et al. <cit.>, Li et al. <cit.>, Mondal et al. <cit.> and Liang et al. <cit.> for molten salts. However, these works often employ large data sets, e.g., with 1,000$+$ DFT calculations, mainly focus on the NVT ensemble, and typically only for one neat IL or molten salt. With the recent development of equivariant MLIPs, such as NequIP/Allegro <cit.> and MACE <cit.>, training a MLIP for a specific application has become significantly more accurate and data efficient, i.e., only $\sim$100 DFT frames are needed. Hence, using equivariant MLIPs <cit.> to simulate ILs should provide an unprecedented accuracy/efficiency trade-off, but there are question which remain to be answered if these MLIPs are to transform IL simulations.

In this work, we aim to answer if a chemically transferable MLIP can be trained, and the accuracy of such a model on a novel system with relatively little training data. For a common salt-in-IL ( in ), we demonstrate that  can learn to be chemically transferable, whilst only being trained on several compositions, but only if these compositions are carefully chosen. Of the routes tested, we find that at least 2 compositions are required, and the best approach is to train on high-entropy Li-salt molar fractions, e.g. 0.4-0.6, as these compositions ensure all interactions are well sampled. Secondly, we synthesize and experimentally characterize a novel ionic liquid, , and train a chemically accurate Allegro model with less than 200 DFT training snapshots, capable of capturing experimental properties of this , with the main limitation being the approximations of our DFT. For this system we also investigate the role of how the structures are generated, and perform a hyperparameter search to optimise the model.

§ METHODS

§.§ Computational

§.§.§ Classical Molecular Dynamics

We use classical molecular dynamics (MD) simulations to generate a diverse set of uncorrelated structures as training data. We study two ILs here, and as the details of the MD simulations are almost identical, we will describe our method and point out where the ILs are different. To generate initial structures, we start by creating a low-density, non-overlapping initial structures by placing the ions randomly on the vertices of a three-dimensional cubic grid <cit.>. To overcome local energy barriers, in search of a lower energy minima, we then use an equilibration routine that comprises of a set of energy minimizations, compression/decompression and annealing stages, broadly based on previous works <cit.> (see  for more details).

All MD simulations have periodic boundary conditions, use a timestep of $\delta t = 1.0$ fs, a velocity-Verlet algorithm to evolve the equation of motions, a Nosé-Hoover barostat (1000$\delta t$ coupling) and thermostat (100$\delta t$ coupling) to enforce pressure and temperature, respectively <cit.>, and are performed using the LAMMPS code <cit.>. The atomic interactions are modeled with OPLS-AA <cit.>, and the parameters obtained from Ref . For the salt-doped IL,  in , to reproduce realistic neighbour densities, we rescale the atomic charges to 80% of the original value, as previously reported <cit.>. Whereas, for the  we rescale the atomic charges by 70% of the original value to match the experimental data better (see  for more details).

For  in  we investigate various Li-salt molar fractions, ranging from 0.0 to 1.0, i.e., from a pure IL to a pure salt, respectively. Each composition has 20 ion pairs in the simulation box, which means that the number of atoms at each composition is different, ranging from 520 at 0.0 Li-salt molar fraction to 160 atoms at 1.0 Li-salt molar fraction. The number of anions remains at 20 for all molar fractions, but the number of each cation changes by 2 for each adjacent molar fraction. For each molar fraction, two initial structures were generated, one for training and validation, and another for testing the trained models. The same number of frames and atoms are in each data set. To generate the data sets, each structure is evolved for 80 ns in a constant number of particles, volume, and temperature (300 K) ensemble. The atomic positions and force components are saved every 0.02 ns, to generate a total of 4000 snapshots for each composition.

To generate a diverse set of uncorrelated structures for , we simulate three different temperatures, from 300-600 K in steps of 150 K, in the NPT ensemble with a pressure of 1 atm, with 10 ion pairs in the simulation box, corresponding to 530 atoms. For every new temperature, we first let the system equilibrate for 1 ns, and then run for 6 ns during which snapshots of the structure are saved every 20 ps to ensure sufficient decorrelation. Thus, we generate a total of 900 snapshots, 300 per temperature.

§.§.§ Density Functional Theory

To obtain accurate energy, force and stress data, we performed density functional theory (DFT) calculations using Quantum Espresso <cit.>. We utilise the PBE exchange-correlation functional <cit.> and D3 <cit.> dispersion correction for van der Waals interactions, and use the ultra-soft pseudo-potentials of Ref. . We use $\Gamma$-point sampling, a cut-off of 50 Ry for the wavefunctions, a cut-off of 500 Ry for the charge density, and a convergence threshold of 10$^{-10}$ Ry.

For , we initially selected 60 frames at each temperature to generate a small DFT training data set. The force and energy values from DFT from these OPLS sampled structures appeared reasonable, examples of which are shown in the , but the diagonal components of the stress tensor were of the order of $\sim$1 GPa. This suggested there could be an issue with the DFT settings, such as the employed PBE exchange-correlation functional not localising the charge. To test the origin of this systematic stress, we performed a number of checks. First, we turned off the D3 correction and performed a variable cell BFGS relaxation with to ensure the system remained condensed. We found the system reached a density of $\sim$1.8 gcm$^{-3}$ after 50 steps, which is slightly higher than the OPLS density and the experimental value (see ). Therefore, the ions appeared to still bind in the absence of a D3 correction, which suggested they remained at least partially charged.

During the course of the variable cell relaxation test, it was noticed that the stress dropped rapidly, even without the cell changing size dramatically. This suggested the origin was not from the cell size, but some structural features of the ions, perhaps from the inconsistency between the potential energy surface of OPLS and PBE+D3. With this in mind, we performed fixed-cell AIMD on each structure for 50 steps (with a $\sim$1 fs timestep) at the corresponding temperatures (using the Verlet algorithm and a tolerance for rescaling velocities of 10 K) and also 10 steps of fixed-cell relaxation for each structure, the results of which are shown in the . For the AIMD trajectories, the energy drops by $\sim$10-20 eV per box (corresponding to $\sim$20-40 meV/atom), the pressure of the box fluctuates, but is distributed around a value close to 0 at the end of the short run, all whilst the forces oscillate, but overall do not dramatically change in magnitude. We also performed 10 steps of fixed-cell relaxations on these starting OPLS-AA frames, and found the energy reduced by 20-40 eV, the pressure reduces close to 0, and the forces are naturally close to 0. Each of these data sets is used to train MLIPs, to see if there is any difference between how the data is generated, the details of which are shown in the . Based on these data-sets, we also perform an iterative training loop to further optimise these MLIPs. We also performed some longer AIMD runs with variable and fixed cell to obtain indications of the equilibrium density at 300 K and to ensure the energy/stress was roughly equilibrated after the 50 AIMD steps.

For  in  we selected 50 frames from each composition and run, 100 in total for each composition, and performed DFT calculations for the energy, force and stress. We also calculated the single atom energies of all the elements in these frames, from only keeping each element in a box with length of $\sim$30 Å.

§.§.§ Machine Learning Interatomic Potentials

We use the state-of-the-art graph neural networks,  and Allegro, to train, validate and run the IL ML IP <cit.>. We refer to readers to Refs.  for an extensive explanation of the network architecture and the importance of equivariant models <cit.>.

For the classical force field data set of  in to perform the chemical transferability test  with  , we use $r_{max} = 5$, $l_{max} = 1$ with irreps 16x0o + 16x0e + 8x1o + 8x1e, with 5 layers, and use a learning rate of 0.05 with a batch size of 5. For each composition we use 600 for training and 50 frames for validation, and when multiple compositions are trained on these data sets are added together to increase the size of each training/validation set. We only train on energy and forces, and do not change the scale of these in the loss function. The test set was an independent 650 frames at each composition.

For the small DFT training data set for  in , we employ Allegro <cit.> to perform the chemical transferability test. We choose $r_{max} = 6$, $l_{max} = 2$, $n_{layers} = 2$, learning rate of 0.0005, batch size of 1, 80-20% split of training validation data (for each model, 50 frames are included for each composition, so the models trained on 1-3 compositions have 50-150 frames in the training data), and find that it is essential to employ single atom energies when training on multiple compositions with DFT data from Quantum Espresso. We do not vary these hyperparameters, as we are not optimising the model to be the best, but are only looking at how energy and force errors propagate across the composition range. The test set was an independent set of 50 frames at each composition.

The results shown for  in the main text are based on a single model (see  for details). We choose $r_{max} = 7$, $l_{max} = 2$, $n_{layers} = 2$, LR = 0.001, batch size of 4, stress weight of 100 (with force and energy of 1), and 80-20% split of training validation data which consists of 180 frames taken from the last frame of the short AIMD trajectories (starting from the OPLS sampled frames) and an additional 40 frames collected from performing volume scans with an Allegro model trained only on the 180 frames (see  for details).

To run these models, we employ LAMMPS with a custom pair style. We use the Nosé-Hoover barostat (1000$\delta t$ coupling) and thermostat (100$\delta t$ coupling) to enforce pressure and temperature with a timestep of $\delta t = 0.5$ fs. To obtain the isothermal compressibility, we first equilibrated the system at 1 bar for 50 ps. Then the pressure was increased by 1000 bar and equilibrated again for 50 ps before running for another 50 ps at this pressure for production; this process was repeated 18 times to obtain how the equilibrium volume changed with pressure, which was then fit to the Tait equation of state to obtain the isothermal compressibility <cit.>. Details of the other settings used for each calculation will accompany the results shown.

§.§ Experimental

§.§.§ Synthesis and Characterization

We synthesized  according to Ref. , which relies on a bicarbonate intermediate reacting with perfluoropentanoic acid to form the desired IL, with concomitant formation of CO$_2$ and H$_2$O (Figure S1). Briefly, 1,1,1,2,2,3,3,4,4,5,5,6,6-tridecafluoro-8-iodooctane and imidazole were reacted to form 1-(3,3,4,4,5,5,6,6,7,7,8,8,8-tridecafluorooctyl) imidazole. This precursor was methylated with methyl iodide to form 1-(3,3,4,4,5,5,6,6,7,7,8,8,8- tridecafluorooctyl)-3-methylimidazolium iodide. A bicarbonate-loaded anion exchange column was then used to generate 1-(3,3,4,4,5,5,6,6,7,7,8,8,8-tridecafluorooctyl)-3-methylimidazolium bicarbonate as an un-isolated intermediate. Reaction of this intermediate with perfluoropentanoic acid in water resulted in the formation of 1-(3,3,4,4,5,5,6,6,7,7,8,8,8-tridecafluorooctyl)-3-methylimidazolium perfluoropentanoate () as a pale yellow liquid. Complete experimental details, including 1H and 19F NMR data consistent with the formation of the target IL, are included in the supporting information of Ref. .

A residual iodide content of 1492 ppm was measured via ion chromatography by Galbraith Laboratories. A residual water content of 3800 ppm was measured via Karl Fisher titration on a Mettler Toledo C20 Coulometer. Note that the Karl Fisher titration was performed on a sample that was dried under vacuum for $\sim$5 h at ambient temperature, and the sample had to be exposed to air prior to performing the titration. Prior to further physical property measurements, the IL was dried at 50 $^{\circ}$C under vacuum for at least 12 h and was subsequently stored in a dry glovebox.

§.§.§ Variable-Temperature Density Measurement

The density of  was measured using an Anton Paar DMA 4100 variable temperature density meter, which has an accuracy of $\pm$0.0001 gcm$^{-3}$. Prior to each measurement, the instrument calibration was verified at 20 $^{\circ}$C and ambient pressure. Briefly, the density of ambient lab air and the density of nanopure water were recorded and compared to the expected values. A verification was deemed successful if the density of air was 0.0012 $\pm$0.0002 gcm$^{-3}$ and the density of nanopure water was 0.9982 $\pm$0.0002 gcm$^{-3}$ at 20 $^{\circ}$C. After calibration verification, the sample cell was briefly dried using a built-in air pump to prepare for sample injection. For each measurement, a pre-dried sample (1.1 mL) was loaded into a 3 mL syringe in a dry glovebox and temporarily protected with a N$_2$ headspace. The syringe was brought out of the dry glovebox, and the sample was immediately injected into the sample cell, with the syringe left in the injection port to prevent movement of the sample or the entry of gas bubbles. The density of the sample was recorded at 2.5 $^{\circ}$C steps between 2.5 $^{\circ}$C and 70 $^{\circ}$C. The sample was then carefully removed using the syringe, and the sample cell was washed with at least 20 mL of nanopure water and dried with an air pump to prepare for the next measurement.

§.§.§ Isothermal Compressibility Measurement

As described in our previous work <cit.>, isothermal compressibility values were obtained using small-angle X-ray scattering (SAXS) measurements performed at beamline 12-ID-B at the Advanced Photon Source at Argonne National Laboratory. Samples were probed using 13.3 keV X-rays, and the sample-to-detector distance was calibrated using a silver behenate standard. Scattered radiation was detected using a Pilatus 2 M detector. Beam defocusing and short exposure times (0.1 s) were used to reduce the chance of sample damage by the beam. The sample holder described previously, which consisted of a quartz capillary glued into a metal ring and connected to tubes to allow injection of the liquid sample, was used to obtain high-quality background and calibration data. Specifically, after measurement of the empty capillary as a background, a measurement of the SAXS pattern of nanopure water was used to obtain an absolute intensity calibration. Subsequently,  was injected into the sample holder and its SAXS pattern was recorded (Figure S2).

The zero-angle scattering intensity of the SAXS pattern was obtained as the y-intercept of a linear fit to the low-angle, featureless region between 0.08 and 0.16 Å$^{-1}$. The known isothermal compressibility of water (4.524 $\times$ 10$^{-10}$ Pa$^{-1}$) <cit.> was used to obtain a scaling factor, $a$, using an expression that relates the zero-angle scattering intensity, $I(0)$, to the isothermal compressibility, $\beta_T = aI(0)/(\eta^2k_BT)$, where $\eta$ is the scattering length density obtained using the molecular formula and density of the sample at 25 $^{\circ}$C, $k_B$ is the Boltzmann constant, and $T$ is temperature in K. This expression, with the value of  obtained from water, can be then used to provide the isothermal compressibility of an unknown liquid such as .

§ RESULTS AND DISCUSSION

§.§ Salt-in-Ionic Liquids - Chemical Transferability

Some of the most valuable properties of ILs are their ability to be mixed together to form “designer fluids” with physiochemical properties that can be precisely tuned <cit.>, and to act as “universal solvents” for reactions or to store gasses <cit.>. Therefore, being able to simulate complex mixtures of ILs and solutes is a necessity. Classical MD and AIMD are highly transferable in chemical space, allowing complex IL mixtures to be studied <cit.>, but they have limitations for IL simulations, as described in the Introduction. Typically, a MLIP is trained and tested on for a specific composition of a mixture, and the potential is often not used outside of the composition range it was trained/tested on. One question that needs to be answered in the context of MLIPs, which has received little quantitative attention, is if a MLIP can be trained to be chemically transferable, and still be accurate enough to reliably use for calculations on compositions that were not trained on. If it can be demonstrated that there is a method to train a MLIP on only few compositions, and it reliably transfers to the entire composition range of interest, this will open up the field of ILs to be accurately simulated with MLIPs. We aim to find out if this is possible, and what composition sampling method is best to achieve this goal.

To cheaply answer this question, we first adopt classical force fields to generate the data to train/test on with  (see Methods for details). While clearly the final model will not possess a greater accuracy than that of the classical force field, here we do not aim for accurate models, but rather to understand the intricacies of applying graph-based machine learning  to complex, molecular ionic systems, such as ILs. This practice allows us to perform a large number of preliminary studies without the high computational cost associated with generating an accurate training set from DFT. Moreover, as the potential energy surface from classical force fields is presumably simpler than DFT, this should set an upper bound on the chemical transferability, i.e., we expect the chemical transferability of a  model trained on DFT data to be worse.

With this goal, we choose a well-known salt-in-IL (SiIL),  in . This SiIL is of interest in energy storage applications because the IL is non-flammable and acts as a medium for the active cations to be transported <cit.>. This SiIL has recently attracted attention because the Li$^+$ has been observed to have a negative transference number at low salt molar fractions, as has been well studied from classical MD <cit.> and experiments <cit.>.

Chemical transferability investigation of  for a typical salt-in-IL,  in . Each row corresponds to a different  model trained on the molar fractions shown by the green triangles. Each number at each molar fraction is the mean absolute error (MAE) on an independent test set for forces, as shown on the left (a) in units of meVÅ$^{-1}$, and energy, on the right (b) in units of meV/atom.

In Fig. <ref> we display the results for our chemical transferability investigation. Each row represents a different  model trained on specific molar fractions (details in the Methods and ), as indicated by the green triangles, with the values at each molar fraction (in each box) being the test error on an independent data set (not in the training/validation set). On the left is the force mean absolute error (MAE) in meVÅ$^{-1}$ and on the right is the energy MAE in meV/atom. In what follows, we will describe the results and implications of Fig. <ref>, starting from the bottom row and working upwards.

If a  model is trained on a single composition at a Li-salt molar fraction of 0.5, we find the force MAE to be 8 meVÅ$^{-1}$ and the energy MAE to be 1 meV/atom at the same mole fraction. As we are using a large amount of data from a classical force field, these error values are considered as the baseline to measure how the test error differs for compositions not in the training data, and we will not consider their absolute values to be significant. For compositions further from the 0.5 molar fraction in the training data, the force MAE remains relatively constant, but the energy MAE increases to values $100-400$ meV/atom at the ends of the composition range. Therefore, training a mixture with only one composition works well for the composition trained on, as is often done for MLIP, but the  model becomes worse the further out-of-distribution the composition becomes.

Another approach could be to train on a couple of compositions at one end of the range, in the hope that several compositions allow better generalisability. The two rows above are examples of this strategy, trained on 2 and 3 compositions, at 0,0.1 and 0,0.1,0.2 Li-salt molar fraction, respectively. Again, the force MAE's are relatively constant across the entire composition range, albeit starting to increase to values larger than 10 meVÅ$^{-1}$ at 0.8-1. This increase in the force errors can be traced back to the increase in force errors of , with the forces of Li$^+$ becoming significantly worse. This is perhaps not surprising as the forces of  have not been sampled well at this low salt molar fraction. For the  model trained on 2 compositions, the energy MAE increases to 30-50 meV/atom at 0.8-1 molar fraction, which is better than the model trained on a single composition, but still corresponds to a $\times 50$ increase. The model trained on 3 compositions between 0-0.2 has similar trends, but the errors are reduced in magnitude, and only correspond to a $\times 6$ increase in energy MAE at 1.0 Li-salt molar fraction.

Finally, if we train on 2 compositions, but not at a one end of the composition range, we can obtain energy and force MAEs that are typically within $\times 2$ of the MAEs of the compositions trained on. Training on the compositions of 0 and 1 gives energy MAEs within $\times 2-3$ over the entire composition range, and the force errors also remain within this limit. The increase in the force error for this model can be traced back to an increase in the error of all elements, which could be rooted in the fact that no Li$^{+}$-[EMIM]$^+$ interactions are included in the data set. Similarly, training on molar fractions of 0.4 and 0.6 have force and energy MAEs which are roughly constant across the entire composition range, i.e., within $\times 2$ (apart from 1.0 Li-salt molar fraction compositions).

Overall, we are able to train a chemical transferable  model that has errors on compositions not trained on that are typically no more than $\times 2$ of the errors of the compositions that are trained on. The best approach, at least for this SiIL, appears to be when at least two compositions are trained on, which are not adjacent to each other and maximise the sampling of all interactions between species. This is demonstrated by the lowest and most consistent errors for the  model trained on molar fractions of 0.4 and 0.6.

To ensure these conclusions are robust, we performed a similar analysis where we took only 50 frames from the training/validation set and 50 frames from the test set and calculated the ab initio energies and forces, and used Allegro to train on. The results of this analysis are shown in the . We find the conclusions drawn from the classical force field are robust, with the chemical transferability test on only a small amount of DFT data to be qualitatively identical. Therefore, such chemical transferability test can be performed with classical force fields first while sampling uncorrelated structures, before proceeding with expensive DFT calculations. In the  we show further examples of Allegro trained on different combinations of compositions.

and Allegro learn directly on the force and energy values, with forces being per atom quantities and the energy being the total energy of the simulation box. We find that the forces predicted by  and Allegro generalise well, even if only 1 composition is sampled, provided interactions between species have been well sampled. This can be achieved from sampling from the middle of the composition range (0.4-0.6), or “high-entropy” configurations, but if we only sample from ends of the composition, the training data does not include enough cation-cation interactions, for example.

In contrast, the generalisation of the energy appears to be significantly more dependent on how compositions are selected. To understand these observations, and generalise the results, we perform a simple analysis. Let's assume that there is just a single training frame for each composition. For the 0.5 composition, let's assume the total energy can be written

\begin{equation}
E_{0.5} = \mu_{c1}N_{c1,0.5} + \mu_{c2}N_{c2,0.5} + \mu_{a}N_{a,0.5},
\end{equation}

where $\mu_j$ and $N_j$ are the energy per ion and the number of ions, respectively. Note this simple, linear equation assumes no ions interact, and just tries to describe the fluid as an ideal mixture, but this can still provide us insight into the energy generalisation. The $N_j$ values are known from the training data, but we are using /Allegro to find the $\mu_j$'s. As we only have 1 equation but 3 unknowns we are under-determined and no unique solution exits. Therefore, /Allegro finds values for $\mu_j$ which satisfy $E_{0.5}$, but these are not necessarily the correct values for the energy per ion, and therefore, we have not learnt the single-ion contributions to the total energy, which means the energy generalises very badly for this case (as more clearly see in the  for the DFT data).

Instead, for example, we train on the best case of 0.4 and 0.6, and because $N_a = 20$ in all frames from the salts sharing the same anions, we can write the equations as

\begin{equation}
E_{0.4} - 20\mu_{a} = E_{0.4}^{'} = \mu_{c1}N_{c1,0.4} + \mu_{c2}N_{c2,0.4},
\end{equation}

\begin{equation}
E_{0.6} - 20\mu_{a} = E_{0.6}^{'} = \mu_{c1}N_{c1,0.6} + \mu_{c2}N_{c2,0.6}.
\end{equation}

From these equations, we can determine the energy of $\mu_{c1}$ and $\mu_{c2}$ up to some arbitrary constant that is incorporated in $\mu_a$, which still ensures the total energy of the composition is correct, even though the individual ion energies are not. If fact, this means the only quantity which is well defined is the energy difference between the cations. Therefore, with only 2 compositions to train on, /Allegro can learn how the energy should change with composition reasonably well. This logic explains why training on 2 compositions instead of 1 significantly helps, but training on 3 instead of 2 does not significantly help with learning the energy to be generalisable.

If this is the case, then we should find that the predictions of the atomic energy of Li$^{+}$, for example, should not be the same across different models trained, but that the energy difference between Li$^+$ and [EMIM]$^+$ should be model independent (if trained on 2 compositions). We confirm with Allegro on DFT data that this is the case (at least approximately), as shown in the . Moreover, we find the models trained on 0.4,0.6 and 0.0,1.0 have almost model-independent values for the energy of each ion, suggesting they have found the “correct” contributions for each ion with only 2 compositions. Whereas, the models trained on one end of the composition range, even when training on 3 compositions, have only learnt the energies up to some constants.

There has recently been significant efforts to develop general MLIPs which work for most biological systems <cit.>, metals <cit.>, crystals <cit.>, or more recently foundational ML IPs trained on the Materials Project data <cit.>. These models have had remarkable success in their applications, but limitations of the models have also been identified <cit.>. We hope the quantitative approach taken here can be further developed and applied to these larger models, to reveal their limitations and permit them to be improved.

§.§ Novel IL - Accuracy

Many ILs are excellent “green” solvents for reactions, owing to their ability to dissolve many molecular an ionic species, their extremely low vapour pressure, low flammability and wide window of operating temperatures/voltages <cit.>. One property in this area that requires extremely accurate energy predictions is the solubility of species, as any error in the energy change is exponentialized <cit.>. For example, calculating the solubility of O$_2$ or CO$_2$ in ILs requires the these species to be inserted into the IL and the energy change calculated. Classical force fields have typically been used to compute solubilities, but the energy of insertion of these molecular species is not at the same level of accuracy as DFT, which means the solubility computed from classical force fields is sometimes erroneous <cit.>. MLIPs promise a way to bridge this gap: they are DFT-accurate, but are still computationally cheap enough to perform these calculations <cit.>.

To work towards this goal, we synthesis and characterise (details in the Methods) a novel IL that is of interest for its O$_2$ solubility <cit.>. Specifically we investigate , the structure of which is shown in Fig. <ref>(a). This IL is similar to several others studied in Ref. . As can be seen in Fig. <ref>(b), this IL has a highly heterogeneous, lamellar-like structure, which creates voids that allow dissolution of O$_2$ <cit.>. We aim to develop and test the accuracy of an Allegro model for , whilst only using a small amount of data and only performing 1 iterative training round at most (as it is expensive to keep training deep learning models), the details of which are given in the Methods and . We validate the predictions of this model against DFT and also the experiments. In Ref. , the solubility of O$_2$ was correlated with the thermal expansion and isothermal compressibility of the ILs, both of which are related to the formation of voids in the IL. Therefore, we test the predictions of our model against the density, thermal expansion and isothermal compressibility.

Allegro model for . (a) - Schematic of , with approximate scale of each ion. (b) - Structure of  showing its heterogeneous nature. (c) - Force parity plot between Allegro and DFT for all the elements and components from a single frame at 50 ps through an Allegro simulation. Force MAE of each element is shown in the legend in units of meVÅ$^{-1}$ (d) - Energy comparison between Allegro, and that calculated by DFT using the Allegro structure, as a function of time. The grey values are also Allegro, but were not used to calculate energy/forces in DFT. (e) - Energy parity plot between Allegro and DFT, half of the frames of which are shown in (d). (f) - Comparison of density as a function of temperature for the experiments, classical force field, and the Allegro model (see  for details).

In the SI we provide an outline of the results of changing certain hyperparameters of Allegro, and the affect of different data-sets (for the structures generated directly from OPLS, or these structures propagated further with AIMD or relaxations). Here we will only show results for one model. Specifically a model with $r_{max} = 7$, $l_{max} = 2$, $n_{layers} = 2$ (referred to as $722$), which proved to be reliable, with a scale factor of 100 on the stress (referred to as $11100$, as the scale of the energy and force are 1), using the data set of 180 frames generated from running short AIMD simulations from the uncorrelated OPLS structures (see Methods and ) with an additional 40 frames of volume scans generated from iterative training.

To test this Allegro model, we run several short 100 ps NPT simulations with 530 atoms to collect 20 structure predictions (separated by $\sim$10 ps), and their corresponding energy and forces. In Fig. <ref>(c) we show a parity plot for the forces predicted by Allegro in comparison to the DFT forces computed from these structures. Overall, excellent agreement can be found, with force MAEs in the range of 20-70 meVÅ$^{-1}$. These errors are practically identical to the MAEs obtained validating the potential. Therefore, the Allegro model for  appears to have ab initio accurate forces for this IL, with only a small amount of DFT data. In Fig. <ref>(d) we show how the energy evolves with the simulation from Allegro, and the frames picked out to test the energy and forces on. In Fig. <ref>(d), we also show the corresponding DFT energy calculated from those structures. Clearly, the DFT energy is higher than the Allegro model for all tested structures, by $\sim$2-4 meV/atom, which is several times the error obtained from the data set that it is tested on, albeit still a respectable value for the amount of data. Importantly, the difference in energy does not get significantly worse with time. In Fig. <ref>(e) we further show this with an energy parity plot, where the DFT values are clearly larger than the Allegro energies. It is likely that this difference between Allegro and DFT originates from long-ranged interactions, as Allegro is a strictly local model <cit.>.

In Fig. <ref>(f) we show the density predictions of our Allegro model for  at various temperatures, calculated from the last half of a 1 ns NPT run with 4240 atoms at 1 bar, and we also show the predictions from the classical force field and the experimental values. The experiment finds the density of  to be slightly less than 1.7 gcm$^{-3}$ at 300 K, decreasing with temperature to 1.64 gcm$^{-3}$ at 340 K. The values from the classical force field are extremely close to the experimental values, in part, from the charge rescaling scheme used to match the experimental density (see ). The density prediction from the Allegro model for  is 1.61 gcm$^{-3}$ at 300 K. The difference between the Allegro model and experiments can be attributed to the employed PBE exchange-correlation functional, as found by others <cit.>. Despite this, the thermal expansion from Allegro, -1.03$\times 10^{-3}$gcm$^{-3}$K$^{-1}$, agrees reasonably well with that from experiments, -1.32$\times 10^{-3}$gcm$^{-3}$K$^{-1}$. Moreover, the isothermal compressibility of  is determined to be 8.6$\times 10^{-10}$ Pa$^{-1}$ at 298 K from experiments, and from fitting the Tait equation of state to the results from our Allegro model, we find a value of 6.46$\times 10^{-10}$ Pa$^{-1}$ at 300 K. While the thermal expansion and isothermal compression are not in perfect agreement with experiments, they are within what can be expected from the employed exchange-correlation functional, as others have found for different ILs and molten salts <cit.>, which appears to be the main limitation.

As discussed in the Introduction, previous works of training MLIPs for ILs and other electrolytes has either involved running long AIMD simulations to collect data, or to use classical force fields to sample uncorrelated structures and then perform several loops of iterative training to improve the accuracy of the ML IP. As training and re-training  and Allegro can be expensive, sometimes comparable to the time taken to perform the DFT calculations, it is preferable to not have to do many iterations of training to expand the training data set. Here we investigated if we could use OPLS force fields to sample uncorrelated structures, but at the same time not perform many iterations of collecting data and retraining. We found the most reliable data for starting a model is if short AIMD trajectories are run from these OPLS structures, and then a single round of iterative training performed. We also found that using a larger $r_{max} = 7$ was essential for getting accurate predictions of the density with Allegro. Our model is able to achieve 2-3 meV/atom accuracy in the energy and 20-70 meVÅ$^{-1}$ accuracy in the forces, with $\sim$200 frames. The methodology outlined here be used to first generate a small, high-quality data set (100-200 frames), from short AIMD trajectories starting from the uncorrelated structures generated with OPLS, which could then be expanded in a single (or few) iterative training loop with structures obtained from running this initial potential to obtain a final, robust potential. This could avoid the costly (re)-training procedure incurred for neural network IPs.

The main limitation of our Allegro model for  is that its predictions do not match experiments exactly, but the deviation we obtain is within an error expected from the employed exchange correlation functional <cit.>. We used the cheap PBE+D3 for testing purposes, but using the training protocol developed we can confidently and efficiently use more accurate exchange-correlation functional, which will be shown in future work.

§ ACKNOWLEDGEMENTS

This work used beamline 12-ID-B at the Advanced Photon Source, a U.S. Department of Energy (DOE) Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory under contract no. DE-AC02-06CH11357. This research was supported under a Multidisciplinary University Research Initiative, sponsored by the Department of the Navy, Office of Naval Research, under Grant N00014-20-1-2418. Work at Harvard by Z.A.H.G., J.H.Y., K.B., A.J., A.M., B.K., and N.M. was supported by Robert Bosch Research North America. J.H.Y. acknowledges funding from the Harvard University Center for the Environment. J.A.M. acknowledges support from the Arnold and Mabel Beckman Foundation through a Beckman Young Investigator grant.