††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Correcting force error-induced underestimation of lattice thermal
conductivity in machine learning molecular dynamics

Xiguang Wu  Guangzhou Key Laboratory of Low-Dimensional Materials and Energy
Storage Devices, School of Materials and Energy, Guangdong University of
Technology, Guangzhou 510006, China  Wenjiang Zhou  Department of Energy and
Resources Engineering, Peking University, Beijing 100871, China  School of
Advanced Engineering, Great Bay University, Dongguan 523000, China  Haikuang
Dong  College of Physical Science and Technology, Bohai University, Jinzhou
121013, China  Penghua Ying  Department of Physical Chemistry, School of
Chemistry, Tel Aviv University, Tel Aviv, 6997801, Israel  Yanzhou Wang  MSP
group, QTF Centre of Excellence, Department of Applied Physics, Aalto
University, FI-00076 Aalto, Espoo, Finland  Bai Song<EMAIL_ADDRESS>Department of Energy and Resources Engineering, Peking University, Beijing
100871, China  Department of Advanced Manufacturing and Robotics, Peking
University, Beijing 100871, China  National Key Laboratory of Advanced
MicroNanoManufacture Technology, Beijing, 100871, China  Zheyong Fan
<EMAIL_ADDRESS>College of Physical Science and Technology, Bohai
University, Jinzhou 121013, China  Shiyun Xiong<EMAIL_ADDRESS>Guangzhou
Key Laboratory of Low-Dimensional Materials and Energy Storage Devices, School
of Materials and Energy, Guangdong University of Technology, Guangzhou 510006,
China

###### Abstract

Machine learned potentials (MLPs) have been widely employed in molecular
dynamics (MD) simulations to study thermal transport. However, literature
results indicate that MLPs generally underestimate the lattice thermal
conductivity (LTC) of typical solids. Here, we quantitatively analyze this
underestimation in the context of the neuroevolution potential (NEP), which is
a representative MLP that balances efficiency and accuracy. Taking crystalline
silicon, GaAs, graphene, and PbTe as examples, we reveal that the fitting
errors in the machine-learned forces against the reference ones are
responsible for the underestimated LTC as they constitute external
perturbations to the interatomic forces. Since the force errors of a NEP model
and the random forces in the Langevin thermostat both follow a Gaussian
distribution, we propose an approach to correcting the LTC by intentionally
introducing different levels of force noises via the Langevin thermostat and
then extrapolating to the limit of zero force error. Excellent agreement with
experiments is obtained by using this correction for all the prototypical
materials over a wide range of temperatures. Based on spectral analyses, we
find that the LTC underestimation mainly arises from increased phonon
scatterings in the low-frequency region caused by the random force errors.

##  I Introduction

Lattice thermal conductivity (LTC) of solids is a crucial physical property in
many applications including thermal management of electronics [1, 2],
thermoelectric energy conversion [3, 4, 5], and thermal barrier coatings [6,
7]. Predicting and engineering LTC [8] is therefore of broad interest.
Nevertheless, challenges abound owing to the presence of complex structures
[9], defects [10], and disorders [10]. Among various approaches to calculating
LTC [11], molecular dynamics (MD) simulation plays a unique role due to its
versatility and its natural inclusion of the full lattice anharmonicity. MD
simulations are widely applicable in crystals, glasses [12], and also liquids
[13]. Two basic categories are commonly used, including equilibrium molecular
dynamics (EMD) base on the Green-Kubo formalism [14, 15] and non-equilibrium
molecular dynamics (NEMD) based on Fourier’s law of heat conduction. Notably,
the homogeneous non-equilibrium molecular dynamics (HNEMD) method, initially
developed by Evans [16] for pairwise interactions and recently generalized to
many-body interactions [17], offers great efficiency for LTC calculations.
However, the applicability and predictive power of MD simulations have long
been limited by the availability and accuracy of empirical interatomic
potentials.

A promising solution to this issue involves constructing machine learned
potentials (MLPs) trained against reference energies, forces, and virial
stresses of diverse atomic structures calculated at the quantum mechanical
level. Many MLPs have been used for thermal conductivity modeling. Enabled by
MLPs, the LTCs of many crystals with strong phonon anharmonicity or disorder
have been successfully obtained through MD simulation driven by MLP (MLMD),
including e.g., amorphous GeTe [18], SnSe [19], PbTe [20], metal-organic
frameworks [21], and PH4AlBr4 [22]. Moreover, with proper quantum corrections,
quantitative agreement with experimental data has also been achieved for
amorphous materials [23, 24] and liquid water [25] over a wide range of
temperatures. Despite these successes, previous works have also shown that for
materials with relatively high LTCs, such as CoSb3 [26] and cubic silicon
(c-Si) [27], the predicted LTCs from MLMD calculations are generally lower
than the experimental values. To the best of our knowledge, this discrepancy
remains to be systematically understood and corrected, which constitutes the
main focus of our present work.

In light of the critical impact of the interatomic forces on the accuracy of
MD simulations, we first evaluate the effect of random forces on LTC using
HNEMD simulations with a Langevin thermostat [28] based on empirical
potentials. A decrease in LTC with increasing level of random forces is
consistently observed in six representative materials: amorphous silicon
(a-Si), c-Si, cubic germanium (c-Ge), Si-Ge alloy, graphene, and
$(10,10)$-carbon nanotube (CNT). Subsequently, we focus on four benchmark
materials including c-Si, GaAs, graphene, and PbTe, and perform LTC
calculations using MLMD. In particular, we employ the neuroevolution potential
(NEP) [29, 30, 31] for its balanced efficiency and accuracy. Similar to
literature results, we observe a consistent underestimation of LTC from the
MLMD simulations, as compared to the experimental values. However, since the
residual force errors of a NEP model and the random forces (white noises) in
the Langevin thermostat both follow a Gaussian distribution, we propose an
approach to correcting the LTC by intentionally introducing different levels
of force noises via the Langevin thermostat and then extrapolating to the
limit of zero force error. This extrapolation successfully corrects the LTCs,
leading to excellent agreement with experimental data for all the materials
considered in a wide range of temperatures. Spectral analyses reveal that the
LTC underestimation before the correction mainly originates from increased
phonon scatterings at low frequencies caused by the force errors.

##  II Methods

###  II.1 Neuroevolution potential

####  II.1.1 The NEP formalism

