Graph Dynamical Networks for Unsupervised Learning of Atomic Scale Dynamics in Materials

TL;DR

This paper introduces graph dynamical networks, an unsupervised learning method that captures atomic scale dynamics from molecular simulations, aiding materials design by analyzing complex local environments.

Contribution

The work presents a novel unsupervised graph neural network approach for learning atomic dynamics across various phases and environments from molecular dynamics data.

Findings

Successfully learned dynamical information for multi-component amorphous materials

Applicable to diverse phases and environments in molecular simulations

Provides an automated tool for analyzing atomic scale dynamics

Abstract

Understanding the dynamical processes that govern the performance of functional materials is essential for the design of next generation materials to tackle global energy and environmental challenges. Many of these processes involve the dynamics of individual atoms or small molecules in condensed phases, e.g. lithium ions in electrolytes, water molecules in membranes, molten atoms at interfaces, etc., which are difficult to understand due to the complexity of local environments. In this work, we develop graph dynamical networks, an unsupervised learning approach for understanding atomic scale dynamics in arbitrary phases and environments from molecular dynamics simulations. We show that important dynamical information can be learned for various multi-component amorphous material systems, which is difficult to obtain otherwise. With the large amounts of molecular dynamics data generated everyday in nearly every aspect of materials design, this approach provides a broadly useful, automated tool to understand atomic scale dynamics in material systems.

Tables

Table 1. Supplementary Table 1: The charge carried by each state in PEO/LiTFSI.

State	0	1	2	3
Charge	+0.040	+0.262	-0.637	+0.889

Table 2. Supplementary Table 2: The diffusion coefficient of each state transition in PEO/LiTFSI at τ = 0.8   ns 𝜏 times 0.8 ns \tau=$0.8\text{\,}\mathrm{n}\mathrm{s}$ . (Unit: 10 − 7 ​   cm 2 / s superscript 10 7 times absent superscript cm 2 s 10^{-7}$\text{\,}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ )

Transition	$j=0$	$j=1$	$j=2$	$j=3$
$i=0$	$3.9\pm 0.2$	$4.6\pm 0.1$	$4.6\pm 0.2$	$4.0\pm 0.2$
$i=1$	$4.5\pm 0.2$	$4.3\pm 0.1$	$4.6\pm 0.1$	$4.4\pm 0.2$
$i=2$	$4.4\pm 0.1$	$4.6\pm 0.1$	$3.5\pm 0.0$	$4.8\pm 0.7$
$i=3$	$3.8\pm 0.1$	$3.9\pm 0.1$	$5.5\pm 0.6$	$2.8\pm 0.2$

Equations

$${\mathbf{x}}_{t+\tau }=\mathbf{F}({\mathbf{x}}_t),$$

$$\mathbf{\chi }({\mathbf{x}}_{t+\tau })\approx {\mathbf{K}}^T\mathbf{\chi }({\mathbf{x}}_t).$$

$${\mathbf{v}}_i^{\prime }=\text{Conv}({\mathbf{v}}_i,{\mathbf{v}}_j,{\mathbf{u}}_{(i,j)}),(i,j)\in \mathcal{G}.$$

$${\mathbf{v}}_{[l]}=\mathrm{Pool}({\mathbf{v}}_{[l,0]},{\mathbf{v}}_{[l,1]},\dots ,{\mathbf{v}}_{[l,k]}).$$

$$\text{Loss}=-\mathrm{VAMP}({\tilde{\mathbf{v}}}_{[l],t},{\tilde{\mathbf{v}}}_{[l],t+\tau }),t\in [0,T-\tau ],l\in [0,n].$$

$$\mathbf{K}(n\tau )={\mathbf{K}}^n(\tau ),n=1,2,\dots .$$

$${\mathbf{u}}_{(i,j)}[t]=\exp (-(d_{(i,j)}-{\mu }_t{)}^2/{\sigma }^2),$$

$${\alpha }_{ij}=\frac{\exp ({\mathbf{z}}_{(i,j)}^{(t-1)}{\mathbf{W}}_{\mathrm{a}}^{(t-1)}+b_{\mathrm{a}}^{(t-1)})}{{\sum }_j\exp ({\mathbf{z}}_{(i,j)}^{(t-1)}{\mathbf{W}}_{\mathrm{a}}^{(t-1)}+b_{\mathrm{a}}^{(t-1)})},$$

$${\mathbf{v}}_i^{(t)}={\mathbf{v}}_i^{(t-1)}+\sum _j{\alpha }_{ij}\cdot g({\mathbf{z}}_{(i,j)}^{(t-1)}{\mathbf{W}}_{\mathrm{n}}^{(t-1)}+{\mathbf{b}}_{\mathrm{n}}^{(t-1)}),$$

$${\mathbf{v}}_{[l]}=\frac{1}{k}\sum _m{\mathbf{v}}_{[l,m]}.$$

$$t_i(\tau )=-\frac{\tau }{\ln |{\lambda }_i(\tau )|},$$

$$g_i(r_{\mathrm{A}})=\frac{1}{{\rho }_i}\frac{\mathrm{d}[n(r_{\mathrm{A}})\cdot p_i]}{4\pi r_{\mathrm{A}}^2\mathrm{d}{r}_{\mathrm{A}}},$$

