# Evaluating the Transferability of Machine-Learned Force Fields for Material
Property Modeling

Shaswat Mohanty111Equal contribution SangHyuk Yoo222Equal contribution
Keonwook Kang  Wei Cai<EMAIL_ADDRESS>Department of Mechanical
Engineering, Stanford University, CA 94305-4040, USA  School of Mechanical
Engineering, Yonsei University, Seoul 03722, Republic of Korea

###### Abstract

Machine-learned force fields have generated significant interest in recent
years as a tool for molecular dynamics (MD) simulations, with the aim of
developing accurate and efficient models that can replace classical
interatomic potentials. However, before these models can be confidently
applied to materials simulations, they must be thoroughly tested and
validated. The existing tests on the radial distribution function and mean-
squared displacements are insufficient in assessing the transferability of
these models. Here we present a more comprehensive set of benchmarking tests
for evaluating the transferability of machine-learned force fields. We use a
graph neural network (GNN)-based force field coupled with the OpenMM package
to carry out MD simulations for Argon as a test case. Our tests include
computational X-ray photon correlation spectroscopy (XPCS) signals, which
capture the density fluctuation at various length scales in the liquid phase,
as well as phonon density-of-state in the solid phase and the liquid-solid
phase transition behavior. Our results show that the model can accurately
capture the behavior of the solid phase only when the configurations from the
solid phase are included in the training dataset. This underscores the
importance of appropriately selecting the training data set when developing
machine-learned force fields. The tests presented in this work provide a
necessary foundation for the development and application of machine-learned
force fields for materials simulations.

###### keywords:

Graph Neural Network , Machine-learned Force Field , X-ray photon correlation
spectroscopy , Optical contrast , Phonon Density of States

††journal: Computer Physics Communications

## Highlights

  * 1.

Development of a suite of benchmarking tests, including X-ray photon
correlation spectroscopy for machine-learned force fields (code available).

  * 2.

Phonon density of states and melting point calculations are necessary tests
for the solid phase.

  * 3.

Training data should include both solid and liquid configurations to
accurately model material behavior.

##  1 Introduction

Molecular dynamics (MD) is a powerful tool for studying the equilibrium and
transport properties of many-body systems with applications in diverse fields
such as materials science [1, 2, 3], polymer chemistry [4, 5, 6], biochemistry
[7] and medical science [8, 9, 10]. In MD simulations, the motion of the atoms
is driven by interatomic forces, which are evaluated using interatomic
potentials. However, there are limitations to traditional molecular dynamics
simulations, including a limited length scale, a limited time scale, and
insufficient force field accuracy. In recent years, machine learning
techniques have been used to address these limitations, with the goal of
accelerating molecular dynamics simulations and improving force field
accuracy. Machine learning has been used to predict atomic trajectories [11,
12] with the ultimate aim of bridging the timescale gap between MD simulations
and experiments [13, 14]. However, these techniques often require retraining
the model for different system sizes, making it challenging to apply to large
systems, which are often necessary for property predictions. As a result,
there is significant interest in developing scalable machine-learning methods
for MD simulations.

Graph-based deep learning techniques have proven effective in developing
accurate and scalable force fields that can be trained on data from small
simulation cells and then applied to larger systems. Graph neural networks
(GNN) can accurately predict atomic forces with a computational cost that
scales linearly with the number of atoms, making them well-suited for large-
scale simulations. Recently, GNN-based force fields have been used to predict
MD trajectories using approaches such as SchNet [15], and DeepMD [11, 16, 17].
These force fields, when trained on high-quality _ab initio_ data, can perform
large-scale simulations at an accuracy comparable to that of _ab initio_ model
but only at a small fraction of the cost [18, 19, 20]. However, before GNN
force fields (or any machine-learned force fields) can be confidently
deployed, their transferability to configurations beyond the training data set
must be established. Currently, most existing benchmark tests on machine-
learned force fields are limited to the radial distribution function and self-
diffusivity in the liquid phase. As we show below, these tests are
insufficient for establishing their reliability for materials modeling.

To this end, we present a series of benchmarking tests, comparing the
predictions of GNN-based force fields [21] to those of the original model that
produces the training data. As an example, the original model is the classical
Lennard-Jones interatomic potential for Argon. The purpose of choosing such a
simple model as a benchmark in this work is to provide a baseline to evaluate
the transferability of machine-learned force fields without incurring
significant computational costs. Using such a simple model, a large amount of
data can be easily generated for training and comparison purposes. The
findings from this study pave the road for future tests for machine-learned
force fields fitted to _ab initio_ data which are much more expensive to
generate. In addition to the radial distribution function and self-
diffusivity, we compare the computational X-ray photon correlation
spectroscopy (XPCS) signals [22], which capture density fluctuations at
different length scales in the liquid phase. We further show that a model
trained exclusively on liquid configurations fails to accurately capture the
vibrational frequency distribution in the solid phase. This deficiency is
fixed only after the model is trained on configurations sampled from both the
liquid and solid phases. These findings raise a concern about the
transferability of the GNN-based force field and underscore the importance of
ensuring that the training data adequately cover the areas of interest in the
configurations space.

The paper is structured as follows. In Section 2 we discuss the preparation of
training data, the GNN model details, and the hyperparameters used for
training the model. Additionally, we present the metrics that will be used to
analyze the performance of the GNN-MD simulations against the MD simulations
using the original model. In Section 3, we present the results of these
benchmark tests. We summarize our findings in Section 4.

##  2 Model details

GNN-based models provide a linearly scalable force field that has the
potential to be used in _ab initio_ molecular dynamics (AIMD) simulations.
However, verifying a model against direct _ab initio_ simulations are
computationally expensive, especially for larger systems. To address this
issue, we will focus on conducting benchmarking tests on a GNN force field
that is trained on forces generated using the classical Lennard-Jones
interatomic potential. This will allow us to analyze the model’s performance
without the added computational cost of using _ab initio_ data. We assume that
the transferability of the GNN model can be assessed when trained on forces
from either classical interatomic potentials or _ab initio_ data, with
comparable results.

###  2.1 Preparation of reference data

We prepared a 256-atom configuration of Argon (Ar), described by the Lennard-
Jones (LJ) interatomic potential used in our earlier work [22], by using
LAMMPS [23].

$\displaystyle V_{\rm LJ}(r)$
$\displaystyle=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$
(1)

