Exploring Li-ion Transport Properties of Li 3 TiCl 6 : A Machine Learning Molecular Dynamics Study

,


I. INTRODUCTION
Addressing the increasing energy demands in the face of climate change concerns requires a sustainable zerocarbon energy future.Rechargeable batteries that are capable of converting electrical energy to chemical energy and vice versa, are pivotal for energy storage in this context.While lithium-ion batteries have proven successful for portable devices, all-solid-state batteries (ASSLBs) present a promising solution for the next generation.ASSLBs offer enhanced safety and a longer lifespan compared to traditional lithium-ion electric vehicle batteries and are aligning with goals for achieving zerocarbon emissions.
In ASSBs, the cathode is a crucial component and it is playing a vital role in determining their overall performance.Particularly, energy storage capacity, voltage, cycle life, safety, cost, and environmental impact, are intricately linked to the cathode material.Among the widely used sulfide-and oxide-based cathode materials, latter take a leading position in energy storage manufacturing.Indeed, such commercial attention is also facilitated by the use of transition metals, including single-LiM O 2 , multi-LiNi 1−x−y Mn x Co y O 2 , and polyanionic-TM-based LiM PO 4 (M =3d-block elements of 3 group; x =y=0.1 to 0.33), [1][2][3][4][5][6][7][8][9][10] which essentially improve the performance of ASSLB devices.Since the projected annual production of Li-ion batteries is expected to reach several terawatt hours, demand for Fe, Co, and Ni based cathode materials increases rapidly [11].However, Co and Ni prove to be expensive.Consequently, ongoing efforts include both simulations and experimental studies, aiming to identify more cost-efficient alternatives for cathodes [12].
Among various studies on cathode materials, a recent research has shown that Li 3 TiCl 6 is both cost-effective and outperforms previous benchmarks for ASSLBs [13].However, the experimental investigation [13] of Li 3 TiCl 6 falls short of providing a comprehensive understanding of the underlying physics and chemistry governing the Li-ion transport mechanism, which is a key factor to determine the performance entire battery.Thus, we turn to atomic simulation techniques to study the underlying Li-ion transport mechanism in Li 3 TiCl 6 .
The conventional atomic simulation based on density functional theory (DFT) is well known for predicting structural, electrochemical, and Li-ion transport properties [12,14] with few hundred atoms.Nevertheless, it requires extensive computational resources and is helpless for large-scale demand posed by Liion intercalation-driven electrochemical studies.To facilitate this challenge, various machine learning (ML) methods are utilized in conjunction with molecular dynamics simulations [15,16].
Therefore, our approach in this work involves the inte- gration of AIMD, machine learning methods based on deep neural networks, and classical molecular dynamics.This integrated approach is collectively referred to as deep learning molecular dynamics (DLMD) simulations [22,23] to investigate the structural and Li-ion transport properties of the Li 3 TiCl 6 cathode.With this, we organize the manuscript as follows: Section II describes detailed simulation techniques, as illustrated in FIG. 1. Section III is composed of results and discussions covering structural parameters, RDF, MSD, self and correlated motion of Li-ion displacement, diffusion coefficient, Li-ion transport mechanism, ionic conductivity, and activation energy of Li 3 TiCl 6 .The calculated values of ionic conductivity and activation energy are in good agreement with those of experimental values [13].Finally, we conclude our results and discussions in Section IV.