$$\mathbb{E}[d^2(t)|i\to j]=\frac{{\sum }_{t^{\prime }}{d}^2(t^{\prime },t^{\prime }+t)p_i(t^{\prime })p_j(t^{\prime }+t)}{{\sum }_{t^{\prime }}{p}_i(t^{\prime })p_j(t^{\prime }+t)},$$

$$D_{ij}(\tau)=\frac{1}{6}\frac{\mathop{}\!\mathrm{d}\mathbb{E}[d^{2}(t)|i\rightarrow j]}{\mathop{}\!\mathrm{d}t}\Bigr{\rvert}_{t=\tau},$$

$${\sigma }_{ij}=\frac{e^2{N}_{\mathrm{Li}}}{V{k}_{\mathrm{B}}T}{\pi }_i{z}_{ij}{K}_{ij}(\tau )D_{ij}(\tau ),$$

$$\frac{{\sigma }_{ij}}{{\sum }_{i,j}{\sigma }_{ij}}.$$

$${\mathbf{K}}_{\mathrm{local}}=\left[\begin{array}{cc}{p}_o& 1-p_o\\ 1-p_t& p_t\end{array}\right],$$

$${\mathbf{K}}_{\mathrm{global}}=\text{setcounter}MaxMatrixCols12\left[\begin{array}{cccccccccccc}{p}_o& 0& 0& 0& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}\\ 0& p_o& 0& 0& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}\\ 0& 0& p_o& 0& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}\\ 0& 0& 0& p_o& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}& \frac{1-p_o}{8}\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& p_t& 0& 0& 0& 0& 0& 0& 0\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& 0& p_t& 0& 0& 0& 0& 0& 0\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& 0& 0& p_t& 0& 0& 0& 0& 0\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& 0& 0& 0& p_t& 0& 0& 0& 0\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& 0& 0& 0& 0& p_t& 0& 0& 0\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& 0& 0& 0& 0& 0& p_t& 0& 0\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& 0& 0& 0& 0& 0& 0& p_t& 0\\ \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& \frac{1-p_t}{4}& 0& 0& 0& 0& 0& 0& 0& p_t\end{array}\right].$$

$${\mathrm{Charge}}_i=1-{\int }_0^{r\ast }{g}_i(r_N)\cdot 4\pi r^2\mathrm{d}r+{\int }_0^{r\ast }{g}_i(r_{Li})\cdot 4\pi r^2\mathrm{d}r,$$

Full text

Introduction

Understanding the atomic scale dynamics in condensed phase is essential for the design of functional materials to tackle the global energy and environmental challenges. Etacheri et al. (2011); Imbrogno and Belfort (2016); Peighambardoust et al. (2010) The performance of many materials, like electrolytes and membranes, depend on the dynamics of individual atoms or small molecules in complex local environments. Despite the rapid advances in experimental techniques Zheng et al. (2016); Yu et al. (2016); Perakis et al. (2016), molecular dynamics (MD) simulations remain one of the few tools for probing these dynamical processes with both atomic scale time and spatial resolutions. However, due to the large amounts of data generated in each MD simulation, it is often challenging to extract statistically relevant dynamics for each atom especially in multi-component, amorphous material systems. At present, atomic scale dynamics are usually learned by designing system-specific description of coordination environments or computing the average behavior of atoms. Wang et al. (2015); Borodin and Smith (2006); Miller III et al. (2017); Getman et al. (2011) A general approach for understanding the dynamics in different types of condensed phases, including solid, liquid, and amorphous, is still lacking.

The advances in applying deep learning to scientific research open new opportunities for utilizing the full trajectory data from MD simulations in an automated fashion. Ideally, we want to trace every atom or small molecule of interest in the MD trajectories, and summarize their dynamics into a linear, low dimensional model that describes how their local environments evolve over time. Recent studies show that combining Koopman analysis and deep neural networks provides a powerful tool to understand complex biological processes and fluid dynamics from data. Li et al. (2017); Lusch et al. (2018); Mardt et al. (2018) In particular, VAMPnets Mardt et al. (2018) develop a variational approach for Markov processes to learn an optimal latent space representation that encodes the long-time dynamics, which enables the end-to-end learning of a linear dynamical model directly from MD data. However, in order to learn the atomic dynamics in complex, multi-component material systems, sharing knowledge learned for similar local chemical environments is essential to reduce the amount of data needed. Recent development of graph convolutional neural networks (GCN) has led to a series of new representations of molecules Duvenaud et al. (2015); Kearnes et al. (2016); Gilmer et al. (2017); Schütt et al. (2017) and materials Xie and Grossman (2018a); Schütt et al. (2018) that are invariant to permutation and rotation operations. These representations provide a general approach to encode the chemical structures in neural networks which shares parameters between different local environments, and they have been used for predicting properties of molecules and materials Duvenaud et al. (2015); Kearnes et al. (2016); Gilmer et al. (2017); Schütt et al. (2017); Xie and Grossman (2018a); Schütt et al. (2018), generating force fields Schütt et al. (2018); Zhang et al. (2018), and visualizing structural similarities Zhou et al. (2018); Xie and Grossman (2018b).