where $\epsilon=0.0103$ eV, $\sigma=3.40$ Å. The density of the system used is
$\rho$ = 0.858 amu/Å3 (0.844 in LJ unit system). The cut-off radius of the
neighbor list is 8.5 Å (2.5 in LJ unit system). Periodic Boundary Conditions
(PBC) were applied in all three directions. The initial positions of atoms
were generated according to a perfect face-centered cubic (FCC) lattice and
the velocities were assigned random values corresponding to a temperature of
10 K. To reach the state of thermal equilibrium, the Nosé-Hoover NPT ensemble
and the NVT ensemble were applied sequentially. The equilibration simulation
was performed for about 2.156 ns, with a timestep of 10.78 fs. Simulations
under both ensembles used a thermostat with a collision frequency of 0.02 /ps,
a chain length of 5, and 5 MTK [24] loops were used. For pressure control
under the NPT ensemble, a barostat with a collision frequency of 0.2 /ps was
used. After implementing this protocol of _initial equilibration_ , the
thermostat temperature (NVT simulation) was increased with a linear ramp from
10 K to 105 K over the 2.156 ns trajectory, during which atomic configurations
were saved every 100 timesteps. This same process is repeated ten times with
different initialization of velocities to generate a total of 20000
configurations. Additionally, the configuration at the end of the _initial
equilibration_ protocol is used as the initial state to run the analysis
trajectory by using LAMMPS for the MD results and the Atomic Simulation
Environment [25] and OpenMM [26] script for the GNN-MD simulation,
respectively.

###  2.2 GNN model

Our benchmarking tests are carried out on the GNN-based force field following
the GNN structure and training procedure of Li et al. [21]. For completeness,
we present a brief description of the GNN model in this section. For the GNN-
based force predictor, the atomic configurations are first converted into
graphs that encode the local atomic environment. Here, each atom acts as a
node, and the vectors towards neighboring atoms within a cut-off radius act as
edges in the graph. The node information includes the atomic position and
force vector. The edge information includes the distance between neighboring
atoms and the unit vector toward the neighbors. A brief description of each
step in the GNN force field prediction is enlisted below. A more detailed
description of the model can be found in [21].

  * 1.

_Encoding layer_
The feature $v_{i}^{(l)}$ is an array of size $h$ and represents the encoded
node feature at the $l^{\rm th}$ message-passing layer. The reference node
encoding, $v_{i}^{(0)}$, is created by passing the atomic species information
$s_{i}$ through an encoding multi-layer perceptron (MLP), $e^{N}$, such that
$v_{i}^{(0)}=e^{N}(s_{i})$. The feature $e_{ij}^{(l)}$ is also an array of
size $h$ represents the encoded edge feature at the $l^{\rm th}$ message-
passing layer. The reference edge feature is also created through an MLP,
$e^{E}$, where $e_{ij}^{(l)}=e^{E}(\mathbf{q}_{ij},d_{ij})$. Here, for a given
pair of atoms[$i$,$j$], the unit vector and the distance between them are
denoted by $\mathbf{q}_{ij}$ and $d_{ij}$, respectively. The superscript $l$
is the index of the message-passing layer, which varies from 0 to $n$ and
denotes how many layers the features have passed through. At the encoding
step, $d_{ij}$ is represented by the radial basis functions, following the
implementation of Schnet [15].

  * 2.

_Message passing layer_
The encoded features, $v_{i}^{(l)}$ and $e_{ij}^{(l)}$, are passed through the
message passing layer $m_{j\rightarrow{i}}^{(l)}$. The equation of the
message-passing layer is as follows,

$m_{j\rightarrow{i}}^{(l)}=\Phi^{(l)}\left(v_{j}^{(l-1)}+e_{ij}^{(l-1)}+v_{i}^{(l-1)}\right)\odot
v_{j}^{(l-1)},\qquad\forall j\in\mathcal{N}(i),$ (2)

where $\Phi^{(l)}$ denotes a MLP with the activation function being the
Gaussian Error Linear Units (GELU) function [27]. As an input of $\Phi^{(l)}$,
the edge features and node features are added up as a vector sum since their
dimensions are the same. $\odot$ refers to elemental multiplication, and $j$
is the index of the set of neighbors $\mathcal{N}(i)$ of the $i^{\rm th}$
node.
Each message from $j$-th atom aggregates into $M_{i}^{(l)}$ and the new node
feature $v_{i}^{(l)}$ is computed by the following equations.

$\displaystyle M_{i}^{(l)}=$ $\displaystyle\sum_{\forall
j\in\mathcal{N}(i)}m_{j\rightarrow{i}}^{(l)},$ (3) $\displaystyle
v_{i}^{(l)}=$
$\displaystyle\,\Theta^{(l)}\left(v_{i}^{(l-1)}+M_{i}^{(l)}\right)+v_{i}^{(l-1)}.$
(4)

Here, $\Theta^{(l)}$ is the node update MLP function at the $l^{\rm th}$
layer. Eq. (4) shows that the node features are updated recursively. However,
the edge features are not updated recursively [21]. Instead, the edge features
at every level are computed directly from the reference edge encoding through
an MLP function, $A^{(l)}$, i.e., $e_{ij}^{(l)}=A^{(l)}(e_{ij}^{(0)})$.

  * 3.

_Decoding layer_
At the last step of the neural network, the updated node feature is decoded
into the interatomic forces $\bm{f}_{i}$ by executing a forward pass through
the decoding MLP.

Although we present our benchmarking tests on a model trained on pair
potential, the GNN model is not constrained to a force field derived from a
pair potential form. The GNN architecture used here is also capable of
learning a force field derived from _ab initio_ data [21].

###  2.3 Details of the training

We set the number of encoding features and latent features as $h=128$ and used
$n=4$ message-passing layers in the GNN. The cut-off radius for converting
atomic configuration to graph is 7.5 Å, which is 1 Å smaller than that in the
LJ potential. The effects of the cutoff radius will be studied in the future.
The loss function for the training contains the L1 distance between the
reference force $\bm{f}_{i}$ and the predicted force $\hat{\bm{f}_{i}}$, as
well as a penalty function, as follows.

$\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\left\lVert\bm{f}_{i}-\hat{\bm{f}_{i}}\right\rVert_{1}+\lambda\left\lVert\frac{1}{N}\sum_{i=1}^{N}\hat{\bm{f}_{i}}\right\rVert_{1},$
(5)

where $N$ is the number of atoms in the reference data. The force on the
center of mass from an interatomic potential is expected to be zero. However,
the GNN model does not compute the forces from the derivative of a potential
function, we are not guaranteed a zero force on the center of mass. We
minimize the magnitude of the total force on the entire system to prevent
divergent dynamics of the center of mass by adding a penalty function to the
loss function. The regularization parameter, $\lambda$, is set to 0.01 during
the training. As described earlier, we used 20,000 atomic configurations as
our training dataset. The dataset is shuffled and split into the train and
test set with a 90:10 ratio. During training, we carry out standard
normalization of the atomic forces such that the mean atomic force is zero and
the variance is unity. We used PyTorch 1.11.0, PyTorch-lightning 1.6.3, DGL
0.8.1 and Scikit-learn 0.24.2 packages to train and test the model. We use the
Adam optimizer with an exponential learning rate scheme (from $3\times
10^{-4}$ to $1\times 10^{-7}$) during the training of the force field. The
total number of epochs is set to 30 and we choose the trained model parameter
(weights and biases) at the last epoch as the force calculator for our GNN-MD
simulation.

###  2.4 Metrics for model performance