In this section, we briefly review the NEP formalism [29, 30, 31]. NEP uses a
feedforward neural network to correlate a local descriptor with the site
energy $U_{i}$ of atom $i$. In a single-hidden-layer neural network comprising
$N_{\rm neu}$ hidden neurons, $U_{i}$ is expressed as:

$U_{i}=\sum_{\mu=1}^{N_{\rm
neu}}\omega_{\mu}^{(1)}\tanh\left(\sum_{\nu=1}^{N_{\rm
des}}\omega_{\mu\nu}^{(0)}q_{\nu}^{i}-b_{\mu}^{(0)}\right)-b^{(1)},$ (1)

where $N_{\rm des}$ is the number of descriptor components, $q_{\nu}^{i}$ is
the $\nu$-th descriptor component of atom $i$, $\omega_{\mu\nu}^{(0)}$,
$\omega_{\mu}^{(1)}$, $b_{\mu}^{(0)}$, and $b^{(1)}$ are the trainable
parameters, and $\tanh(x)$ is the nonlinear activation function in the hidden
layer.

The descriptor vector in NEP includes radial and angular components. The
radial components $q_{n}^{i}$ $(0\leq{n}\leq n_{\rm max}^{\rm R})$ are defined
as

$q_{n}^{i}=\sum_{j\neq{i}}g_{n}(r_{ij}),$ (2)

where $r_{ij}$ is the distance between atoms $i$ and $j$ and $g_{n}(r_{ij})$
are a set of radial functions, each of which is formed by a linear combination
of Chebyshev polynomials. The angular components include the so-called
$n$-body ($n\geq 3$) correlations. For example, the $3$-body ones $q_{nl}^{i}$
$(0\leq{n}\leq n_{\rm max}^{\rm A}$, $1\leq{l}\leq l_{\rm max})$ are defined
as

$q_{nl}^{i}=\frac{2l+1}{4\pi}\sum_{j\neq{i}}\sum_{k\neq{i}}g_{n}(r_{ij})g_{n}(r_{ik})P_{l}(\cos\theta_{ijk}).$
(3)

Here, $P_{l}$ is the Legendre polynomial and $\theta_{ijk}$ is the angle
formed by the $ij$ and $ik$ bonds. Note that the radial functions
$g_{n}(r_{ij})$ for the radial and angular descriptor components can have
different cutoff radii, which are denoted as $r_{\rm c}^{\rm R}$ and $r_{\rm
c}^{\rm A}$, respectively. The free parameters are optimized using the
separable natural evolutionary strategy [32] by minimizing a loss function
that is a weighted sum of the root-mean-square errors (RMSEs) of energy,
force, and virial stress, for $N_{\rm gen}$ generations with a population size
of $N_{\rm pop}$. The hyperparameters used for all the materials considered in
this work are listed in Table S1.

####  II.1.2 Training datasets

For c-Si, GaAs, graphene, and PbTe, we generate datasets through density
functional theory (DFT) calculations using the vasp [33] with the ion-electron
interactions described by the projector-augmented wave method [33, 34]. For
GaAs, the Perdew-Zunger functional with the local density approximation [35]
is used to describe the exchange-correlation of electrons, while the Perdew-
Burke-Ernzerhof functional with the generalized gradient approximation [36] is
used for the other materials. The cutoff energy is 400 eV for PbTe and 600 eV
for the other materials. The k-point mesh is $4\times 4\times 4$ for c-Si,
$2\times 2\times 2$ for GaAs and PbTe, and $4\times 4\times 1$ for graphene.
The energy convergence threshold is $10^{-6}$ eV for c-Si and graphene and
$10^{-8}$ eV for GaAs and PbTe.

The dataset for each materials consists of structures from _ab initio_
molecular dynamics (AIMD) simulations (called AIMD structures below) possibly
supplemented by those from random cell deformations and atom displacements
(called perturbation structures below). For c-Si, there are 900 AIMD
structures sampled at various temperatures (100 K to 1000 K) and strain states
(unstrained, uniaxial strains of $\pm 1\%$ and $\pm 2\%$, biaxial strains of
$\pm 0.5\%$ and $\pm 1\%$) and 70 perturbation structures. Each c-Si structure
has 64 atoms. For GaAs, there are 197 AIMD structures sampled at various
temperatures (100 K to 900 K) in the $NPT$ ensemble and 99 perturbation
structures up to ±4% strains. Each GaAs structure has 250 atoms. For graphene,
there are 700 AIMD structures sampled at various temperatures (100 K to 1000
K) and strain states (unstrained, biaxial strains of $\pm 0.5\%$, $\pm 1\%$,
and $2\%$). Each graphene structure has 72 atoms. For PbTe, there are 60 AIMD
structures sampled from 100 to 1100 K with fixed cell and 64 perturbation
structures up to ±4% strains. Each PbTe structure has 216 atoms.

After obtaining the datasets, we used the gpumd package [37] (the `nep`
executable) to train the NEP models. The parity plots and accuracy metrics are
detailed in Figs. S1-S4. Force test errors will be further discussed and used
in Sec. III.2.

###  II.2 Thermal conductivity calculation using MD

####  II.2.1 The HNEMD method

We use the efficient HNEMD method [17] for many-body potentials to calculate
the LTCs. In HNEMD, an external driving force on each atom $i$

$\mathbf{F}_{i}^{\rm ext}=\mathbf{F}_{\rm e}\cdot\mathbf{W}_{i}$ (4)

is applied during the simulation. Here, $\mathbf{F_{\rm e}}$ is the driving
force parameter (of the dimension of inverse length) and [29, 31]

$\mathbf{W}_{i}=\sum_{j\neq
i}\mathbf{r}_{ij}\otimes\frac{\partial{U_{j}}}{\partial\mathbf{r}_{ji}}$ (5)

is the virial tensor of atom $i$, where $U_{j}$ is the site energy of atom
$j$, $\mathbf{r}_{ij}\equiv\mathbf{r}_{j}-\mathbf{r}_{i}$, $\mathbf{r}_{i}$
being the position of atom $i$. The driving force parameter should be large
enough to ensure a large signal-to-noise ratio and be small enough to maintain
the system in the linear-response regime. In the linear-response regime, the
LTC tensor $\kappa_{\mu\nu}$ can be calculated from the following relation
[17]:

$\frac{\langle J_{\mu}(t)\rangle_{\rm
ne}}{TV}=\sum_{\nu}\kappa_{\mu\nu}{F_{\rm e}}^{\nu},$ (6)