In this work, we develop a deep learning architecture, Graph Dynamical Networks (GDyNets), that combines Koopman analysis and graph convolutional neural networks to learn the dynamics of individual atoms in material systems. The graph convolutional neural networks allow for the sharing of knowledge learned for similar local environments across the system, and the variational loss developed in the VAMPnets is employed to learn a linear model for atomic dynamics. Our method focuses on the modeling of local atomic dynamics instead of global dynamics. This significantly improves the sampling of the atomic dynamical processes because a typical material system includes a large number of atoms or small molecules moving in structurally similar but distinct local environments. We demonstrate this distinction using a toy system that shows global dynamics can be exponentially more complex than local dynamics. Then, we apply this method to two realistic material systems – silicon dynamics at solid-liquid interfaces and lithium ion transport in amorphous polymer electrolytes – to demonstrate the new dynamical information one can extract for such complex materials and environments. Given the enormous amount of MD data generated in materials research, we believe the broad applicability of this method could help uncover important new physical insights from atomic scale dynamics that may have otherwise been overlooked.

Results

Koopman analysis of atomic scale dynamics. In materials design, the dynamics of target atoms, like the lithium ion in electrolytes and the water molecule in membranes, provide key information to material performance. We describe the dynamics of the target atoms and their surrounding atoms as a discrete process in MD simulations,

[TABLE]

where $\bm{x}_{t}$ and $\bm{x}_{t+\tau}$ denote the local configuration of the target atoms and their surrounding atoms at time steps $t$ and $t+\tau$, respectively. Note that Eq. (1) implies that the dynamics of $\bm{x}$ is Markovian, i.e. $\bm{x}_{t+\tau}$ only depends on $\bm{x}_{t}$ not the configurations before it. This is exact when $\bm{x}$ includes all atoms in the system, but an approximation if only neighbor atoms are included. We also assume that each set of target atoms follow the same dynamics $\bm{F}$. These are valid assumptions since 1) most interactions in materials are short-range, 2) most materials are either periodic or have similar local structures, and we could test them by validating the dynamical models using new MD data which we will discuss later.

The Koopman theory Koopman (1931) states that there exists a function $\bm{\chi}(\bm{x})$ that maps the local configuration of target atoms $\bm{x}$ into a lower dimensional feature space, such that the non-linear dynamics $\bm{F}$ can be approximated by a linear transition matrix $\bm{K}$,

[TABLE]

The approximation becomes exact when the feature space has infinite dimensions. However, for most dynamics in material systems, it is possible to approximate it with a low dimensional feature space with a sufficiently large $\tau$ due to the existence of characteristic slow processes. The goal is to identify such slow processes by finding the feature map function $\bm{\chi}(\bm{x})$.

Learning feature map function with graph dynamical networks. In this work, we use graph convolutional neural networks (GCN) to learn the feature map function $\bm{\chi}(\bm{x})$. GCN provides a general framework to encode the structure of materials that is invariant to permutation, rotation, and reflection Xie and Grossman (2018a); Schütt et al. (2018). As shown in Fig. 1, for each time step in the MD trajectory, a graph $\mathcal{G}$ is constructed based on its current configuration with each node $\bm{v}_{i}$ representing an atom and each edge $\bm{u}_{i,j}$ representing a bond connecting nearby atoms. We connect $M$ nearest neighbors considering periodic boundary condition while constructing the graph, and a gated architecture Xie and Grossman (2018a) is used in GCN to reweigh the strength of each connection (see Supplementary Note 1 for details). Note that the graphs are constructed separately for each step, so the topology of each graph may be different. Also, the 3-dimensional information is preserved in the graphs since the bond length is encoded in $\bm{u}_{i,j}$. Then, each graph is input to the same GCN to learn an embedding for each atom through graph convolution (or neural message passing Gilmer et al. (2017)) that incorporates the information of its surrounding environments.

[TABLE]

After $K$ convolution operations, information from the $K$th neighbors will be propagated to each atom, resulting in an embedding $\bm{v}_{i}^{(K)}$ that encodes its local environment.

To learn a feature map function for the target atoms whose dynamics we want to model, we focus on the embeddings learned for these atoms. Assume that there are $n$ sets of target atoms each made up with $k$ atoms in the material system. For instance, in a system of 10 water molecules, $n=10$ and $k=3$. We use the label $\bm{v}_{[l,m]}$ to denote the $m$th atom in the $l$th set of target atoms. With a pooling function Xie and Grossman (2018a), we can get an overall embedding $\bm{v}_{[l]}$ for each set of target atoms to represent its local configuration,

[TABLE]

Finally, we build a shared output layer with a Softmax activation function to map the embeddings $\bm{v}_{[l]}$ to a feature space $\tilde{\bm{v}}_{[l]}$ with a pre-determined dimension. This is the feature space described in Eq. (2), and we can select an appropriate dimension to capture the important dynamics in the material system. The Softmax function used here allows us to interpret the feature space as a probability over several states Mardt et al. (2018). Below, we will use the term “number of states” and “dimension of feature space” interchangeably.

To minimize the errors of the approximation in Eq. (2), we compute the loss of the system using a VAMP-2 score Wu and Noé (2017); Mardt et al. (2018) that measures the consistency between the feature vectors learned at timesteps $t$ and $t+\tau$,