II. SIMULATION METHODOLOGIES
A. AIMD simulation DFT-based calculations were employed to optimize the crystal structure of Li 3 TiCl 6 using the Vienna Ab−initio Simulation Package (VASP) [24,25].The initial structure of Li 3 TiCl 6 in monoclinic cell with C2/m space group for VASP was adopted from experimental results [13].The projector augmented wave (PAW) formalism described the valence electrons of Li, Cl, and Ti atoms using plane wave-based wave functions was employed [26].The structure optimization, involving the minimization of ground state energy, utilized the generalized gradient approximation of Perdew and Wang method with different Uparameters [27,28].A set of U − J values, 0 and 4eV were chosen to taking into account of strong correlation effect of Ti-3d electrons [29,30].A kinetic cutoff energy of 450 eV was set to enhance calculation accuracy.The Brillouin zone of the supercell with 160 and 320 atoms was sampled with 4×2×4 and 4×2×2 k-meshes for ionic optimization, respectively.Ionic and electronic optimizations were alternately performed until the forces on each ion reached less than ±10 meV/ Å.
The designed simulation cells were optimized and subjected to AIMD simulations.The time step for the AIMD simulations was set to 1 fs, and the simulations were continued for a duration considered separately within the canonical NVT ensemble to generate the dataset for the DLP model.The DLP was developed using the deep neural network method implemented in DeePMD-kit (v2.2.7) [31].The deep neural network algorithm in DeePMD is designed using the TensorFlow Python library [32].The deep potential-smooth edition (DeepPot-SE), an end-to-end machine learning-based potential energy surface (PES) model, was employed with a cutoff radius of 7 Å to include more neighbour atoms in the featurization process.Indeed, it efficiently represents the PES of a wide variety of systems with the accuracy of ab − initio models [33].
In the process of constructing DLP, local coordinate matrix, R and local atomic environment matrix, {R ij } N i=1 are represented as shown in Eq. ( 1) and ( 2).
where r i = (x i , y i , z i ) which contains 3 Cartesian coordinates of atom and N is total number of atoms.And R can be transformed into local environment matrices as where j and N Rc (i) are indexes and number of neighbors of i th atom within the cut-off radius r c = 7.0 Å, and r ji ≡ r j − r i is defined as the relative coordinate.
An embedding neural network with three layers, each containing 32, 64, and 128 neurons, was used to convert the local atomic coordinates (R ∈ R N ×3 ) into atomic descriptors Di = Di {R ij } N i=1 that preserves the structural symmetries of the system [31,34].Descriptor values were input into a fitting neural network with three layers, each comprising 256 neurons, that maps descriptors to atomic energies E i [31,34], and thus the total energy and force of each atom in the system are calculated by Eq. (3) and Eq. ( 4).
The hyperbolic tangent (tanh) activation function was applied in the neural network, introducing non-linearity to effectively train the intricate atomic descriptor data D i .The training process consisted of 10 4 steps, utilizing mean squared error for both energies and forces as the loss function at each training step.The optimization was facilitated by the Adam optimizer, initiating with a learning rate of 10 −3 and concluding at 10 −8 , with a decay parameter set to 5000.With these specified parameters, a dataset comprising 1.3 million training samples and 10000 test samples from diverse trajectories were employed to construct the DLP model.The calculated loss function parameters such as mean absolute error and root mean square error during validation of developed DLP model is tabulated in Table-1.Accuracy of the predicted data of force and energy is more than 99 % and data points of predicted versus trained are shown in FIG.1b to FIG. 1e.

C. Deep learning molecular dynamic simulation
Similarly to AIMD conditions, The simulations of DLMD was conducted in the NVT ensemble using the LAMMPS simulation package [35] coupled with the DeepMD plugin [34].The simulation cell was designed to accommodate 20,000 atoms.Simulations were carried out at six different temperatures-298, 313, 328, 333, 358, and 373 K-controlled by the Nose-Hoover thermostat.
Prior to NVT simulation, the energy minimization was conducted using conjugate gradient algorithm.The temperature damping parameter was set to 100 fs, and a uniform integration time-step of 1 fs was employed for all simulations, extending over a total duration of 5.5 ns.Computing time for AIMD and DLMD simulation for 10 ps per atom on one computing core was calculated and plotted in FIG.1g and FIG.1g, respectively and it reveals that DLMD is 3730 times faster than AIMD simulation.The unit-cell of AIMD simulations were visualized using VESTA software [36] and the trajectories of AIMD and DLMD simulations were visualized using OVITO software.[37]