Here we list the static and dynamic metrics that will be used to compare the
performance of the GNN-MD simulations against the MD simulations. To examine
the structure of the liquid phase we will study the pair distribution
function. To examine the dynamics of the liquid phase, we evaluate the self-
diffusivity and carry out the computational XPCS analysis on the MD
trajectory. The relations to compute the properties such as the radial
distribution function, the structure factor, and the computational XPCS
signal, are described only briefly in A.2 and A.3 (see [22] for more details),
since the focus of this paper is to establish the accuracy of the GNN-MD in
capturing the dynamics and structure of the simulation. To evaluate the model
performance in the solid phase, we compute the phonon density of states
(PDOS), and the melting point by using the interface method. Further details
on the key equations we solve to obtain the numerical results are given in A.

##  3 Numerical Results

The MD and GNN-MD simulations discussed in this section are carried out on a
256-atom and a 4000-atom system to test the scalability and transferability of
the GNN force field which was trained on 256-atom configurations. The training
samples for the 256 and 4000 atoms system are obtained from MD simulations
using LAMMPS [23]. While in the previous work [21] only liquid configurations
are included in the training data, here we train two models: Model A is
trained only using liquid configurations and Model B is trained using both
liquid and solid configurations (see Section 2.1 for the procedure of
preparing these configurations). The results presented in Sections 3.1, 3.2,
and 3.3, corresponds to Model B. However, Model A gives similar results for
the tests in these sections; hence they are omitted here for brevity. On the
other hand, Model A and Model B give different results for some of the tests
involving the solid phase, and both will be presented in Section 3.4.

The initial configuration for the testing MD and GNN-MD simulations are
prepared using the _initial relaxation_ protocol. The final trajectory for the
analysis is then run using the same NVT ensemble as in the final relaxation
step. However, the total simulation time extends to 539 ps with a timestep of
10.78 fs. The simulation frames are stored at 107.8 fs intervals for our
analysis. In comparison to MD simulations using the LJ potential, the GNN-MD
is $10^{1}$ times slower, because the atomic forces from the LJ potential can
be computed very quickly. We reiterate that the focus of this work is to use
the LJ potential as an example to test the transferability of the GNN force
field in capturing the static and dynamic properties of the material system
that it is trained on.

###  3.1 Accuracy of force prediction

We test the atomic forces predicted atomic forces against the calculated
atomic forces from the Lennard-Jones potential for 1,000 sampled
configurations of our 256-atom simulation. These configurations are different
from the ones used for the training of the GNN force field. We see that the
predicted forces are well correlated with the actual forces across all atoms
and all configurations, as shown in Fig. 1 (parity plot). The coefficient of
determination, R2, is almost $1.0$ ($1-$R${}^{2}\approx 10^{-4}$) which shows
how well correlated the predicted force is to the actual forces computed from
LAMMPS by using the Lennard-Jones force field. The generalizability of the
force field to a 4000-atom system is captured by the parity plot shown in Fig.
1, where the quality of the force prediction remains qualitatively the same.

Figure 1: Parity plot for the predicted forces and the actual forces for the
(a) 256 atom and (b) 4000 atom system. (c) Distribution of the absolute error
of the predicted forces for the 256 atom and 4000 atom system over the
configurations of analysis.

Figure 1 shows the distribution of the absolute error of the predicted forces,
for the 256-atom (distribution of atomic force components over 1,000
configurations) and the 4000-atom (distribution of atomic force components
over 100 configurations) system. The standard deviation of the force error is
0.15 meV/Å. We see that the absolute error distribution is the same
irrespective of the system size, which indicates the error of the trained
model and is independent of the size of the system. This establishes the
quality of the trained model consistent with previous report [21].

###  3.2 Structure of liquid phase

Here we examine the orientation-averaged radial distribution function, $g(r)$,
as defined in A.2. We average the $g(r)$ over 100 configurations that are
sampled every 5.34 ps from our simulation. Figure 2 shows that the $g(r)$
obtained from the MD and GNN-MD simulation are in close agreement for the
(small) system it is trained on, confirming the previous observation [21]. We
further find that $g(r)$ between the MD and GNN-MD simulation agree well with
each other even for the larger simulation cell containing 4000 atoms, as shown
in Fig. 2. We also observe that the orientation-averaged structure factor,
$S(q)$ (not shown here), which can be obtained from the Fourier transform of
$g(r)$, agrees between MD and GNN-MD simulations, for both the 256-atom and
4000-atom systems.

Figure 2: The $g(r)$ over 100 configurations for the (a) 256 atoms and (b)
4000 atoms MD and GNN-MD simulation.

###  3.3 Dynamics of liquid phase

Given that atomic trajectories in MD simulations are chaotic, small
differences in atomic forces can lead to large divergences in the atomic
trajectories at a later time. Therefore, we don’t expect the atomic
trajectories in MD and GNN-MD to agree well with each other. However, for the
GNN-based force field to be useful, the statistical properties of the atomic
trajectories need to be consistent with those of the original MD model. A
commonly tested statistical property is the self-diffusivity obtained from the
mean-squared displacement (MSD) as described in A. Figure 1 shows that the MSD
predicted by the GNN model is consistent with that from the original MD model.
For the 256-atom system, the MSD from GNN-MD agrees very well with the MD
model for the first $150$ ps, after which the difference grows somewhat
larger, as shown in Fig. 3. For the 4000-atom system, the MSD from GNN-MD
agrees very well with the MD model for the first $300$ ps, after which the
difference also grows somewhat larger, as shown in Fig. 3.

Figure 3: The mean-squared displacement (MSD) for the (a) 256 atoms and the
(b) 4000 atoms from the MD and GNN-MD simulation over the 539 ps trajectory at
100 K.

The MSD curves from the entire ($\sim$ 539 ps) trajectories are fitted to
straight (dashed) lines to obtain the predicted self-diffusivity $D$ (see Fig.
3). The agreement between MD and GNN-MD is within $10\%$ for the 256-atom
system, and within $5\%$ for the 4000-atom system, at the temperature of 100
K. The level of agreement for the 4000-atom system is maintained at all
temperatures from 95 K to 110 K, and appears to improve with increasing
temperature, as shown in Table 1.

Temperature |  $D$ ($\mu$m2/s)
---|---
MD | GNN-MD | Relative Error ($\%$)
95 K | 3254.70 | 3096.14 | 4.87
100 K | 3673.08 | 3802.89 | 3.53
110 K | 5010.12 | 4979.92 | 0.60
Table 1: The self-diffusivity calculated from the MSD for the 4000 atom system
at different temperatures using the MD and GNN-MD simulation. The MSD for all
temperatures can be found in Fig. 11.

In addition to the bulk diffusivity, we also test the GNN-based force field by
carrying out the XPCS analysis. This allows us to quantify the density
fluctuations arising from atomic motion at different length scales. Figure 4
shows two computed XPCS speckle patterns from snapshots of the MD and GNN-MD
trajectories, respectively. The speckle patterns are computed for all saved
configurations ($\sim$ 107.8 fs apart) in the atomic trajectories so that the
time correlation at each pixel can be obtained.