[TABLE]

This means that a single VAMP-2 score is computed over the whole trajectory and all sets of target atoms. The entire network is trained by minimizing the VAMP loss, i.e. maximizing the VAMP-2 score, with the trajectories from the MD simulations.

Hyperparameter optimization and model validation. There are several hyperparameters in the GDyNets that need to be optimized, including the architecture of GCN, the dimension of the feature space, and lag time $\tau$. We divide the MD trajectory into training, validation, and testing sets. The models are trained with trajectories from the training set, and a VAMP-2 score is computed with trajectories from the validation set. The GCN architecture is optimized according to the VAMP-2 score similar to ref. Xie and Grossman (2018a).

The accuracy of Eq. (2) can be evaluated with a Chapman-Kolmogorov (CK) equation,

[TABLE]

This equation holds if the dynamic model learned is Markovian, and it can predict the long-term dynamics of the system. In general, increasing the dimension of feature space makes the dynamic model more accurate, but it may result in overfitting when the dimension is very large. Since a higher feature space dimension and a larger $\tau$ make the model harder to understand and contain less dynamical details, we select the smallest feature space dimension and $\tau$ that fulfills the CK equation within statistical uncertainty. Therefore, the resulting model is interpretable and contains more dynamical details. More about the effects of feature space dimension and $\tau$ can be found in refs. Mardt et al. (2018); Wu and Noé (2017).

Local and global dynamics in the toy system. To demonstrate the advantage of learning local dynamics in material systems, we compare the dynamics learned by the GDyNet with VAMP loss and a standard VAMPnet with fully connected neural networks that learns global dynamics for a simple model system using the same input data. As shown in Fig. 2(a), we generated a $200\text{\,}\mathrm{n}\mathrm{s}$ MD trajectory of lithium atom moving in a face-centered cubic (FCC) lattice of sulfur atoms at a constant temperature, which describes an important lithium ion transport mechanism in solid-state electrolytes Wang et al. (2015). There are two different sites for the lithium atom to occupy in a FCC lattice, tetrahedral sites and octahedral sites, and the hopping between the two sites should be the only dynamics in this system. As shown in Fig. 2(b-d), after training and validation with the first $100\text{\,}\mathrm{n}\mathrm{s}$ trajectory, the GDyNet correctly identified the transition between the two sites with a relaxation timescale of $42.3\text{\,}\mathrm{p}\mathrm{s}$ while testing on the second $100\text{\,}\mathrm{n}\mathrm{s}$ trajectory, and it performs well in the CK test. In contrast, the standard VAMPnet, which inputs the same data as the GDyNet, learns a global transition with a much longer relaxation timescale at $236\text{\,}\mathrm{p}\mathrm{s}$, and it performs much worse in the CK test. This is because the model views the four octahedral sites as different sites due to their different spatial location. As a result, the transition between these identical sites are learned as the slowest global dynamics.

It is theoretically possible to identify the faster local dynamics from a global dynamical model when we increase the dimension of feature space (Supplementary Fig. 1. However, when the size of system increases, the number of slower global transitions will increase exponentially, making it practically impossible to discover important atomic scale dynamics within a reasonable simulation time. In addition, it is possible in this simple system to design a symmetrically invariant coordinate to include the equivalence of the octahedral and tetrahedral sites. But in a more complicated multi-component or amorphous material system, it is difficult to design such coordinates that take into account the complex atomic local environments. Finally, it is also possible to reconstruct global dynamics from the local dynamics. Since we know how the 4 octahedral and 8 tetrahedral sites are connected in a FCC lattice, we can construct the 12 dimensional global transition matrix from the 2 dimensional local transition matrix (see Supplementary Note 2 for details). We obtain the slowest global relaxation timescale to be $531\text{\,}\mathrm{p}\mathrm{s}$, which is close to the observed slowest timescale of $528\text{\,}\mathrm{p}\mathrm{s}$ from the global dynamical model in Supplementary Fig. 1. Note that the timescale from the two-state global model in Fig. 2 is less accurate since it fails to learn the correct transition. In sum, the built-in invariances in GCN provides a good tool to reduce the complexity of learning atomic dynamics in material systems.

Silicon dynamics in solid-liquid interface. To evaluate performance of the GDyNets with VAMP loss for a more complicated system, we study the dynamics of silicon atoms at a binary solid-liquid interface. Understanding the dynamics at interfaces is notoriously difficult due to the complex local structures formed during phase transitions. Sastry and Angell (2003); Angell (2008) As shown in Fig. 3(a), an equilibrium system made of two crystalline Si {110} surfaces and a liquid Si-Au solution is constructed at the eutectic point ($629\text{\,}\mathrm{K}$, 23.4% Si Ryu and Cai (2010)) and simulated for $25\text{\,}\mathrm{n}\mathrm{s}$ using MD. We train and validate a four-state model using the first $12.5\text{\,}\mathrm{n}\mathrm{s}$ trajectory, and use it to identify the dynamics of Si atoms in the last $12.5\text{\,}\mathrm{n}\mathrm{s}$ trajectory. Note that we only use the Si atoms in the liquid phase and the first two layers of {110} surfaces as the target atoms (Fig. 3(b)). This is because the Koopman models are optimized for finding the slowest transition in the system, and including additional solid Si atoms will result in a model that learns the slower Si hopping in the solid phase which is not our focus.

In Fig. 3(b, c), the model identified four states that are crucial for the Si dynamics at the solid-liquid interface – liquid Si at the interface (state 0), solid Si (state 1), solid Si at the interface (state 2), and liquid Si (state 3). These states provide a more detailed description of the solid-liquid interface structure than conventional methods. In Supplementary Fig. 2, we compare our results with the distribution of the $q_{3}$ order parameter of the Si atoms in the system, which measures how much a site deviates from a diamond-like structure and is often used for studying Si interfaces Wang et al. (2017). We learn from the comparison that 1) our method successfully identifies the bulk liquid and solid states, and learns additional interface states that cannot be obtained from $q_{3}$; 2) the states learned by our method are more robust due to access to dynamical information, while $q_{3}$ can be affected by the accidental order structures in the liquid phase; 3) $q_{3}$ is system specific and only works for diamond-like structures, but the GDyNets can potentially be applied to any material given the MD data.