A. Radial distribution
The monoclinic simulation cells, comprising 160 and 320 atoms, were employed to create three distinct Ti occupancies at the 2a and 4g sites, aiming to reproduce the experimentally reported structure of Li 3 TiCl 6 [13] (FIG.2).The three different occupancies, along with the corresponding experimental Ti occupancy, are summarized in Table-2.
TABLE II.Calculated lattice parameters of monoclinic Li3TiCl6 with different Ti occupancy at 4h and 2a sites.0.875Ti@2a 0.063Ti@4g 0.851Ti@2a 0.075Ti@4g 0.844Ti@2a 0.078Ti@4g 0.754Ti@2a 0.123Ti@4g 0.750Ti@2a 0.125Ti@4g Lattice constant, a( Å) Subsequently, the designed cells underwent optimization, and their structural parameters closely matched experimental values [13] for DFT+U calculation with U=4.0 eV.The slight variations, when compared to their experimental values, may be attributed to the small change in the Ti occupancy within the designed simulation cells as well as the overestimation of GGA functional that used in the DFT simulation.Based on the structural optimization results, we selected the Li 3 TiCl 6 structure with Ti-sites at 0.754 on 2a and 0.123 on 4g, which represents experimentally annealed structure at 300 o C and exhibited a maximum Li-ion conductivity of 1.04 mS/cm at room temperature [13] .
Atomic trajectories produced by AIMD and DLMD at 298 K were compared in the framework of radial distribution function (RDF), calculated by the equation given below: where, N , r ij , and ∆r represent the total number of atoms in the radius r, the distance between atoms i and j, and the width of each bin, respectively.
Generally, RDFs between pairs of atoms in the same material, estimated using AIMD and DLMD, may provide valuable confirmation of reliability, accuracy, and Comparison of radial distribution function of Li3TiCl6 analyzed based on AIMD with U=4.0 eV and DLMD trajectories at 298K.consistency in capturing the structural features of the system under investigation.Thus, results of the RDF analyzed using two theoretical approximations for Cl-Ti, Cl-Li, Ti-Li, Li-Li, and Cl-Cl pairs of atoms are presented in FIG. 3. As can be clearly seen, all the peak positions, shapes, and heights of RDF obtained by means of DLMD exhibit good agreement with AIMD results, even at larger distance, confirming the ability of DLMD to predict other static properties with high-level accuracy.
Considering the RDF results more carefully, one can observe multiple peaks along the large separation distances which indicates a relatively high degree of atomic structuring in Li 3 TiCl 6 .For Cl-Ti and Cl-Li pairs, the first peak position is nearly the same at a distance of 2.53 Å, while other pair combinations are situated at larger distances of 3.5 Å.At the same time, the height of the Li-Cl peak is greater than that of Ti-Cl, which reveals that a relatively higher probability of Li-Cl interactions compared to Ti-Cl.This is evident when comparing the peak heights of Li-Li, Ti-Li, and Cl-Cl pairs.

B. Diffusion coefficients
The diffusion of Li-ions was determined from DLMD trajectories, where the positions of Li atoms, denoted as {r i (t)} N i=1 , are tracked as a function of time t for all N Li atoms.This involves calculating the mean square displacement of Li atoms, MSD Total [38,39] as follows.
We compute uncorrelated motion of MSD as follows.
From Eq. ( 6) and Eq. ( 7), distinct part of MSD is calculated by When MSD distinct < 0 describes negative correlation between Li-ions and, hence, negatively affect the Liion transport.Based on Equations (6-8), we calculated MSD self and MSD distinct and collect the results in FIG. 4. It is seen that MSD distinct values for Li 3 TiCl 6 at different temperatures show negative correlation (or anticorrelation) between Li atoms during their dynamics.To quantify this observation the diffusion coefficient was calculated on the basis of MSD using where d is dimensionless quantity which is equal to 3 for three dimensional transport.Since, Li 3 TiCl 6 has anticorrelated Li-ion dynamics, calculated D self is highly overestimated about 14 times of D T otal , which is 1.6 × 10 −14 cm 2 /s at 298K with DLP generated from DFT+U = 4.00 eV.
One way to understand this phenomenon is by analyzing Haven's ratio (H R ) [40] which is defined as the ratio of the self-diffusion coefficient of Li-ions (D self ) to the total diffusion coefficient (D Total ) using Eq. 9.
Here, D self represents the individual motion of ions which is independent of overall charge transport.On the other hand, D T otal is the diffusion coefficient calculated from the material's conductivity by considering the collective motion of ions that contributes to overall charge transport.In the case of uncorrelated Li-ion dynamics, D Total ≡ D self , and H R becomes 1.If H R is less than 1, D distinct should have been positively correlated with Li-ion dynamics, typically observed in liquid and glass electrolytes.[39,40] In our case, since H R is greater than one, we propose that Li-ion dynamics may involve an interstitial and/or inter-layer transport mechanism.