Figure 4: The speckle patterns of scattered intensity $I(\bm{q})$ on a
spherical slice ($81\times 81$ pixels) in the $\bm{q}$-space for (a) MD and
(b) GNN-MD simulation.

Figure 5 shows that the decay of the time correlation of the speckle averaged
over all wave-vectors close to a given $q=|\bm{q}|$, agrees well between MD
and GNN-MD, for both the 256-atom system and the 4000-atom system. If the
system dynamics is purely diffusive, then $g_{2}(q,\tau)$ should decay
exponentially.

Figure 5: The $g_{2}(q,\tau)$ at $q=1.844\pm 0.029$ Å-1 over the first 3.5574
ps for the (a) 256 atoms and the (b) 4000 atoms MD and GNN-MD simulation at
100 K.

Figure 5 compares the time correlation for $q=1.844\pm 0.029$ Å-1 (probing the
length scale corresponding to the first nearest neighbor). We have also
carried out the same computation for all wave-vector magnitudes from
$q=0.461\pm 0.029$ Å-1 to $q=1.844\pm 0.029$ Å-1 and observed close agreement
between MD and GNN-MD models. In all cases, we see that the $g_{2}(q,t)$
decays as a single exponential after the initial sub-diffusive dynamics
associated with caging effects [22], as is expected for liquid Ar.

Figure 6: The $\Gamma(q)$ as a function of $q^{2}$ for the (a) 256 atoms and
the (b) 4000 atoms MD and GNN-MD simulation at 100 K.

From $g_{2}(q,\tau)$, we extract the rate of exponential decay $\Gamma(q)$ in
the long time limit, which reveals the dispersion relation in the liquid [22].
Figure 6 shows that the $\Gamma(q)$ function obtained from MD and GNN-MD are
in good agreement. Interestingly, the agreement appears better in the
$4000$-atom system than in the $256$-atom system on which the GNN-base force
field is trained. Additionally, we found that the dispersion relation from the
GNN-MD simulation agrees well with MD simulation results for all temperatures
of 95 K, 100 K, and 110 K (see Fig. 13).

Figure 7: The optical contrast $\beta_{\Delta}(q)$ as a function of X-ray
pulse width ($\Delta_{\tau}$) for the (a) 256 atoms and the (b) 4000 atoms MD
and GNN-MD simulation for $\lvert\bm{q}\rvert=1.844\pm 0.029$ Å-1 at 100 K.

In addition to the dynamics obtained from the XPCS speckle fluctuations, we
also test the statistical distribution of the X-ray speckles by evaluating the
optical contrast in the $q$-range $\lvert\bm{q}\rvert=1.844\pm 0.029$ Å-1 for
varying pulse widths, $\Delta_{\tau}$, of the incident X-ray. Figure 7 shows
that the optical contrast obtained from MD and GNN-MD models are in close
agreement. The optical contrast $\beta_{\Delta}(q)$ is expected to decrease
with increasing pulse width, $\Delta_{\tau}$, where the time at which
$\beta_{\Delta}(q)$ decreases to half of its maximum value gives an estimate
of the correlation time. Table 2 lists the correlation time obtained at three
temperatures, where the estimates from MD and GNN-MD models agree well with
each other (time resolution of 0.11 ps).

Temperature | Correlation time (ps)
---|---
MD | GNN-MD | Relative Error ($\%$)
95 K | 2.16 | 2.25 |  $\sim$ 4.17
100 K | 2.05 | 2.05 |  $\sim$ 0.0
110 K | 1.83 | 1.83 |  $\sim$ 0.0
Table 2: The correlation time (time taken for $\beta_{\Delta}(q)$ to drop to
0.5) for the 4000 atom system at different temperatures using the MD and GNN-
MD simulation. The $\beta_{\Delta}(q)$ as a function of $\Delta_{\tau}$ for
all temperatures can be found in Fig. 14.

###  3.4 Solid phase and melting

In order to test the model performance in the solid phase, we performed MD and
GNN-MD simulations of a face-centered cubic (FCC) lattice at 20 K, which is
below the melting temperature. The resulting pair distribution function $g(r)$
captures the amplitude of thermal vibration in the solid phase. Figure 8 shows
that the GNN-MD captures the $g(r)$ of the solid phase accurately even for
Model A, which is the GNN-based force field trained exclusively on liquid
configurations. Similar results for the predicted $g(r)$ of the solid phase
are observed using Model B (not shown).

Figure 8: The predicted pair distribution function $g(r)$ of the solid phase
at 20 K using Model A, which is trained only using liquid configurations.

To further test the capability of the model to capture solid phase
configurations, we estimate the melting point of the system by using the
solid-liquid interface method [28]. The solid half is initialized as an FCC
lattice and is equilibrated at 20 K, whereas, the liquid half is equilibrated
at 100 K to yield the initial configuration, as shown in Fig. 9.

Figure 9: (a) The solid-liquid interface used for the melting point estimation
simulation and (b) the interface velocity at different temperatures for Model
A and Model B.

The solid-liquid interface is maintained at 100 K, 61.36 K, and 20 K across
three different simulations. The interface will move into the liquid or solid
phase depending on whether the temperature is below or above the melting
point. To quantify these transitions we estimate the velocity of the solid-
liquid interface. Our results show a linear dependence of the interface
velocity with the simulation temperature, as shown in Fig. 9. Here the melting
point corresponds to the temperature where the interface velocity is zero. We
obtain the melting point from the regular MD simulation to be 55.2 K, whereas
the melting point from the GNN-MD is: 56.4 K from Model A and 55.8 K from
Model B. (The differences here are below the error bar $\sim 1.34$ K.) This
demonstrates that the GNN force field is able to capture the solid-liquid
phase transition with sufficient accuracy. It is somewhat surprising that this
level of agreement can be reached even for Model A, which is trained on liquid
configurations only.

To see whether Model A can capture all the properties of the solid phase
without being trained on any solid configurations, we performed further tests.
First, we examine the distribution of the eigenfrequencies of a perfect FCC
crystal. The eigenfrequencies can be obtained by first diagonalizing the
Hessian matrix, $\mathcal{H}$, and then dividing the eigenvalues by the atomic
mass $m$ and taking the square root. The distribution of the eigenfrequencies
provides an estimate of the phonon density of states (PDOS) in the crystal, as
shown in Fig 10. For this analysis, we use a simulation cell with 2048 atoms,
such that the $\mathcal{H}$ is defined over 6144 degrees of freedom. The
details of the calculations can be found in A.4.

Figure 10: The phonon density of states obtained from Lennard-Jones force
field computation and the GNN force field trained on: (a) exclusively liquid
configurations (Model A), and (b) solid and liquid configurations (Model B).