In addition, important dynamical processes at the solid-liquid interface can be learned with the model. Remarkably, the model identified the relaxation process of the solid-liquid transition with a timescale of $538\text{\,}\mathrm{n}\mathrm{s}$ (Fig. 3(d, e)), which is one order of magnitude longer than the simulation time of $12.5\text{\,}\mathrm{n}\mathrm{s}$. This is because the large number of Si atoms in the material system provide an ensemble of independent trajectories that enable the identification of rare events Pande et al. (2010); Chodera and Noé (2014); Husic and Pande (2018). The other two relaxation processes corresponds to the transitions of solid Si atoms at the interface ($73.2\text{\,}\mathrm{n}\mathrm{s}$) and liquid Si atoms at interface ($2.26\text{\,}\mathrm{n}\mathrm{s}$), respectively. These processes are difficult to obtain with conventional methods due to the complex structures at solid-liquid interfaces, and the results are consistent with our understanding that the former solid relaxation is significantly slower than the latter liquid relaxation. Finally, the model performs excellently in the CK test on predicting the long-term dynamics.

Lithium ion dynamics in polymer electrolytes. Finally, we apply GDyNets with VAMP loss to study the dynamics of lithium ions (Li-ions) in solid polymer electrolytes (SPEs), an amorphous material system composed of multiple chemical species. SPEs are candidates for next generation battery technology due to their safety, stability, and low manufacturing cost, but they suffer from low Li-ion conductivity compared with liquid electrolytes. Meyer (1998); Hallinan Jr and Balsara (2013) Understanding the key dynamics that affect the transport of Li-ions is important to the improvement of Li-ion conductivity in SPEs.

We focus on the state-of-the-art Hallinan Jr and Balsara (2013) SPE system – a mixture of poly(ethylene oxide) (PEO) and lithium bis-trifluoromethyl sulfonimide (LiTFSI) with $\mathrm{Li/EO}=0.05$ as shown in Fig. 4(a). Five independent $80\text{\,}\mathrm{n}\mathrm{s}$ trajectories are generated to model the Li-ion transport at $363\text{\,}\mathrm{K}$. We train a four-state GDyNet with one of the trajectories, and use the model to identify the dynamics of Li-ions in the remaining four trajectories. The model identified four different solvation environments, i.e. states, for the Li-ions in the SPE. In Fig. 4(b), the state 0 Li-ion has a population of $50.6\pm 0.8\%$, and it is coordinated by a PEO chain on one side and a TFSI anion on the other side. The state 1 has a similar structure as state 0 with a population of $27.3\pm 0.4\%$, but the Li-ion is coordinated by a hydroxyl group on the PEO side rather than an oxygen. In state 2, the Li-ion is completely coordinated by TFSI anion ions, which has a population of $15.1\pm 0.4\%$. And the state 3 Li-ion is coordinated by PEO chains with a population of $7.0\pm 0.9\%$. Note that the structures in Fig. 4(b) only show a representative configuration for each state. We compute the element-wise radial distribution function (RDF) for each state in Supplementary Fig. 3 to demonstrate the average configurations, which is consistent with above description. We also analyze the total charge carried by the Li-ions in each state considering their solvation environments in Fig. 4(c) (see Supplementary Note 3 and Supplementary Table 1 for details). Interestingly, both state 0 and state 1 carry almost zero total charge in their first solvation shell due to the one TFSI anion in their solvation environments.

We further study the transition between the four Li-ion states. Three relaxation processes are identified in the dynamical model as shown in Fig. 4(d, e). By analyzing the eigenvectors, we learn that the slowest relaxation is a process involving the transport of a Li-ion into and out of a PEO coordinated environment. The second slowest relaxation happens mainly between state 0 and state 1, corresponding to a movement of the hydroxyl group. The transitions from state 0 to states 2 and 3 constitute the last relaxation process, as state 0 can be thought of an intermediate state between state 2 and state 3. The model performs well in CK tests (Fig. 4(f)). Relaxation processes in the PEO/LiTFSI systems have been extensively studied experimentally Mao et al. (2000); Do et al. (2013), but it is difficult to pinpoint the exact atomic scale dynamics related to these relaxations. The dynamical model learned by GDyNet provides additional insights into the understanding of Li-ion transport in polymer electrolytes.