where $\langle J_{\mu}(t)\rangle_{\rm ne}$ represents a non-equilibrium
ensemble average of the heat current, $T$ is the system temperature, and $V$
is the system volume. The heat current for the NEP model has been derived to
be [29, 31]

$\mathbf{J}=\sum_{i}\mathbf{W}_{i}\cdot\mathbf{v}_{i},$ (7)

where $\mathbf{v}_{i}$ is the velocity of atom $i$.

The HNEMD formalism also allows for an efficient calculation of the frequency-
resolved LTC $\kappa(\omega)$ via the following relation [17]:

$\frac{2}{VT}\int_{-\infty}^{+\infty}e^{i\omega
t}K^{\mu}(t)dt=\sum_{\nu}\kappa_{\mu\nu}(\omega){F_{\rm e}}^{\nu},$ (8)

where

$K^{\mu}(t)=\sum_{i}\sum_{\nu}\langle
W_{i}^{\mu\nu}(0)v_{i}^{\nu}(t)\rangle_{\rm ne}$ (9)

is the virial-velocity correlation function.

####  II.2.2 Thermostats in HNEMD simulations

HNEMD simulations are normally performed in the $NVT$ ensemble realized by
using a global thermostat such as the Nosé-Hoover chain (NHC) [38] or the
Bussi-Donadio-Parrinello [39] thermostat. In contrast, a local thermostat such
as the Langevin thermostat [28] is avoided because it can introduce (white)
noises through random forces, leading to the following equations of motion:

$\frac{d\bm{\mathrm{r}}_{i}}{dt}=\frac{\bm{\mathrm{p}}_{i}}{m_{i}},\hskip
28.45274pt\frac{d\bm{\mathrm{p}}_{i}}{dt}=\bm{\mathrm{F}}_{i}-\frac{\bm{\mathrm{p}}_{i}}{\tau_{T}}+\bm{\mathrm{f}}_{i}.$
(10)

Here, $\tau_{T}$ is a time parameter, $\mathbf{r}_{i}$, $\mathbf{p}_{i}$,
$m_{i}$ are respectively the position, momentum, and mass of atom $i$,
$\mathbf{F}_{i}$ is the force on atom $i$ resulting from the interatomic
potential, and $\mathbf{f}_{i}$ is the random force on atom $i$. Each
component of the random force forms a Gaussian distribution with zero mean and
a variance of

$\sigma^{2}_{\rm L}=\frac{2k_{\rm B}Tm}{\tau_{T}\Delta t},$ (11)

where $m$ is the average atom mass in the system, $k_{\rm B}$ is the Boltzmann
constant and $\Delta t$ is the integration time step. The random forces can
affect the dynamics of the system and thus time-correlation properties such as
the heat current autocorrelation function, leading to reduced LTC as compared
to the case of using a global thermostat. Clearly, a smaller $\tau_{T}$ gives
a larger random force variance and a stronger reduction of the LTC. We will
demonstrate this effect using examples.

####  II.2.3 MD simulation details

All the MD simulations are performed using the gpumd package [37] (the `gpumd`
executable), with a time steps of 1 fs. For all the materials, we use a
sufficiently large simulation cell to eliminate finite-size effects. In MD
simulations with empirical potentials, the simulation cells contain $32\,768$
atoms for Si-Ge alloy, c-Ge, c-Si, and a-Si, $15\,416$ atoms for graphene, and
$16\,280$ atoms for the $(10,10)$-CNT. The a-Si samples are prepared by
employing a melt-quench-anneal process, first equilibrating at 2000 K for 10
ns, then quenching down to 300 K during 30 ns, and finally annealing at 300 K
for 10 ns. In MD simulations with NEP models, the simulation cells contain
$13\,824$, $8000$, $16\,000$, and $36\,000$ atoms for c-Si, GaAs, graphene,
and PbTe, respectively. For each material, we first equilibrate the system in
the $NPT$ ensemble (with a target pressure of zero) for 2 ns and $NVT$
ensemble for another 2 ns, and then calculate the LTC in the $NVT$ ensemble
during a production time of 10 to 20 ns. For each material at each
temperature, three to five independent runs are performed to improve the
statistical accuracy and obtain an error estimate.

##  III Results and discussion

###  III.1 Thermal conductivity underestimation in MLMD

Figure 1: Comparison of $\kappa$ for c-Si from NEP-MD simulations and
experimental measurements [40, 41, 42]. Here, the NHC thermostat is used.
Error bars are smaller than the symbol sizes for the calculated values.

To begin with, we take c-Si as an example to demonstrate the thermal
conductivity underestimation in MLMD simulations, using NEP as a
representative MLP. As illustrated in Fig. 1, the calculated LTC values from
300 K to 700 K are consistently lower than experimental measurements,
especially at low temperatures. For instance, at 300 K, MLMD simulations yield
a LTC of $102\pm 1$ W/m-K. While this is more accurate than the value of 240
W/m-K as obtained from a Stillinger-Weber potential [43], it is still
approximately 32% lower than the experimental value of about 150 W/m-K [40].
An EMD simulation based on the Gaussian approximation potential also reported
a lower-than-experiment value of 121 W/m-K [27]. Further, a similar trend of
underestimation is observed at other temperatures for c-Si [27] and in the
case of CoSb3 [26].

###  III.2 Role of force noises in reducing LTC

Table 1: RMSEs $\sigma_{\rm mlp}$ of force prediction for the four NEP models at various temperatures. $T$ (K) |  $\sigma_{\rm mlp}$ (meV/Å)
---|---
c-Si | GaAs | graphene | PbTe
300 | 16.7 | 16.6 | 29.2 | 27.0
400 | 21.3 | 19.8 | 30.1 | 29.9
500 | 28.3 | 23.5 | 32.5 | 34.4
600 | 30.1 | 26.8 | 36.6 | 37.1
700 | 41.6 | 30.2 | 42.4 | 42.1
Figure 2: Force error distribution for c-Si at _T_ = 300 K. Fitting to the
Gaussian distribution yields a coefficient of determination $R^{2}=0.995$.
Figure 3: Inverse LTC ($1/\kappa$) as a function of the random force variance
$\sigma_{\rm L}$ of the Langevin thermostat (see Eq. 11) for (a) a-Si and Si-
Ge alloy, (b) c-Ge and c-Si, and (c) graphene and $(10,10)$-CNT at 300 K. The
hollow and filled symbols are the results from the NHC and Langevin
thermostats, respectively. The solid lines represent linear fits to the
Langevin data only.  Figure 4: Inverse LTC ($1/\kappa$) from NEP-MD
simulations as functions of the total force error $\sigma_{\rm total}$ at
different temperatures for (a) c-Si, (b) GaAs, (c) graphene, and (d) PbTe.
Solid lines indicate linear fits and the points of intersection at
$\sigma_{\rm total}=0$ correspond to the corrected LTC values.