The Hessian matrix computed using the Lennard-Jones potential
($\mathcal{H}_{\rm LJ}$) is perfectly symmetric whereas the one computed from
the GNN force field ($\tilde{\mathcal{H}}_{\rm GNN}$) is only approximately
symmetrical. Here we use the ∼ sign to indicate that the matrix is not exactly
symmetrical. Before computing the eigenfrequencies we need to symmetrize the
Hessian matrix by $\mathcal{H}_{\rm GNN}^{\rm
s}=\frac{1}{2}\left(\tilde{\mathcal{H}}_{\rm GNN}+\tilde{\mathcal{H}}_{\rm
GNN}^{\rm T}\right)$. Figure 10 shows that the PDOS computed from GNN Model A
(trained on liquid configurations only) exhibits an appreciable difference
from the reference LJ model. The agreement with the reference LJ model becomes
significantly improved for GNN Model B (trained on both liquid and solid
configurations), as shown in Fig. 10.

The lack of perfect symmetry in the Hessian matrix is a type of error in the
GNN-based force field considered here, in which the atomic forces are not
obtained from the partial derivatives of a potential energy function. To
quantify this error, we define a _symmetricity_ measure $\mathcal{S}$ for any
square matrix (see A.4) such that $\mathcal{S}=1$ if the matrix is perfectly
symmetric.

Configuration | Training Data
---|---
Liquid (Model A) | Solid and Liquid (Model B)
Solid | 0.785 | 0.949
Liquid | 0.984 | 0.984
Table 3: The symmetricity, $\mathcal{S}(\tilde{\mathcal{H}}_{\rm GNN})$, of
the Hessian matrix of GNN Model A (trained on liquid configuration only) and
Model B (trained on both liquid and solid configurations) when tested on a
solid configuration (at $\sim 20$ K) and a liquid configuration (at $\sim 100$
K).

Table 3 shows the symmetricity of the Hessian matrix from Model A and Model B
evaluated on a solid (perfect FCC crystal) and on a randomly chosen liquid
configuration. The symmetricity on the liquid configuration is 0.984 for both
Model A and Model B, which is considered acceptable in this study. The
symmetricity on the solid configuration is only 0.785 for Model A, which is
correlated with its poor prediction of the PDOS of the perfect crystal (Fig.
10). On the other hand, after the solid configurations are included in the
training set, Model B gives a symmetricity of 0.949 for the solid
configuration, which is correlated with the improved prediction of the PDOS of
the perfect crystal (Fig. 10). Therefore, our results demonstrate the need for
a comprehensive set of tests to evaluate the transferability of machine-
learned force fields, and the importance of including both liquid and solid
configurations in the training dataset.

##  4 Conclusions

This paper discusses the comprehensive lists of tests needed to check the
transferability of a neural network force field. We carry out our benchmarking
tests on the GNN-based force field presented in [21]. Our analysis is based on
a reference system of liquid Argon by modeling it as a Lennard-Jones fluid.
The typical tests for a force field include validating the radial distribution
function and mean-squared displacement. However, in this paper, we present
additional tests to validate the dynamics from the GNN force field by carrying
out the computational XPCS analysis. These tests are carried out at different
system sizes and temperatures to further establish the transferability of the
GNN force field.

In addition to testing the liquid behavior, we present a few solid-phase tests
such as melting point estimation, which are not usually validated for machine-
learned force fields. Notably, most solid and liquid phase tests were passed
by the force field trained on just liquid configurations, however, estimating
the phonon density of states required solid configurations in the training
dataset. This form of dataset engineering is necessary for the force field to
be able to capture both solid and liquid behavior. In summary, these
benchmarking tests provide a more comprehensive check on the transferability
of the GNN-based force fields.

On having established the transferability of the GNN force field, efforts can
now be made to train the model from AIMD calculations for systems that are not
well modeled by traditional MD force fields. This resulting GNN-AIMD promises
the force field accuracy of AIMD and can be used for larger atomic systems
which is not possible using traditional AIMD (due to computational expense)
and traditional MD (due to the inaccuracy in the force field).

## Conflict of Interest

The authors declare no competing interests.

## Data Access

All data is available in the main text and the Supplementary appendices.
Further information about the computation can be obtained on request from the
corresponding author. The library for the testbed of machine-learned force
fields can be found at TB-MLFF. The code on the XPCS and XSVS analysis is
obtained from the C-XPCS library.

## Acknowledgements

S.M. and W.C. acknowledge support from the Precourt Pioneering Project of
Stanford University. S.Y. was supported by Korea Institute for Advancement of
Technology (KIAT) grant funded by the Korea Government (MOTIE) (P0017304,
Human Resource Development Program for Industrial Innovation). K.K.
acknowledges support from the National Research Foundation of Korea (NRF)
grant funded by the Korean government (MSIT) (NRF- 2022R1A2C2011266).