C. Transport mechanism
Furthermore, we analyzed the Li-ion path to confirm the inter-and intra-layer diffusion of the ions.Among all the Li atoms in the trajectory simulated at 298 K, we arbitrarily selected two Li atoms to track the ionic motion over time (in FIG. 5).Figures 5a and 5b clearly illustrate the Li atom (gray-colored dots) moving within the Li layer through interstitial hopping.Simultaneously, interlayer motion between Li-and Ti-layers is also observed through interstitial migration, as shown in Fig. 5d and  5e.The movements of ionic movement corresponding to intra-and inter-layer migration are shown in Fig. 5c and 5f.More interestingly, Fig. 5f demonstrates that inter-layer Li-ion migration reaches the third cation layer at 4.7 ns through interstitial sites, proving multiple-site hopping of Li-ion.Since the Li 3 TiCl 6 crystal structure has partially occupied Ti-2a as well as Ti-4g, the structure possesses inherent voids that can also possibly act as hopping sites for Li-ion migration.Therefore, our study confirms that Li 3 TiCl 6 has an inter-and intra-layer interstitial hopping-based Li-ion transport mechanism.

D. Ionic conductivity
The ionic conductivity, denoted as σ, can be determined from the self-diffusion coefficient using the Nernst-Einstein equation [41]: Here, n represents the ion density of Li, e is the elementary electron charge, Z denotes the valence of Li, and k B is the Boltzmann constant.The calculated ionic conductivity values from diffusion coefficients align well with experimental values [13], as demonstrated in FIG. 6.Additionally, when considering the temperature dependence of ionic conductivity in solid-state electrolytes, hightemperature ionic conductivities obtained from DLMD simulations can be leveraged to estimate the activation barrier of the electrolytes at lower temperatures using the Arrhenius relationship: σ = σ 0 e −Ea/(kBT ) .The calculated activation energy is 0.29 eV, closely matching the experimental value of 0.32 eV.In addition, Li-ion migra-

IV. CONCLUSION
We performed large-scale and MLFF-driven molecular dynamics simulations to investigate the Li-ion transport mechanism in cation-disordered Li 3 TiCl 6 cathode at six different temperatures, ranging from 298K to 373K.Deep neural network method along with data generated by AIMD simulation were used to build a high-fidelity MLFF.Predicting accuracy of atomic forces, energy, and structure by our trained MLFF was confirmed with set of new AIMD data and corresponding RDF.The calculated self and distinct part of Li-ion MSD reveal that the Li-ions are involved in anti-correlation motion that was rarely reported for solid-state materials.
In the same way, analysis of trajectory from DLMD infers that the Li-ion transportation occurs through interstitial hopping which was confirmed by intra-and interlayer Li-ion movement as a function of simulation time.The temperature dependent ionic conductivity and, thus, activation barrier values for Li 3 TiCl 6 demonstrate a decreasing trend with temperature, aligning with typical behavior of ionic conductors.Moreover, activation energy of 0.29 eV, which is in close agreement with experimental result, matches well with inter-layer ionic diffusion barrier calculated by DFT along [110] crystallographic direction.Overall, the combination of machinelearning methods and AIMD simulations elucidates the complex Li-ion electrochemical properties of the Li 3 TiCl 6 cathode by significantly reducing computational time.Hence, our work strongly suggests that the MLFF using deep neural networks could be promising for studying large-scale complex materials.

FIG. 1 .
FIG. 1.(a) Schematic diagram of the protocol of DLMD development (see Section 2 for the details).(b-d) DLMD predicted forces corresponding to AIMD forces for Cl (b) , Li (c) and Ti (d) elements, as well as the respective forces along x, y, and z directions.(e) DLMD predicted energies corresponding to AIMD energies.(e-f) Simulation time (Ct) of AIMD(e) and DLMD(f) approaches, respectively.

FIG. 5 .
FIG. 5. Illustration of intra-layer Li-ion movement (gray balls) along Li3TiCl6 (a) [010], and (b) [110] view-point directions.(c) Time-evolution of Li atom translations corresponding to (a) and (b) as a function of time.Illustration of inter-layer Liion movement (gray balls) along Li3TiCl6 (d) [010], and (e) [110] view-point directions.(f) Translation distance of Li atoms corresponding to (d) and (e) as a function of time.

TABLE I .
The calculated loss function parameter expressed by mean absolute error (MAE) and root mean square error (RMSE) of deep learning potential model