To understand the underestimation of the LTC from MLMD simulations, we notice
that a MLP usually has a certain level of error for force prediction compared
to the reference data. The RMSEs $\sigma_{\rm mlp}$ of force prediction for
the four materials we considered at different temperatures are presented in
Table 1.

A crucial observation is that the force errors follow a Gaussian distribution,
as shown in Fig. 2 for the example of c-Si at 300 K. This distribution is the
same as that for the random forces in the Langevin thermostat, i.e., a
Gaussian distribution with zero mean and a certain variance. Based on this
similarity, an understanding of the underestimation of the LTC by MLMD
simulations can thus be obtained by studying the effect of the Langevin
thermostat on the LTC. One could directly use a NEP model, but due to the
lower computational cost of empirical potentials, we first use the Tersoff
empirical potential [44] to study this effect.

Figure 3 shows the inverse LTC ($1/\kappa$) at 300 K as a function of
$\sigma_{\rm L}$ for six representative materials, including a-Si, Si-Ge
alloy, c-Ge, c-Si, graphene, and $(10,10)$-CNT. As expected, $1/\kappa$
increases with increasing $\sigma_{\rm L}$, which indicates a stronger effect
of the random forces in the Langevin thermostat in reducing the calculated
LTC. Notably, for all the six materials, $1/\kappa$ exhibits a linear
relationship with $\sigma_{\rm L}$. This suggests that the intrinsic LTC
without the influence of the random forces in the Langevin thermostat can be
obtained by extrapolating to $\sigma_{\rm L}=0$. Indeed, the extrapolated
values align well with the results from HNEMD simulations based on the NHC
thermostat that does not involve random forces, with the largest relative
error being $<1.5\%$ (see Table S2).

The linear relation between $1/\kappa$ and $\sigma_{\rm L}$ can be justified
based on the kinetic theory of phonons and Matthiessen’s rule. Taking the
random forces in the Langevin thermostat as an extra source of phonon
scattering, we have

$\frac{1}{\kappa}=\frac{1}{\kappa_{0}}+\frac{1}{1/3Cv_{\rm g}\Lambda_{\rm
L}},$ (12)

where $\kappa$ and $\kappa_{0}$ are the LTCs with and without the influence of
the random forces, respectively, $C$ is the heat capacity, $v_{\rm g}$ is the
phonon group velocity, and $\Lambda_{\rm L}$ is the phonon mean free path
resulting from the random forces in the Langevin thermostat. Under first-order
approximation with sufficiently small $\sigma_{\rm L}$, $1/\Lambda_{\rm L}$
should be proportional to $\sigma_{\rm L}$, which brings Eq. 12 to

$\frac{1}{\kappa}=\frac{1}{\kappa_{0}}+\beta\sigma_{\rm L},$ (13)

which gives the observed linear relation between $1/\kappa$ and $\sigma_{\rm
L}$ with $\beta$ being a slope parameter.

Figure 5:  Corrected (using the Langevin thermostat and exptrapolation) and
uncorrected (using the NHC thermostat) $\kappa$ as a function of temperature
for (a) c-Si, (b) GaAs, (c) graphene, and (d) PbTe. Experimental values are
from Ref. [40, 41, 42] (c-Si), Ref. [45, 46, 47, 48, 49] (GaAs), Ref. [50, 51,
52, 53] (graphene), and Ref. [54, 55] (PbTe).

###  III.3 Correction of LTC in MLMD

Based on the results above, we can understand why MLMD usually underestimates
the LTC, particularly for high-$\kappa$ materials. According to the linear
relation between the inverse LTC and the random force variance, we can devise
a method to correct the underestimation of LTC due to the force errors in
MLMD. To this end, we note that both the force errors in MLMD and the random
forces in the Langevin thermostat follow a Gaussian distribution, and when
they are present simultaneously, a new set of force errors are created with a
larger variance given by

${\sigma_{\rm total}}^{2}={\sigma_{\rm L}}^{2}+{\sigma_{\rm mlp}}^{2},$ (14)

according to the properties of Gaussian distribution. Therefore, we can
intentionally introduce extra force errors by using MLP-based HNEMD
simulations with the Langevin thermostat. The LTC $\kappa_{0}$ without any
force errors (including the force errors of the MLP) can be obtained by an
extrapolation based on the following relation:

$\frac{1}{\kappa}=\frac{1}{\kappa_{0}}+\beta\sigma_{\rm total},$ (15)

where $\kappa$ is the LTC of a material calculated by using MLMD with a
certain force error variance $\sigma_{\rm mlp}$ and the Lagevin thermostat
with a certain random force variance $\sigma_{\rm L}$. The linear relation
between $1/\kappa$ and $\sigma_{\rm total}$ is unfailingly confirmed in Fig. 4
for the four representative materials we considered in a wide range of
temperatures, whose LTCs span three orders of magnitude.

In Fig. 5, we compare the uncorrected and corrected LTCs from MLMD simulations
with experimental results for c-Si, GaAs, graphene, and PbTe. In all the
systems, the uncorrected LTCs are consistently lower than the experimental
results in the entire temperature range due to the presence of force errors in
the MLPs. Remarkably, once the force errors in the MLPs are eliminated via our
extrapolation scheme, the LTCs closely approach the experimental data at all
the temperatures studied. For graphene, the corrected LTCs slightly exceed the
measured values but remain within the experimental uncertainties. This minor
discrepancy could arise from factors such as isotope scattering and finite-
size effects in the experimental setups [50, 51, 52, 53], which generally lead
to reduced LTCs.

Figure 6: Spectral LTC $\kappa(\omega)$ from NEP-MD simulations using the NHC
and Langevin thermostats for (a) c-Si and (b) PbTe, both at 300 K.