## References

  * Agrawal and Choudhary [2019] Ankit Agrawal and Alok Choudhary.  Deep materials informatics: Applications of deep learning in materials science.  _MRS Communications_ , 9(3):779–792, 2019.
  * Choudhary et al. [2022] Kamal Choudhary, Brian DeCost, Chi Chen, Anubhav Jain, Francesca Tavazza, Ryan Cohn, Cheol Woo Park, Alok Choudhary, Ankit Agrawal, Simon JL Billinge, et al.  Recent advances and applications of deep learning methods in materials science.  _npj Computational Materials_ , 8(1):1–26, 2022\.
  * Huang and Ling [2020] Liyuan Huang and Chen Ling.  Practicing deep learning in materials science: An evaluation for predicting the formation energies.  _Journal of Applied Physics_ , 128(12):124901, 2020.
  * Ethier et al. [2021] Jeffrey G Ethier, Rohan K Casukhela, Joshua J Latimer, Matthew D Jacobsen, Allen B Shantz, and Richard A Vaia.  Deep learning of binary solution phase behavior of polystyrene.  _ACS Macro Letters_ , 10(6):749–754, 2021.
  * Chandrasekaran et al. [2020] Anand Chandrasekaran, Chiho Kim, Shruti Venkatram, and Rampi Ramprasad.  A deep learning solvent-selection paradigm powered by a massive solvent/nonsolvent database for polymers.  _Macromolecules_ , 53(12):4764–4769, 2020.
  * Shayeganfar and Shahsavari [2021] Farzaneh Shayeganfar and Rouzbeh Shahsavari.  Deep learning method to accelerate discovery of hybrid polymer-graphene composites.  _Scientific Reports_ , 11(1):1–9, 2021.
  * Sino et al. [2021] Nezar Ismet Sino, Rabah N Farhan, and Mohammed E Seno.  Review of deep learning algorithms in computational biochemistry.  In _Journal of Physics: Conference Series_ , volume 1804, page 012135\. IOP Publishing, 2021.
  * Bhatt et al. [2021] Chandradeep Bhatt, Indrajeet Kumar, V Vijayakumar, Kamred Udham Singh, and Abhishek Kumar.  The state of the art of deep learning models in medical science and their challenges.  _Multimedia Systems_ , 27(4):599–613, 2021.
  * Suzuki [2017] Kenji Suzuki.  Overview of deep learning in medical imaging.  _Radiological physics and technology_ , 10(3):257–273, 2017.
  * Shen et al. [2017] Dinggang Shen, Guorong Wu, and Heung-Il Suk.  Deep learning in medical image analysis.  _Annual review of biomedical engineering_ , 19:221–248, 2017.
  * Wang et al. [2018] Han Wang, Linfeng Zhang, Jiequn Han, and E Weinan.  Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics.  _Computer Physics Communications_ , 228:178–184, 2018.
  * Zheng et al. [2021] Tianze Zheng, Weihao Gao, and Chong Wang.  Learning large-time-step molecular dynamics with graph neural networks.  _arXiv preprint arXiv:2111.15176_ , 2021.
  * Do et al. [2022] Hung N Do, Jinan Wang, Apurba Bhattarai, and Yinglong Miao.  Glow: A workflow integrating gaussian-accelerated molecular dynamics and deep learning for free energy profiling.  _Journal of Chemical Theory and Computation_ , 18(3):1423–1436, 2022.
  * Perez et al. [2009] Danny Perez, Blas P Uberuaga, Yunsic Shim, Jacques G Amar, and Arthur F Voter.  Accelerated molecular dynamics methods: introduction and recent developments.  _Annual Reports in computational chemistry_ , 5:79–98, 2009\.
  * Schütt et al. [2018] Kristof T Schütt, Huziel E Sauceda, P-J Kindermans, Alexandre Tkatchenko, and K-R Müller.  Schnet–a deep learning architecture for molecules and materials.  _The Journal of Chemical Physics_ , 148(24):241722, 2018.
  * Marcolongo et al. [2020] Aris Marcolongo, Tobias Binninger, Federico Zipoli, and Teodoro Laino.  Simulating diffusion properties of solid-state electrolytes via a neural network potential: performance and training scheme.  _ChemSystemsChem_ , 2(3):e1900031, 2020.
  * Balyakin and Rempel [2020] IA Balyakin and AA Rempel.  Machine learning interatomic potential for molten tizrhfnb.  In _AIP Conference Proceedings_ , volume 2313, page 030037. AIP Publishing LLC, 2020.
  * Smith et al. [2020] Justin S Smith, Roman Zubatyuk, Benjamin Nebgen, Nicholas Lubbers, Kipton Barros, Adrian E Roitberg, Olexandr Isayev, and Sergei Tretiak.  The ani-1ccx and ani-1x data sets, coupled-cluster and density functional theory properties for molecules.  _Scientific data_ , 7(1):1–10, 2020.
  * Dajnowicz et al. [2022] Steven Dajnowicz, Garvit Agarwal, James M Stevenson, Leif D Jacobson, Farhad Ramezanghorbani, Karl Leswing, Richard A Friesner, Mathew D Halls, and Robert Abel.  High-dimensional neural network potential for liquid electrolyte simulations.  _The Journal of Physical Chemistry B_ , 2022.
  * Mohanty et al. [2023] Shaswat Mohanty, James Stevenson, Andrea Browning, Leif Jacobson, Karl Leswing, Mathew Halls, and Mohammad Atif Faiz Afzal.  Development of scalable and generalizable machine learned force field for polymers.  2023\.  URL https://doi.org/10.26434/chemrxiv-2023-4z2bw.
  * Li et al. [2022] Zijie Li, Kazem Meidani, Prakarsh Yadav, and Amir Barati Farimani.  Graph neural networks accelerated molecular dynamics.  _The Journal of Chemical Physics_ , 156(14):144103, 2022.
  * Mohanty et al. [2022] Shaswat Mohanty, Christopher B. Cooper, Hui Wang, Mengning Liang, and Wei Cai.  Computational approaches to model x-ray photon correlation spectroscopy from molecular dynamics.  _Modelling and Simulation in Materials Science and Engineering_ , 8 2022.
  * Thompson et al. [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  _Comp. Phys. Comm._ , 271:108171, 2022.
  * Martyna and Tuckerman [1999] Glenn J Martyna and Mark E Tuckerman.  A reciprocal space based method for treating long range interactions in ab initio and force-field-based calculations in clusters.  _The Journal of chemical physics_ , 110(6):2810–2821, 1999.
  * Larsen et al. [2017] Ask Hjorth Larsen, Jens JØrgen Mortensen, Jakob Blomqvist, Ivano E. Castelli, Rune Christensen, Marcin Dułak, Jesper Friis, Michael N. Groves, BjØrk Hammer, Cory Hargus, Eric D. Hermes, Paul C. Jennings, Peter Bjerre Jensen, James Kermode, John R. Kitchin, Esben Leonhard Kolsbjerg, Joseph Kubal, Kristen Kaasbjerg, Steen Lysgaard, Jón Bergmann Maronsson, Tristan Maxson, Thomas Olsen, Lars Pastewka, Andrew Peterson, Carsten Rostgaard, Jakob SchiØtz, Ole Schütt, Mikkel Strange, Kristian S. Thygesen, Tejs Vegge, Lasse Vilhelmsen, Michael Walter, Zhenhua Zeng, and Karsten W. Jacobsen.  The atomic simulation environment—a python library for working with atoms.  _Journal of Physics: Condensed Matter_ , 29:273002, 6 2017\.
  * Eastman et al. [2017] Peter Eastman, Jason Swails, John D. Chodera, Robert T. McGibbon, Yutong Zhao, Kyle A. Beauchamp, Lee Ping Wang, Andrew C. Simmonett, Matthew P. Harrigan, Chaya D. Stern, Rafal P. Wiewiora, Bernard R. Brooks, and Vijay S. Pande.  Openmm 7: Rapid development of high performance algorithms for molecular dynamics.  _PLOS Computational Biology_ , 13:e1005659, 7 2017.
  * Hendrycks and Gimpel [2016] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus), 2016.  URL https://arxiv.org/abs/1606.08415.
  * Zhu et al. [2021] Li-Fang Zhu, Jan Janssen, Shoji Ishibashi, Fritz Körmann, Blazej Grabowski, and Jörg Neugebauer.  A fully automated approach to calculate the melting temperature of elemental crystals.  _Computational Materials Science_ , 187:110065, 2021.
  * Prince and Wilson [2004] Edward Prince and Arthur James Cochran Wilson.  _International tables for crystallography_ , volume 100.  Kluwer., 2004.
  * Chapman et al. [2006] Henry N Chapman, Anton Barty, Michael J Bogan, Sébastien Boutet, Matthias Frank, Stefan P Hau-Riege, Stefano Marchesini, Bruce W Woods, Saša Bajt, W Henry Benner, et al.  Femtosecond diffractive imaging with a soft-x-ray free-electron laser.  _Nature Physics_ , 2(12):839–843, 2006.
  * Jakeman [1974] E Jakeman.  Photon correlation.  In _Photon correlation and light beating spectroscopy_ , pages 75–149. Springer, 1974.

##  Appendix A Analysis Details

###  A.1 Mean-squared Displacement