Implications to lithium ion conduction. The state configurations and dynamical model allow us to further quantify the transitions that are responsible for the Li-ion conduction. In Fig. 5, we compute the contribution from each state transition to the Li-ion conduction using the Koopman model at $\tau=$0.8\text{,}\mathrm{n}\mathrm{s}$$. First, we learn that the majority of conduction results from transitions within the same states ($i\rightarrow i$). This is because the transport of Li-ions in PEO is strongly coupled with segmental motion of the polymer chains Borodin and Smith (2006); Diddens et al. (2010), in contrast to the hopping mechanism in inorganic solid electrolytes Bachman et al. (2015). In addition, due to the low charge carried by state 0 and state 1, the majority of charge conduction results from the diffusion of states 2 and 3, despite their relatively low populations. Interestingly, the diffusion of state 2, a negatively charged species, accounts for $\sim 40\%$ of the Li-ion conduction. This provides an atomic scale explanation to the recently observed negative transference number at high salt concentration PEO/LiTFSI system Pesko et al. (2017).

Discussion

We have developed a general approach, Graph Dynamical Networks (GDyNets), to understand the atomic scale dynamics in material systems. Despite being widely used in biophysics Husic and Pande (2018), fluid dynamics Mezić (2013), and kinetic modeling of chemical reactions Georgiev et al. (2005); Valor et al. (2013); Miller and Klippenstein (2006), Koopman models, (or Markov state models Husic and Pande (2018), master equation methods Buchete and Hummer (2008); Sriraman et al. (2005)) have not been used in learning atomic scale dynamics in materials from MD simulations except for a few examples in understanding solvent dynamics Gu et al. (2013); Hamm (2016); Schulz et al. (2018). Our approach also differs from several other unsupervised learning methods Cubuk et al. (2016); Nussinov et al. (2016); Kahle et al. (2019) by directly learning a linear Koopman model from MD data. Many crucial processes that affect the performance of materials involve the local dynamics of atoms or small molecules, like the dynamics of lithium ions in battery electrolytes Funke (1993); Xu (2004), the transport of water and salt ions in water desalination membranes Corry (2008); Cohen-Tanugi and Grossman (2012), the adsorption of gas molecules in metal organic frameworks Rowsell et al. (2005); Li et al. (2009), among many other examples. With the improvement of computational power and continued increase in the use of molecular dynamics to study materials, this work could have broad applicability as a general framework for understanding the atomic scale dynamics from MD trajectory data.

Compared with the Koopman models previously used in biophysics and fluid dynamics, the introduction of graph convolutional neural networks enables parameter sharing between the atoms and an encoding of local environments that is invariant to permutation, rotation, and reflection. This symmetry facilitates the identification of similar local environments throughout the materials, which allows the learning of local dynamics instead of exponentially more complicated global dynamics. In addition, it is easy to extend this method to learn global dynamics with a global pooling function Xie and Grossman (2018a). However, a hierarchical pooling function is potentially needed to directly learn the global dynamics of large biological systems including thousands of atoms. It is also possible to represent the local environments using other symmetry functions like smooth overlap of atomic positions (SOAP) Bartók et al. (2013), social permutation invariant (SPRINT) coordinates Pietrucci and Andreoni (2011), etc. By adding a few layers of neural networks, a similar architecture can be designed to learn the local dynamics of atoms. However, these built-in invariances may also cause the Koopman model to ignore dynamics between symmetrically equivalent structures which might be important to the material performance. One simple example is the flip of ammonia molecule – the two states are mirror symmetric to each other so the GCN will not be able to differentiate them by design. This can potentially be resolved by partially break the symmetry of GCN based on the symmetry of the material systems.

The graph dynamical networks can be further improved by incorporating ideas from both the fields of Koopman models and graph neural networks. For instance, the auto-encoder architecture Lusch et al. (2018); Wehmeyer and Noé (2018); Ribeiro et al. (2018) and deep generative models Wu et al. (2018) start to enable the direct generation of future structures in the configuration space. Our method currently lacks a generative component, but this can potentially be achieved with a proper graph decoder Jin et al. (2018); Simonovsky and Komodakis (2018). Furthermore, transfer learning on graph embeddings may reduce the number of MD trajectories needed for learning the dynamics Sultan and Pande (2017); Altae-Tran et al. (2017).

In summary, graph dynamical networks present a general approach for understanding the atomic scale dynamics in materials. With a toy system of lithium ion transporting in a face-centered cubic lattice, we demonstrate that learning local dynamics of atoms can be exponentially easier than global dynamics in material systems with representative local structures. The dynamics learned from two more complicated systems, solid-liquid interfaces and solid polymer electrolytes, indicate the potential of applying the method to a wide range of material systems and understanding atomic dynamics that are crucial to their performances.

Methods

Construction of the graphs from trajectory. A separate graph is constructed using the configuration in each time step. Each atom in the simulation box is represented by a node $i$ whose embedding $\bm{v}_{i}$ is initialized randomly according to the element type. The edges are determined by connecting $M$ nearest neighbors whose embedding $\bm{u}_{(i,j)}$ is calculated by,