The need for LTC correction is more pronounced in materials with higher LTCs
and at lower temperatures, highlighting the more significant impact of force
errors when anharmonic phonon scattering is weak. This explains why the amount
of correction is large for graphene that is one of the most thermally
conductive material, intermediate for c-Si and GaAs that have intermediate
LTCs, and small for PbTe that has low LTC. Furthermore, the spectral LTC
results in Fig. 6 show that the force errors mainly reduce $\kappa(\omega)$ in
the low-frequency region. With increasing force errors, $\kappa(\omega)$ in
the low-frequency region is more and more reduced. This further supports the
large effect of the force errors in high-LTC materials, which usually have
large $\kappa(\omega)$ in the low-frequency region. Therefore, MLMD
simulations remain largely accurate for low-LTC materials, such as PbTe [20],
a-Si [23], amorphous SiO2 [24], and liquid water [25].

##  IV Conclusions

In summary, our systematic investigation revealed that the underestimation of
lattice thermal conductivity commonly observed in the literature are primarily
due to force fitting errors in machine learned potentials. Using empirical
potentials and Langevin thermostat we demonstrated that introducing random
forces on atoms can significantly reduce the lattice thermal conductivity,
supporting our hypothesis. These random forces act as an additional source of
phonon scattering, thereby reducing the lattice thermal conductivity.
Employing the kinetic theory of phonons and Matthiessen’s rule, we established
a linear extrapolation formula to estimate the thermal conductivity in the
absence of random forces. The validity of the extrapolation scheme was tested
using empirical potentials on various materials, including a-Si, Si-Ge alloys,
c-Si, c-Ge, graphene, and CNT.

We established that the force errors in machine-learned potentials follow a
Gaussian distribution, akin to the distribution of random forces in the
Langevin thermostat. This similarity inspired us to intentionally introduce
extra force noises via the Langevin thermostat and then extrapolate to the
limit of zero force error. The extrapolated results show excellent agreement
with experimental data over a broad temperature range for all the materials
studied. Spectral thermal conductivity analyses further indicate that the
underestimation of the lattice thermal conductivity is mainly due to increased
acoustic phonon scatterings caused by the force errors. Our findings provide a
clear explanation for the underestimated thermal conductivity often observed
in molecular dynamics simulations based on machine learned potentials. The
method of correcting this underestimation we developed will significantly
enhance the applicability of machine learned potentials in the prediction of
lattice thermal conductivity.

Data availability: All the training datasets and the trained NEPs models are
freely available at https://gitlab.com/brucefan1983/nep-data.

###### Acknowledgements.

X.W. and W.Z. contributed equally. This work is supported by the National
Natural Science Foundation of China (Grant No.12174276 and No. 52076002), the
New Cornerstone Science Foundation through the XPLORER PRIZE. The Center of
Campus Network and Modern Educational Technology of Guangdong University of
Technology, and the High-performance Computing Platform of Peking University
are acknowledged for providing computational resources and technical support
for this work. P.Y. is supported by the Israel Academy of Sciences and
Humanities & Council for Higher Education Excellence Fellowship Program for
International Postdoctoral Researchers.