Mean squared displacement (MSD) is a metric used in statistical mechanics to
measure the diffusive motion of the material system being simulated. The MSD
is the most commonly used metric to measure the spatial deviation as a
function of time. The MSD is an ensemble average and is computed using the
position ($\bm{r}\in\mathcal{R}^{d}$ where $d$ is the dimensionality of the
position),

${\rm
MSD}(\tau)=\langle(\bm{r}(\tau)-\bm{r}(0))^{2}\rangle=\frac{1}{N}\sum_{i=1}^{N}(\bm{r}^{i}(\tau)-\bm{r}^{i}(0))^{2}.$
(6)

For a purely diffusive motion, such as the dynamics of Brownian motion, the
MSD is related to the diffusivity, $D$, by,

${\rm MSD}(\tau)=2dD\tau.$ (7)

For the liquid simulations, the MSD varies linearly with time after the
initial slow dynamics (attributed to the sub-diffusive motion brought about by
the caging effect). As a consequence, the MSD is analyzed for $t>1$ ps.

###  A.2 Basics of structural analysis

An important structural measure of a collection of atoms is the pair
distribution function, $g(\bm{r},t)$, defined by,

$g(\bm{r},t)=\frac{1}{N\rho_{0}}\sum_{i}\sum_{\begin{subarray}{c}j\\\ j\neq
i\end{subarray}}\delta(\bm{r}-(\bm{r}_{i}(t)-\bm{r}_{j}(t))),$ (8)

where $\rho_{0}$ is the bulk atomic density of the sample, $N$ is the total
number atoms, and $\bm{r}_{i}$ (and $\bm{r}_{j}$) is the atomic position in
the simulation box. Although the convention is to exclude the correlation
between the atom and itself (i.e. $i=j$) from the definition of $g(\bm{r},t)$,
in the following it is more convenient to introduce an alternative definition,
$\tilde{g}(\bm{r},t)$, where such a constraint is removed.

$\tilde{g}(\bm{r},t)=\frac{1}{N\rho_{0}}\sum_{i}\sum_{j}\delta(\bm{r}-(\bm{r}_{i}(t)-\bm{r}_{j}(t)))=g(\bm{r},t)+\frac{1}{\rho_{0}}\,\delta(\bm{r})$
(9)

When all atoms are of the same type, the Fourier transform of
$\rho_{0}\,\tilde{g}(\bm{r},t)$ is the structure factor,
$\tilde{S}(\bm{q},t)$,

$\tilde{S}(\bm{q},t)=\frac{1}{N}\sum_{i}\sum_{j}{\rm e}^{-{\rm
i}\bm{q}\cdot[\bm{r}_{i}(t)-\bm{r}_{j}(t)]},$ (10)

where $\bm{q}$ is the scattered wave-vector. When the sample contains atoms of
different types, the definition of the structure factor is generalized to the
following

$S(\bm{q},t)=\frac{1}{\sum_{j}f_{j}(\bm{q})^{2}}\sum_{i}\sum_{j}f_{i}(\bm{q})f_{j}(\bm{q})\,{\rm
e}^{-{\rm i}\bm{q}\cdot[\bm{r}_{i}(t)-\bm{r}_{j}(t)]}.$ (11)

Here $f_{i}(\bm{q})$ is the X-ray atomic form factor of atom $i$. The X-ray
atomic form factor can be obtained by the Fourier transform of the electron
density field for a given type of atom. For many elements, the atomic form
factor can be well parameterized by a sum of Gaussians [29].

###  A.3 Basics of XPCS and scattering statistics

We then consider the time correlation function of the intensity of the XPCS
speckles that we observe. The normalized auto-correlation function helps in
studying the dynamics of the simulation. The auto-correlation is denoted by
$g_{2}(\tau)$, since it is a second order correlation and is given by,

$\displaystyle g_{2}(\bm{q},\tau)$ $\displaystyle=\langle
I(\bm{q},t)I(\bm{q},t+\tau)\rangle/\langle I\rangle^{2}$ (12) $\displaystyle
I(\bm{q},\tau)$ $\displaystyle\propto S(\bm{q},\tau)$ (13)

Another metric that helps us in studying the dynamic characteristics for the
simulation is the intermediate scattering function, $f(\bm{q},\tau)$,
described as,