[TABLE]

where $\mu_{t}=t\cdot$0.2\text{,}\mathrm{\SIUnitSymbolAngstrom} for $t=0,1,...,K$, $\sigma=$0.2\text{\,}\mathrm{\SIUnitSymbolAngstrom}, and $d_{(i,j)}$ denotes the distance between $i$ and $j$ considering the periodic boundary condition. The number of nearest neighbors $M$ is 12, 20, and 20 for the toy system, Si-Au binary system, and PEO/LiTFSI system, respectively.

Graph convolutional neural network architecture details. The convolution function we employed in this work is similar to those in refs. Xie and Grossman (2018a, b) but features an attention layer Velickovic et al. (2017). For each node $i$, we first concatenate neighbor vectors from last iteration $\bm{z}_{(i,j)}^{(t-1)}=\bm{v}_{i}^{(t-1)}\oplus\bm{v}_{j}^{(t-1)}\oplus\bm{u}_{(i,j)}$, then we compute the attention coefficient of each neighbor,

[TABLE]

where $\bm{W}_{\mathrm{a}}^{(t-1)}$ and $b_{\mathrm{a}}^{(t-1)}$ denotes the weights and biases of the attention layers and the output $\alpha_{ij}$ is a scalar number between 0 and 1. Finally, we compute the embedding of node $i$ by,

[TABLE]

where $g$ denotes a non-linear ReLU activation function, and $\bm{W}_{\mathrm{n}}^{(t-1)}$ and $\bm{b}_{\mathrm{n}}^{(t-1)}$ denotes weights and biases in the network.

The pooling function computes the average of the embeddings of each atom for the set of target atoms,

[TABLE]

Determination of the relaxation timescales. The relaxation timescales represent the characteristic timescales implied by the transition matrix $\bm{K}(\tau)$, where $\tau$ denotes the lag time of the transition matrix. By conducting an eigenvalue decomposition for $\bm{K}(\tau)$, we could compute the relaxation timescales as a function of lag time by,

[TABLE]

where $\lambda_{i}(\tau)$ denotes the $i$th eigenvalue of the transition matrix $\bm{K}$. These $t_{i}(\tau)$ are plotted in Figs. 2(b), 3(d), and 4(d) for each material system. If the dynamics of the system is Markovian, i.e. Eq. (6) holds, one can prove that the relaxation timescales $t_{i}(\tau)$ will be constant for any $\tau$. Mardt et al. (2018); Wu and Noé (2017) Therefore, we select a smallest $\tau*$ from Figs. 2(b), 3(d), and 4(d) to obtain a dynamical model that is Markovian and contains most dynamical details. We then compute the relaxation timescales using this $\tau*$ for each material system, and these timescales remain constant for any $\tau>\tau*$.

State-weighted radial distribution function. Radial distribution function (RDF) describes how particle density varies as a function of distance from a reference particle. RDF is usually determined by counting the neighbor atoms at different distances over MD trajectories. We calculate the RDF of each state by weighting the counting process according to the probability of the reference particle being in state $i$,

[TABLE]

where $r_{\mathrm{A}}$ denotes the distance between atom A and the reference particle, $p_{i}$ denotes the probability of reference particle being in state $i$, and $\rho_{i}$ denotes the overall density of state $i$.

Analysis of Li-ion conduction. We first compute the expected mean-squared-displacement of each transition at different $t$ using the Bayesian rule,

[TABLE]

where $p_{i}(t)$ is the probability of state $i$ at time $t$, and $d^{2}(t^{\prime},t^{\prime}+t)$ is the mean-squared-displacement between $t^{\prime}$ and $t^{\prime}+t$. Then, the diffusion coefficient of each transition $D_{i\rightarrow j}(\tau)$ at the lag time $\tau$ can be calculated by,

[TABLE]

which is shown in Supplementary Table 2.

Finally, we compute the contribution of each transition to Li-ion conduction with Koopman matrix $\bm{K}(\tau)$ using the cluster Nernst-Einstein equation France-Lanord and Grossman (2018),

[TABLE]

where $e$ is the elementary charge, $k_{\mathrm{B}}$ is the Boltzmann constant, $V$, $T$ are the volume and temperature of the system, $N_{\mathrm{Li}}$ is the number of Li-ions, $\pi_{i}$ is the stationary distribution population of state $i$, and $z_{ij}$ is the averaged charge of state $i$ and state $j$. The percentage contribution is computed by,

[TABLE]

Lithium diffusion in the FCC lattice toy system. The molecular dynamics simulations are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Plimpton (1995), as implemented in the MedeA®Med (2018) simulation environment. A purely repulsive interatomic potential in the form of a Born-Mayer term was used to describe the interactions between Li particles and the S sublattice, while all other interactions (Li–Li and S–S) are ignored. The cubic unit cell includes one Li atom and four S atoms, with a lattice parameter of $6.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, a large value allowing for a low barrier energy. $200\text{\,}\mathrm{n}\mathrm{s}$ MD simulations are run in the canonical ensemble (nVT) at a temperature of $64\text{\,}\mathrm{K}$, using a timestep of 1 fs, with the S particles frozen. The atomic positions, which constituted the only data provided to the GDyNet and VAMPnet models, are sampled every $0.1\text{\,}\mathrm{p}\mathrm{s}$. In addition, the energy following the Tet-Oct-Tet migration path was obtained from static simulations by inserting Li particles on a grid.