## References

  * Moore and Shi [2014] A. L. Moore and L. Shi, Emerging challenges and materials for thermal management of electronics, Materials Today 17, 163 (2014).
  * Kang  _et al._ [2021] J. S. Kang, M. Li, H. Wu, H. Nguyen, T. Aoki, and Y. Hu, Integration of boron arsenide cooling substrates into gallium nitride devices, Nature Electronics 4, 416 (2021).
  * Snyder and Toberer [2008] G. J. Snyder and E. S. Toberer, Complex thermoelectric materials, Nature materials 7, 105 (2008).
  * Zhao  _et al._ [2014] L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid, and M. G. Kanatzidis, Ultralow thermal conductivity and high thermoelectric figure of merit in SnSe crystals, Nature 508, 373 (2014).
  * Zhou  _et al._ [2022] W. Zhou, Y. Dai, J. Zhang, B. Song, T.-H. Liu, and R. Yang, Effect of four-phonon interaction on phonon thermal conductivity and mean-free-path spectrum of high-temperature phase SnSe, Applied Physics Letters 121, 112202 (2022).
  * Perepezko [2009] J. H. Perepezko, The hotter the engine, the better, Science 326, 1068 (2009).
  * Vaßen  _et al._ [2010] R. Vaßen, M. O. Jarligo, T. Steinke, D. E. Mack, and D. Stöver, Overview on advanced thermal barrier coatings, Surface and Coatings Technology 205, 938 (2010).
  * Qian  _et al._ [2021] X. Qian, J. Zhou, and G. Chen, Phonon-engineered extreme thermal conductivity materials, Nature Materials 20, 1188 (2021).
  * Tadano and Tsuneyuki [2018] T. Tadano and S. Tsuneyuki, Quartic anharmonicity of rattlers and its effect on lattice thermal conductivity of clathrates from first principles, Physical Review Letters 120, 105901 (2018).
  * Hanus  _et al._ [2021] R. Hanus, R. Gurunathan, L. Lindsay, M. T. Agne, J. Shi, S. Graham, and G. Jeffrey Snyder, Thermal transport in defective and disordered materials, Applied Physics Reviews 8, 031311 (2021).
  * Gu  _et al._ [2021] X. Gu, Z. Fan, and H. Bao, Thermal conductivity prediction by atomistic simulation methods: Recent advances and detailed comparison, Journal of Applied Physics 130, 210902 (2021).
  * Lee  _et al._ [1991] Y. H. Lee, R. Biswas, C. M. Soukoulis, C. Z. Wang, C. T. Chan, and K. M. Ho, Molecular-dynamics simulation of thermal conductivity in amorphous silicon, Physical Review B 43, 6573 (1991).
  * Vogelsang  _et al._ [1987] R. Vogelsang, C. Hoheisel, and G. Ciccotti, Thermal conductivity of the lennard-jones liquid by molecular dynamics calculations, The Journal of chemical physics 86, 6371 (1987).
  * Green [1954] M. S. Green, Markoff random processes and the statistical mechanics of time‐dependent phenomena. II. Irreversible processes in fluids, The Journal of Chemical Physics 22, 398 (1954).
  * Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems, Journal of the physical society of Japan 12, 570 (1957).
  * Evans [1982] D. J. Evans, Homogeneous NEMD algorithm for thermal conductivity—Application of non-canonical linear response theory, Physics Letters A 91, 457 (1982).
  * Fan  _et al._ [2019] Z. Fan, H. Dong, A. Harju, and T. Ala-Nissila, Homogeneous nonequilibrium molecular dynamics method for heat transport and spectral decomposition with many-body potentials, Physical Review B 99, 064308 (2019).
  * Sosso  _et al._ [2012] G. C. Sosso, D. Donadio, S. Caravati, J. Behler, and M. Bernasconi, Thermal transport in phase-change materials from atomistic simulations, Phys. Rev. B 86, 104301 (2012).
  * Liu  _et al._ [2021] H. Liu, X. Qian, H. Bao, C. Y. Zhao, and X. Gu, High-temperature phonon transport properties of SnSe from machine-learning interatomic potential, Journal of Physics: Condensed Matter 33, 405401 (2021).
  * Cheng  _et al._ [2023] R. Cheng, X. Shen, S. Klotz, Z. Zeng, Z. Li, A. Ivanov, Y. Xiao, L.-D. Zhao, F. Weber, and Y. Chen, Lattice dynamics and thermal transport of pbte under high pressure, Phys. Rev. B 108, 104306 (2023).
  * Ying  _et al._ [2023] P. Ying, T. Liang, K. Xu, J. Zhang, J. Xu, Z. Zhong, and Z. Fan, Sub-Micrometer Phonon Mean Free Paths in Metal–Organic Frameworks Revealed by Machine Learning Molecular Dynamics Simulations, ACS Applied Materials & Interfaces 15, 36412 (2023).
  * Du  _et al._ [2023] P.-H. Du, C. Zhang, T. Li, and Q. Sun, Low lattice thermal conductivity with two-channel thermal transport in the superatomic crystal $\mathrm{P}{\mathrm{h}}_{4}\mathrm{Al}{\mathrm{br}}_{4}$, Physical Review B 107, 155204 (2023).
  * Wang  _et al._ [2023] Y. Wang, Z. Fan, P. Qian, M. A. Caro, and T. Ala-Nissila, Quantum-corrected thickness-dependent thermal conductivity in amorphous silicon predicted by machine learning molecular dynamics simulations, Physical Review B 107, 054303 (2023).
  * Liang  _et al._ [2023] T. Liang, P. Ying, K. Xu, Z. Ye, C. Ling, Z. Fan, and J. Xu, Mechanisms of temperature-dependent thermal transport in amorphous silica from machine-learning molecular dynamics, Phys. Rev. B 108, 184203 (2023).
  * Xu  _et al._ [2023] K. Xu, Y. Hao, T. Liang, P. Ying, J. Xu, J. Wu, and Z. Fan, Accurate prediction of heat conductivity of water by a neuroevolution potential, The Journal of Chemical Physics 158, 204114 (2023).
  * Korotaev  _et al._ [2019] P. Korotaev, I. Novoselov, A. Yanilkin, and A. Shapeev, Accessing thermal conductivity of complex compounds by machine learning interatomic potentials, Physical Review B 100, 144308 (2019).
  * Qian  _et al._ [2019] X. Qian, S. Peng, X. Li, Y. Wei, and R. Yang, Thermal conductivity modeling using machine learning potentials: application to crystalline and amorphous silicon, Materials Today Physics 10, 100140 (2019).
  * Bussi and Parrinello [2007] G. Bussi and M. Parrinello, Accurate sampling using langevin dynamics, Physical Review E 75, 056707 (2007).
  * Fan  _et al._ [2021] Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen, and T. Ala-Nissila, Neuroevolution machine learning potentials: Combining high accuracy and low cost in atomistic simulations and application to heat transport, Physical Review B 104, 104309 (2021).
  * Fan [2022] Z. Fan, Improving the accuracy of the neuroevolution machine learning potential for multi-component systems, Journal of Physics: Condensed Matter 34, 125902 (2022).
  * Fan  _et al._ [2022] Z. Fan, Y. Wang, P. Ying, K. Song, J. Wang, Y. Wang, Z. Zeng, K. Xu, E. Lindgren, J. M. Rahm, A. J. Gabourie, J. Liu, H. Dong, J. Wu, Y. Chen, Z. Zhong, J. Sun, P. Erhart, Y. Su, and T. Ala-Nissila, GPUMD: A package for constructing accurate machine-learned potentials and performing highly efficient atomistic simulations, The Journal of Chemical Physics 157, 114801 (2022).
  * Schaul  _et al._ [2011] T. Schaul, T. Glasmachers, and J. Schmidhuber, High dimensions and heavy tails for natural evolution strategies, in  _Proceedings of the 13th Annual Conference on Genetic and Evolutionary Computation_, GECCO ’11 (Association for Computing Machinery, New York, NY, USA, 2011) p. 845–852.
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical Review B 59, 1758 (1999).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Physical Review B 23, 5048 (1981).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Fan  _et al._ [2017] Z. Fan, W. Chen, V. Vierimaa, and A. Harju, Efficient molecular dynamics simulations with many-body potentials on graphics processing units, Computer Physics Communications 218, 10 (2017).
  * Martyna  _et al._ [1992] G. J. Martyna, M. L. Klein, and M. Tuckerman, Nosé–Hoover chains: The canonical ensemble via continuous dynamics, The Journal of Chemical Physics 97, 2635 (1992).
  * Bussi  _et al._ [2007] G. Bussi, D. Donadio, and M. Parrinello, Canonical sampling through velocity rescaling, The Journal of Chemical Physics 126, 014101 (2007).
  * Glassbrenner and Slack [1964] C. J. Glassbrenner and G. A. Slack, Thermal conductivity of silicon and germanium from 3 K to the melting point, Physical Review 134, A1058 (1964).
  * Shanks  _et al._ [1963] H. R. Shanks, P. D. Maycock, P. H. Sidles, and G. C. Danielson, Thermal conductivity of silicon from 300 to 1400 K, Physical Review 130, 1743 (1963).
  * Maycock [1967] P. Maycock, Thermal conductivity of silicon, germanium, III–V compounds and III–V alloys, Solid-State Electronics 10, 161 (1967).
  * Volz and Chen [2000] S. G. Volz and G. Chen, Molecular-dynamics simulation of thermal conductivity of silicon crystals, Physical Review B 61, 2651 (2000).
  * Tersoff [1988] J. Tersoff, Empirical interatomic potential for silicon with improved elastic properties, Physical Review B 38, 9902 (1988).
  * Amith  _et al._ [1965] A. Amith, I. Kudman, and E. F. Steigmeier, Electron and phonon scattering in GaAs at high temperatures, Physical Review 138, A1270 (1965).
  * Blakemore [1982] J. S. Blakemore, Semiconducting and other major properties of gallium arsenide, Journal of Applied Physics 53, R123 (1982).
  * Vega-Flick  _et al._ [2019] A. Vega-Flick, D. Jung, S. Yue, J. E. Bowers, and B. Liao, Reduced thermal conductivity of epitaxial GaAs on Si due to symmetry-breaking biaxial strain, Physical Review Materials 3, 034603 (2019).
  * Inyushkin  _et al._ [2003] A. V. Inyushkin, A. N. Taldenkov, A. Y. Yakubovsky, A. V. Markov, L. Moreno-Garsia, and B. N. Sharonov, Thermal conductivity of isotopically enriched 71GaAs crystal, Semiconductor Science and Technology 18, 685 (2003).
  * Holland [1964] M. G. Holland, Phonon scattering in semiconductors from thermal conductivity studies, Physical Review 134, A471 (1964).
  * Chen  _et al._ [2011] S. Chen, A. L. Moore, W. Cai, J. W. Suk, J. An, C. Mishra, C. Amos, C. W. Magnuson, J. Kang, L. Shi, _et al._ , Raman measurements of thermal transport in suspended monolayer graphene of variable sizes in vacuum and gaseous environments, ACS nano 5, 321 (2011).
  * Chen  _et al._ [2012] S. Chen, Q. Wu, C. Mishra, J. Kang, H. Zhang, K. Cho, W. Cai, A. A. Balandin, and R. S. Ruoff, Thermal conductivity of isotopically modified graphene, Nature Materials 11, 203 (2012).
  * Lee  _et al._ [2011] J.-U. Lee, D. Yoon, H. Kim, S. W. Lee, and H. Cheong, Thermal conductivity of suspended pristine graphene measured by raman spectroscopy, Physical Review B 83, 081419 (2011).
  * Cai  _et al._ [2010] W. Cai, A. L. Moore, Y. Zhu, X. Li, S. Chen, L. Shi, and R. S. Ruoff, Thermal transport in suspended and supported monolayer graphene grown by chemical vapor deposition, Nano letters 10, 1645 (2010).
  * Androulakis  _et al._ [2010] J. Androulakis, I. Todorov, D.-Y. Chung, S. Ballikaya, G. Wang, C. Uher, and M. Kanatzidis, Thermoelectric enhancement in PbTe with K or Na codoping from tuning the interaction of the light- and heavy-hole valence bands, Physical Review B 82, 115209 (2010).
  * Pei  _et al._ [2011] Y. Pei, X. Shi, A. LaLonde, H. Wang, L. Chen, and G. J. Snyder, Convergence of electronic bands for high performance bulk thermoelectrics, Nature 473, 66 (2011).