$f(\bm{q},\tau)=\frac{1}{N}\left\langle\sum_{i\neq j}e^{({\rm
i}\bm{q}(r_{i}(0)-r_{j}(\tau))}\right\rangle_{t}$ (14)

For a simple diffusive process following Brownian motion, the intermediate
scattering function reduces to a single exponential,
$F(\bm{q},\tau)=\frac{f(\bm{q},\tau)}{f(\bm{q},0)}=e^{-\Gamma\tau}$, where
$\Gamma=D\bm{q}^{2}$ [30].

Since the time auto-correlation ($g_{2}(\bm{q},\tau)$) is related to the
$f(\bm{q},\tau)$ by the Siegert relation [31], it can also be represented as a
single exponential under the assumption of Brownian dynamics,

$\displaystyle g_{2}(\bm{q},\tau)$ $\displaystyle=1+\beta(\bm{q})\lvert
F(\bm{q},\tau)\rvert^{2}$
$\displaystyle=1+\beta(\bm{q})\exp[-2\Gamma(\bm{q})\tau]$ (15) $\displaystyle
g_{2}(\bm{q},\tau)-1$
$\displaystyle=G(\bm{q},\tau)\propto\exp[-2\Gamma(\bm{q})\tau]$ (16)

where $\Gamma(\bm{q})$ is the relaxation rate and the factor of 2 arises from
the Siegert relation which relates the intensity auto-correlation function to
the electric field correlation function [31].

Furthermore, the auto-correlation of the speckle at $\bm{q}$ is denoted by
$g_{2}(\bm{q},\tau)$ and it is non-dimensionalized by the square of its mean
at $\tau=0$, as given below for infinitesimally short pulses,

$g_{2}(\bm{q},\tau)=\frac{\langle
I(\bm{q},t)\,I(\bm{q},t+\tau)\rangle_{t}}{\langle
I(\bm{q},t)\rangle_{t}^{2}}=\frac{\langle
S(\bm{q},t)\,S(\bm{q},t+\tau)\rangle_{t}}{\langle
S(\bm{q},t)\rangle_{t}^{2}},$ (17)

where $\langle\cdot\rangle_{t}$ represents the time-averaged value of the
enclosed entity. Experimentally, the X-ray pulses always have a finite
duration. Hence the intensity should be replaced by the time-averaged X-ray
intensity, $I_{\Delta}(\bm{q})$, over the exposure duration of
$\Delta_{\tau}$,

$I_{\Delta}(\bm{q})=\int_{t}^{t+\Delta_{\tau}}I(\bm{q},t)\,{\rm d}t.$ (18)

In this case,

$g_{2}(\bm{q},\tau)=\frac{\langle
I_{\Delta}(\bm{q},t)\,I_{\Delta}(\bm{q},t+\tau)\rangle_{t}}{\langle
I_{\Delta}(\bm{q},t)\rangle_{t}^{2}}.$ (19)

For a X-ray speckle pattern obtained from a single X-ray pulse, the optical
contrast $\beta({q})$ is defined by the variance of the intensity divided by
the square of its mean [22]. While $\beta(q)$ is given by the scattering
intensity distribution $I(\bm{q})$ from an infinitesimally short pulse, we
denote the optical contrast from an X-ray pulse of finite duration
${\Delta}_{t}$ as $\beta_{\Delta}(q)$ where,

$\beta_{\Delta}(q)=\frac{\langle I_{\Delta}(\bm{q})^{2}\rangle_{q}-\langle
I_{\Delta}(\bm{q})\rangle_{q}^{2}}{\langle
I_{\Delta}(\bm{q})\rangle_{q}^{2}},$ (20)

where $\langle\cdot\rangle_{q}$ represents the average over all detector
pixels that satisfy $q-dq/2\leq\lvert\bm{q}\rvert<q+dq/2$, for a small $dq$.
For an ergodic and isotropic system, the distribution of pixel intensity at
around a particular $\bm{q}$ is the same in $\bm{q}$-space and $t$, so that

$\langle I_{\Delta}(\bm{q})\rangle_{q}\approx\langle
I_{\Delta}(\bm{q})\rangle_{t},\quad\langle
I_{\Delta}(\bm{q})^{2}\rangle_{q}\approx\langle
I_{\Delta}(\bm{q})^{2}\rangle_{t}.$ (21)

This means that we can also define an optical contrast $\beta_{0}(\bm{q})$
from the time variation of the intensity at a single $\bm{q}$, i.e,

$\beta_{0}(\bm{q})=g_{2}(\bm{q},\tau=0)-1.$ (22)

For an ergodic and isotropic system, $\beta_{\Delta}(q)$ and
$\beta_{0}(\bm{q})$ should be equal. The detailed discussion on the XPCS
relations and the statistics of X-ray speckles, the computational algorithm
and its implementation can be found in [22].

###  A.4 Hessian matrix and phonon density of states

To estimate the phonon density of states (PDOS) we examine the eigenvalues of
the Hessian matrix, $\mathcal{H}$, of our relaxed system. The $\mathcal{H}$ is
defined as,

$\displaystyle\mathcal{H}_{ij}$ $\displaystyle=\frac{\partial}{\partial
x_{j}}\left(\frac{\partial U}{\partial x_{i}}\right),$ (23)
$\displaystyle=-\frac{\partial{f}_{i}}{\partial x_{j}},$ (24)

where $U$ is the potential energy and ${f}_{i}$ is the atomic force component
corresponding to the $i^{\rm th}$ degree of freedom ($i$ goes from 1 to $3N$).
We can approximate this derivative numerically by a central difference scheme,

$\mathcal{H}_{ij}=\frac{{f}_{i}(x_{j}-\Delta x)-{f}_{i}(x_{j}+\Delta
x)}{2\Delta x},$ (25)

where $\Delta x=3.405\times 10^{-6}$ Å($10^{-6}$ in LJ unit). If $\lambda_{n}$
represents the eigenvalues of $\mathcal{H}$ then the eigenfrequencies
$\nu_{n}$ can be obtained from

$\nu_{n}=\frac{1}{2\pi}\sqrt{\frac{\lambda_{n}}{m}}.$ (26)

The distribution of the $\nu_{n}$ gives us the PDOS in arbitrary units. The
Hessian matrix $\mathcal{H}$ is expected to be symmetric for a conservative
force field that is derived from an interatomic potential. Consistency or
symmetricity requires,

$\frac{\partial{f}_{i}}{\partial x_{j}}=\frac{\partial{f}_{j}}{\partial
x_{i}}.$ (27)

The Hessian matrix computed using the Lennard-Jones potential
($\mathcal{H}_{\rm LJ}$) is perfectly symmetric whereas the one computed from
the GNN force field ($\tilde{\mathcal{H}}_{\rm GNN}$) is only nearly
symmetric. We define a measure for symmetricity, $\mathcal{S}$, for a matrix
$A$ as

$\mathcal{S}(A)=\frac{||A_{\rm s}||_{2}-||A_{\rm a}||_{2}}{||A_{\rm
s}||_{2}+||A_{\rm a}||_{2}},$ (28)

where,

$\displaystyle A_{\rm s}$ $\displaystyle=\frac{A+A^{T}}{2},$ (29)
$\displaystyle A_{\rm a}$ $\displaystyle=\frac{A-A^{T}}{2}.$ (30)

The symmetricity is such that $-1\leq\mathcal{S}(A)\leq 1$, where
$\mathcal{S}(A)=1$ for a perfectly symmetric matrix and $\mathcal{S}(A)=-1$
for a perfectly anti-symmetric matrix.

##  Appendix B Analysis at different temperatures

In this appendix, we show the results corroborating our findings that have
been presented in main text. We present simulation results corresponding to
the 4000 atoms case which is run at 95 K, 100 K, and 110 K.

Figure 11: The mean-squared displacement (MSD) for the (a) 256 atoms at 100 K
and the 4000 atoms at (b) 100 K, (c) 95 K, and (d) 110 K, MD and GNN-MD
simulation over the 539 ps trajectory.

We see that the MSD is in close agreement between the MD and GNN-MD
simulations, as shown in Fig. 11, irrespective of the temperature at which the
simulations are carried out.

Figure 12: The $g_{2}(q,\tau)$ at $q=1.844\pm 0.029$ Å-1 over the first $\sim
3.56$ ps for the (a) 256 atoms (25 tracked atoms) at 100 K and the 4000 atoms
(45 tracked atoms) at (b) 100 K, (c) 95 K and (d) 110 K, MD and GNN-MD
simulation.

We observe that the decay in $g_{2}(q,\tau)$ between the MD and the GNN-MD
simulation is also in agreement at all three temperatures.

Figure 13: The $\Gamma(q)$ as a function of $q^{2}$ for the (a) 256 atoms at
100 K and the 4000 atoms at (b) 100 K, (c) 95 K, and (d) 110 K, MD and GNN-MD
simulation.

We examine the $\Gamma(q)$ as a function of $q^{2}$ for different temperatures
and see a very close agreement between the MD and GNN-MD simulations.

Figure 14: The $\beta_{\Delta}(q)$ as a function of $\Delta_{\tau}$ for the
(a) 256 atoms at 100 K and the 4000 atoms at (b) 100 K, (c) 95 K and (d) 110
K, MD and GNN-MD simulation for $\lvert\bm{q}\rvert=1.844\pm 0.029$ Å-1.

In addition to the previous analysis, we also examine the XSVS analysis by
plotting the optical contrast, $\beta_{\Delta}(q)$, as a function of the
exposure time $\Delta_{\tau}$. The decay in the optical contrast is also in
agreement irrespective of the temperature of the simulation which shows that
the GNN force field approximates the Lennard-Jones potential sufficiently to
capture the dynamics on liquid phase.