**Silicon dynamics at solid-liquid interface. **The molecular dynamics simulation for the Si-Au binary system was carried out in LAMMPS Plimpton (1995), using the modified embedded-atom method interatomic potential Wang et al. (2017); Ryu and Cai (2010). A sandwich like initial configuration was created, where Si-Au liquid alloy was placed in the middle, contacting with two {110} orientated crystalline Si thin films. $25\text{\,}\mathrm{n}\mathrm{s}$ MD simulations are run in the canonical ensemble (nVT) at the eutectic point (629K, 23.4% Si Ryu and Cai (2010)), using a time step of 1 fs. The atomic positions, which constituted the only data provided to the GDyNet model, are sampled every 20 ps.

Scaling of the algorithm. The scaling of the algorithm is $\mathcal{O}(NMK)$, where $N$ is the number of atoms in the simulation box, $M$ is the number of neighbors used in graph construction, and $K$ is the depth of the neural network.

Data availability. The MD simulation trajectories of the toy system, the Si-Au binary system, and the PEO/LiTFSI system are available at https://archive.materialscloud.org/2019.0017.

Code availability. GDyNets is implemented using TensorFlow Abadi et al. (2016) and the code for the VAMP loss function is modified on top of ref. Mardt et al. (2018). The code is available from https://github.com/txie-93/gdynet.

Additional information

Acknowledgements This work was supported by Toyota Research Institute. Computational support was provided by Google Cloud, the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and the Extreme Science and Engineering Discovery Environment, supported by National Science Foundation grant number ACI-1053575.

Author contributions. T.X. developed the software and performed the analysis. A.F.L. and Y.W. performed the molecular dynamics simulations. T.X., A.F.L., Y.W., Y.S.H., and J.C.G. contributed to the interpretation of the results. T.X. and J.C.G. conceived the idea and approach presented in this work. All authors contributed to the writing of the paper.

Competing interests. The authors declare no competing interests.

Materials & Correspondence. Correspondence and requests for materials should be addressed to J.C.G. (jcg@mit.edu).

Supplementary Information: Graph Dynamical Networks for Unsupervised Learning of Atomic Scale Dynamics in Materials

Xie et al.

Supplementary Notes

Supplementary Note 1: gated architecture in graph convolutional neural networks. The gated architecture in GCN is important for learning the local environments in a complex material system. In general, it is challenging to define bonds, i.e. the topology of the graph, in materials to capture the non-covalent interactions which may affect the atomic dynamics. We resolve this challenge by introducing the gated architecture in GCN to reweigh the strength of each connection. In the graph construction, we connect each atom with its $M$ nearest neighbors, where the $M$ is large enough to include non-covalent interactions. During the training, the gated architecture (or the attention layer as described in the methods section) automatically learns a weight factor that is related to the bond length and atom types. This weight factor reweighs the importance of the $M$ nearest neighbors to center atom. As a result, we can learn a representation of the local environments in complex materials using the graphs constructed by the simple nearest neighbor approach.

Supplementary Note 2: computation of global dynamics from local dynamics in the toy system. To compute the global dynamics from local dynamics, we first assume that the transition matrix of the local Koopman model has the form,

[TABLE]

where $p_{o}$ and $p_{t}$ denotes the probability of the lithium atom staying in the octahedral and tetrahedral sites, respectively. Since there are 4 octahedral sites and 8 tetrahedral sites that are connected to each other in the FCC lattice, we can write the transition matrix of the global Koopman model as,

[TABLE]

By computing the eigenvalues of $\bm{K}_{\mathrm{global}}$, we could obtain the relaxation timescales and understand the global dynamics of the toy system. There are two major reasons for the discrepancy between the computed and observed global Koopman model: 1) the amount of MD data is not large enough to capture of full global dynamics of the lithium atom by directly learning a global dynamical model; 2) the probability of lithium atom transporting to nearby sites of the same type is not strictly zero at a given $\tau$, so the $p_{o}$ and $p_{t}$ in $\bm{K}_{\mathrm{local}}$ and the zero terms in $\bm{K}_{\mathrm{global}}$ are approximate.

Supplementary Note 3: determination of charge for each state in PEO/LiTFSI. The charge carried by each state is determined by computing the charge integral within the first solvation shell of the Li-ions. We perform a Gaussian curve fit using the state-weighted radial distribution function of nitrogen in Supplementary Figure 3(c), since the nitrogen is the center of the TFSI anion. We assume the edge of the first solvation shell as the mean of the Gaussian curve plus 3 sigma. The charge carried by each state is computed by integrating the charge density within the first solvation using,

[TABLE]

where $r^{*}$ denotes the edge of the first solvation shell, and $g_{i}(r_{N})$ and $g_{i}(r_{Li})$ denote the state-weighted radial distribution functions of nitrogen and lithium, respectively. The resulting charge for each state is summarized in Supplementary Table 1.

Supplementary Tables

Supplementary Figures