1Center for Computational Health, IBM Research, NY, USA;
2IBM Research, Nairobi, Kenya

# Disease Progression Modeling Workbench 360

Parthasarathy Suryanarayanan1 Prithwish Chakraborty1 Piyush Madan1 Kibichii
Bore2 William Ogallo2 Rachita Chandra1 Mohamed Ghalwash1 Italo Buleje1 Sekou
Remy2 Shilpa Mahatma1 Pablo Meyer1 Jianying Hu1

##  1 Background

Disease Progression Modeling (DPM) [1] aims to characterize the progression of
a disease and its comorbidities over time using a wide range of analytics
models including disease staging [2], patient trajectory analytics [3],
prediction [4], and time-to-event estimations [5] for key disease-related
events. DPM has applications throughout the healthcare ecosystem, from
providers (e.g., decision support for patient staging), to payers (e.g., care
management), and pharmaceutical companies (e.g., clinical trial enrichment).
But the complexity of building effective DPM models can be a road-block for
their rapid experimentation and adoption. Some of this is addressed by
standardization of data model and tooling for data analysis and cohort
selection [6]. However, there are still unmet needs to facilitate the
development of advanced machine learning techniques such as deep learning with
additional requirements such as experiment tracking and reproducibility [7].
Furthermore, to accelerate DPM research, a data scientist’s available tools
should include a framework for deploying models as cloud-ready microservices
for rapid prototyping and dissemination [8].

Figure 1: DPM360 Component View

In this demonstration, we introduce Disease Progression Modeling Workbench 360
(DPM360) open source project (https://ibm.github.io/DPM360/). DPM360 is an
easy-to-install system to help research and development of DPM models (Figure
1). It manages the entire modeling life cycle, from data analysis (e.g, cohort
identification) to machine learning algorithm development and prototyping.
DPM360 augments the advantages of data model standardization and tooling
(OMOP-CDM, Athena, ATLAS) provided by the widely-adopted OHDSI [9] initiative,
with a powerful machine learning training framework, and a mechanism for rapid
prototyping through automatic deployment of models as containerized services
into a cloud environment. This enables a quicker and flexible implementation
and validation of the models.

##  2 Methods

The architecture, shown in Figure 2 has four main components.

Figure 2: DPM360 Architecture

(1) Lightsaber: an extensible training framework which provides blueprints for
the development of disease progression models (DPM). It is designed ground up
using state-of-the art open source tools [10, 11] to provide a simple modular
and unified model training framework to support some of the common use cases
for DPM. Lightsaber contains four key modules:

  * •

data ingestion modules to support standardized methods of ingesting data

  * •

model trainers to support standardized model training incorporating best
practices

  * •

metrics to support pre-built DPM problem specific model evaluation

  * •

in-built model tracking and support for post-hoc model evaluation by
integrating with a Model Registry.

Users can select specific modules and integrate them into their modeling
workflow. Lightsaber also comes with a reusable library of state-of-the-art
machine and deep learning algorithms for DPM (e.g. LSTM [12] for in-hospital
mortality predictions).

Lightsaber integrates naturally with ATLAS using a client called Lightsaber
Client for ATLAS (LCA), enabling automated extraction of features from the
OMOP CDM model, thus complementing the ease and flexibility of defining
standardized cohorts using ATLAS graphical user interface with the ability to
quickly develop deep learning algorithms for DPM in Lightsaber using Python.
LCA can be configured with the cohort details, covariate settings, model
training settings for Lightsaber to extract the right set of features in
formats currently supported in the OHDSI stack (see FeatureExtraction and
PatientLevelPrediction R packages via the Rpy2 interface). Additionally, the
LCA uses custom queries and algorithms to extract and transform complex time
series features into formats required for DPM in Lightsaber. For each feature
extraction process, a YAML configuration file is automatically generated. This
file specifies outcomes, covariate types, and file locations of the extracted
feature files. Thus, Lightsaber allows a user to concentrate just on the logic
of their model as it takes care of the rest.

(2) Tracking provenance of all aspects of model building is essential for
trust and reproducibility - thus experiments ran using Lightsaber are
automatically tracked in a Model Registry including model parameters, problem
specific metrics, and model binaries allowing the identification of algorithms
and parameters that result in the best model performance.

(3) The Service Builder component automatically converts registered models in
Model Registry into microservices [13], through the usage of hooks that listen
for production ready models in the registry and thereafter start the model
packaging execution pipeline. The pipeline includes extraction of model and
its dependencies from the registry, containerization, and deployment in the
target cluster (Kubernates or OpenShift). Upon successful model deployment, a
callback function updates model metadata in the registry with deployment
status and model access endpoint. Using this endpoint, potential users (data
scientist or product manager) can interact with the model, now deployed as a
microservice, though a Swagger based interface.

(4) The Installer component installs the fully functional DPM360, including
OHDSI tools, Model Registry, and Service Builder into a Kubernetes or
OpenShift cluster using Helm charts. Each of these components are run as
services within the cluster. The implementation also uses Persistent Volume
Claims to persist the data (e.g: model artifacts, ATLAS database files etc.).

##  3 Results

A detailed description of the target use cases, architecture and road map is
available in our public GitHub repository at https://ibm.github.io/DPM360/
which contains our initial implementation. To verify the integration of the
overall system, we reproduced the results from clinical prediction benchmarks
from MIMIC III [14]. Specifically, we defined and extracted cohorts related to
in-hospital mortality in critical care settings from an OMOP CDM version of
the MIMIC III dataset [15]. Using ATLAS, we defined our target cohort
consisting of adult patients who have been hospitalized for the first time for
at least two days, and who have at least one measurement recorded within the
first 48 hours, and our outcome cohort of adult patients who died in hospital
within 30 days of the first admission. We subsequently used the LCA to extract
the features for the prediction task, and then trained a model using the
Lightsaber framework. We obtained comparable metrics to the benchmark, thus
proving the soundness of the overall system. Also, independently, the
Lightsaber component has been successfully used for rapid experimentation and
analysis of real-world data [4, 3].

##  4 Conclusion

In the coming months, we plan to execute our published road-map for DPM360 and
report on its effectiveness at improving collaborative research and rapid
commercialization from several ongoing DPM collaborations.

## Acknowledgements

We would like to thank Divya Pathak and Daby Sow for their guidance, and
Sundar Saranthan for his contributions.

## References

  * 1. Wang X, Qi J, Yang Y, Yang P.  A Survey of Disease Progression Modeling Techniques for Alzheimer’s Diseases.  In: 2019 IEEE 17th International Conference on Industrial Informatics (INDIN). vol. 1. IEEE; 2019. p. 1237–1242.
  * 2. Sun Z, Ghosh S, Li Y, Cheng Y, Mohan A, Sampaio C, et al.  A probabilistic disease progression modeling approach and its application to integrated Huntington’s disease observational data.  JAMIA open. 2019;2(1):123–130.
  * 3. Dey S, Bose A, Chakraborty P, Ghalwash M, Saenz AG, Ultro F, et al.. Impact of Clinical and Genomic Factors on SARS-CoV2 Disease Severity. Cold Spring Harbor Laboratory Press; 2021.  Accepted at AMIA 2021 Annual Symposium. To appear, early access at https://www.medrxiv.org/content/10.1101/2021.03.15.21253549v1.
  * 4. Chakraborty P, Codella J, Madan P, Li Y, Huang H, Park Y, et al.. Blending Knowledge in Deep Recurrent Networks for Adverse Event Prediction at Hospital Discharge ; 2021.  Presented at AMIA 2021 Virtual Informatics Summit. To appear, early access at https://arxiv.org/abs/2104.04377.
  * 5. Liu B, Li Y, Sun Z, Ghosh S, Ng K.  Early prediction of diabetes complications from electronic health records: A multi-task survival analysis approach.  In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 32; 2018. .
  * 6. Hripcsak G, Ryan PB, Duke JD, Shah NH, Park RW, Huser V, et al.  Characterizing treatment pathways at scale using the OHDSI network.  Proceedings of the National Academy of Sciences. 2016;113(27):7329–7336.
  * 7. McDermott M, Wang S, Marinsek N, Ranganath R, Ghassemi M, Foschini L.  Reproducibility in machine learning for health.  arXiv preprint arXiv:190701463. 2019.
  * 8. Mäkinen S, Skogström H, Laaksonen E, Mikkonen T.  Who Needs MLOps: What Data Scientists Seek to Accomplish and How Can MLOps Help?  arXiv preprint arXiv:210308942. 2021.
  * 9. Hripcsak G, Duke JD, Shah NH, Reich CG, Huser V, Schuemie MJ, et al.  Observational Health Data Sciences and Informatics (OHDSI): opportunities for observational researchers.  Studies in health technology and informatics. 2015;216:574.
  * 10. Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al.  Scikit-learn: Machine Learning in Python.  Journal of Machine Learning Research. 2011;12:2825–2830.
  * 11. Falcon W, et al.  PyTorch Lightning.  GitHub Note: https://githubcom/PyTorchLightning/pytorch-lightning. 2019;3.
  * 12. Gers FA, Schmidhuber J, Cummins F.  Learning to forget: Continual prediction with LSTM.  Neural computation. 2000;12(10):2451–2471.
  * 13. Dragoni N, Giallorenzo S, Lafuente AL, Mazzara M, Montesi F, Mustafin R, et al.  Microservices: yesterday, today, and tomorrow.  Present and ulterior software engineering. 2017:195–216.
  * 14. Harutyunyan H, Khachatrian H, Kale DC, Ver Steeg G, Galstyan A.  Multitask learning and benchmarking with clinical time series data.  Scientific Data. 2019 Jun;6(1).  Available from: http://dx.doi.org/10.1038/s41597-019-0103-9.
  * 15. Johnson AE, Pollard TJ, Shen L, Li-Wei HL, Feng M, Ghassemi M, et al.  MIMIC-III, a freely accessible critical care database.  Scientific data. 2016;3(1):1–9.