# Predicting Properties of Periodic Systems from Cluster Data: A Case Study of
Liquid Water

Viktor Zaverkin  Institute for Theoretical Chemistry, University of Stuttgart,
Pfaffenwaldring 55, 70569 Stuttgart, Germany  David Holzmüller  Institute for
Stochastics and Applications, University of Stuttgart, Pfaffenwaldring 57,
70569 Stuttgart, Germany  Robin Schuldt  Institute for Theoretical Chemistry,
University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany
Johannes Kästner<EMAIL_ADDRESS>Institute for Theoretical
Chemistry, University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart,
Germany

###### Abstract

The accuracy of the training data limits the accuracy of bulk properties from
machine-learned potentials. For example, hybrid functionals or wave-function-
based quantum chemical methods are readily available for cluster data but
effectively out-of-scope for periodic structures. We show that local, atom-
centred descriptors for machine-learned potentials enable the prediction of
bulk properties from cluster model training data, agreeing reasonably well
with predictions from bulk training data. We demonstrate such transferability
by studying structural and dynamical properties of bulk liquid water with
density functional theory and have found an excellent agreement with
experimental as well as theoretical counterparts.

##  I  Introduction

The ability to perform simulations that allow for bond formation and bond
breaking in aqueous environments is necessary to understand several chemical
processes.[1] While first-principles sampling can provide some insight, it is
usually limited by its high computational cost. Approximate methods allow for
simulations taking into account the large size of conformational and chemical
space of interest. Combining molecular dynamics simulations with machine-
learned potentials (MLPs), which have been proven to learn any complex non-
linear relationship, enables us to extend the simulation time scales and
system sizes. These factors are pivotal for obtaining the relevant time
correlation functions to converge various dynamical properties in solution,
ranging from diffusion coefficients and reorientation times to reaction rates
and spectroscopy.

The advent of machine learning methods in computational chemistry and physics
made it feasible to construct efficient interatomic potentials with accuracy
on par with reference ab initio methods.[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24] Neural networks (NNs),
especially, have been proven to approximate any non-linear functional
relationship.[25] This fact has promoted their popularity in computational
chemistry and materials science.[26, 27, 28, 29, 30, 31, 32, 33, 34] NNs were
initially applied to represent potential energy surfaces (PESs) of small
atomistic systems[35, 36] and were later extended to high-dimensional
systems.[2] Once trained, the computational cost of machine-learned potentials
(MLPs) based on NNs is independent of the number of data points used for
training as opposed to kernel-based models.[4, 5, 6, 7]

Figure 1: Schematic illustration of conventional and investigated
physicochemical property prediction path for a periodic system. Path A: An
interatomic potential for a periodic system is created using periodic
reference data. Path B: An interatomic potential for a periodic system is
built using cluster reference data.

Typically, MLPs are trained following path A in Fig. 1, i.e. models supposed
to predict periodic systems use periodic reference data for training and vice
versa. Special reference ab initio methods might be required for some
applications, e.g., hybrid functionals such as B3LYP[37, 38] or wave-function-
based approaches, like coupled-cluster methods. For these, the generation of
training data sets for periodic systems can be computationally prohibitive.
Taking into account the local, atom-centred architecture of most machine
learning approaches employed in computational chemistry and physics, it is, in
principle, possible to train the models on cluster data and use them at the
inference step for the prediction of periodic systems, see path B in Fig. 1.
Such an approach aims to reduce the computational cost for the construction of
highly accurate machine-learned interatomic potentials for periodic systems.

Previous work[39, 40] employed an architecture that computes total energies
and, thus, atomic forces from a charge equilibration process.[41, 42, 43]
However, here we envisage studying the direct transferability of properties
obtained employing the locality assumption from small clusters to bulk
materials, rather than using them for an intermediate step to get the global
charge density determined by atomic environment-dependent
electronegativities.[39, 40] Recently, the transferability of interatomic MLPs
has attracted much interest by the community,[44, 45, 46, 47, 48, 49]
resulting in an excellent generalisation ability of some MLPs on the out-of-
sample configurations. A few cases[44, 48, 49] demonstrate the transferability
of models trained on smaller molecular fragments to larger atomistic systems.
In Ref. 44, in particular, NN-based MLPs, for which the total energy of a
system is built from the environment-dependent atomic energies, have been
compared to fragmentation-based approaches.[50, 51, 52, 53, 54]

In this work, we investigate the transferability of MLPs trained on local
reference data (water clusters) to a periodic liquid water system.
Specifically, we use the open-source package Gaussian moment neural network
(GM-NN) for constructing MLPs and applying them in atomistic simulations. The
GM-NN source code is available free-of-charge from gitlab.com/zaverkin_v/gmnn.
Note that no difference between MLPs for water with and without explicit long-
range electrostatics has been found for a cutoff radius of $6.35$ Å. [26]
Therefore, all models used in this work are built employing a cutoff of $6.5$
Å but do not contain the explicit long-range electrostatics.

To assess the transferability of cluster models, we thoroughly benchmark them
comparing the mean absolute errors and respective error distributions to the
respective results obtained by employing models trained on periodic data. A
more rigorous evaluation of the transferability of local models is achieved by
running long molecular dynamics (MD) simulations in canonical (NVT)
statistical ensemble at 300 K. As a baseline for the assessment, we consider
structural properties such as oxygen–oxygen, oxygen–hydrogen, and
hydrogen–hydrogen radial distribution functions. Additionally, we compute the
water self-diffusion coefficient employing local and periodic MLPs. Besides
the NVT simulations, the equilibrium density of liquid water at ambient
conditions (at a temperature of 300 K and a pressure of 1 bar) has been
estimated by running NPT (constant pressure) simulations. All results are
carefully compared to those obtained from training on bulk data and to
experiment.

While we apply NN-based interatomic potentials to the example of liquid water,
many conventional force fields exist to perform accurate simulations of the
latter. One of the recent examples is the TL4P water model [55] which has been
designed to describe the condensed-phase and spectroscopic properties of
liquid water correctly. [56] However, there are other non-conventional force
fields, [56] which allow building bulk water models from small clusters. [57,
58, 59, 60, 61, 62, 63] For a detailed overview of non-conventional force
fields, see elsewhere. [56] Here, we compare the estimated properties of bulk
liquid water by employing cluster models to the fragment-based coupled-cluster
approach [62] and the well-established TIP4P water model.[64] However, the
main goal of this work is to investigate the transferability of MLPs trained
on cluster data to periodic systems. Therefore, a detailed comparison of their
accuracy in terms of predicted energy and forces is out of the scope of this
work and the reader is referred to, e.g., Ref. 44.

##  II  Computational Details

###  II.1  Data Sets

This work demonstrates that atom-centred machine-learned potentials trained on
cluster data can be transferred to periodic systems without loss of
generalisation capability. For this purpose, we employ two kinds of data sets.
The first data set contains structures, energies, and atomic forces of a
periodic water box with 64 water molecules (192 atoms) at the revPBE-D3/PAW
level of theory.[65, 66, 67] This data set was generated by running an ab
initio molecular dynamics (AIMD) simulation at a temperature of 400 K.[23] We
refer to this data set and the respective interatomic potentials as
revPBE-D3-bulk. The respective data set is available free-of-charge from Ref.
68.

The second kind of data contains larger water clusters extracted from the
N-ASW data set, created to study nitrogen atoms’ adsorption and desorption
dynamics on top of amorphous solid water (ASW).[29] The N-ASW data set is
available free-of-charge from Ref. 69. In total, we extracted 1520 structures
from the original N-ASW data set: 800 structures with 30 H2O sampled at 150 K,
400 structures with 30 H2O sampled at 800 K, 200 structures with 60 H2O
sampled at 300 K, 100 structures with 93 H2O sampled at 300 K, 10 structures
with 109 H2O sampled at 300 K, and 10 structures with 126 H2O sampled at 300
K.

To allow for a comparison with the models trained on the revPBE-D3-bulk data
set, energies and atomic forces of the respective 1520 configurations were
recalculated at the revPBE-D3/def2-TZVP level of theory[65, 66, 67, 70, 71]
using Turbomole 7.4[72] within ChemShell.[73, 74] We refer to the respective
interatomic potentials and the data set as revPBE-D3-cluster. In this work, we
are especially interested in using hybrid DFT functionals for the construction
of training data which is typically infeasible for periodic structures with a
reasonable effort, even for liquid water. Therefore, we compute energies and
atomic forces at the B3LYP-D3/def2-TZVP level of theory,[37, 38, 67, 70, 71]
popular in the computational chemistry community and refer to the respective
data set as B3LYP-D3-cluster. In addition to the mentioned functionals, we
construct a data set employing the BLYP-D3/def2-TZVP level of theory[37, 75,
67, 70, 71] and refer to it as BLYP-D3-cluster. The respective data sets are
available free-of-charge from Ref. 76.

###  II.2  Construction of MLPs

To study the transferability of atom-centred machine-learned potentials (MLP)
trained on cluster data to periodic structures, we employ molecular dynamics
(MD) simulations in combination with the Gaussian moment neural network (GM-
NN) proposed in Refs. 15, 24. GM-NN potentials are atom-centered interatomic
potentials for which the total energy $E$ of a system
$S=\\{\mathbf{r}_{i},Z_{i}\\}_{i=1}^{N_{\mathrm{at}}}$, with
$\mathbf{r}_{i}\in\mathbb{R}^{3}$ being the Cartesian coordinates and
$Z_{i}\in\mathbb{Z}$ being the atomic numbers, is decomposed into its atomic
contributions $E_{i}$

$E\left(S,\bm{\theta}\right)=\sum_{i=1}^{N_{\mathrm{at}}}E_{i}\left(\mathbf{G}_{i},\bm{\theta}\right).$
(1)

GM-NN models use neural networks (NNs) to map symmetry-preserving local atomic
descriptors $\mathbf{G}_{i}$, Gaussian moments (GMs), to auxiliary atomic
energy and include both the geometric and alchemical information about the
atomic species of both the central and neighbour atoms. Only a single NN has
to be trained for all atomic energy contributions, in contrast to using an
individual NN for each species as frequently done in the literature.[2]

In this work, we use a neural network with two hidden layers consisting of
$512$ nodes each and an input dimension of $360$ (number of invariant
features). We minimise the combined loss function

$\begin{split}\mathcal{L}_{E,\mathbf{F}}\left(\bm{\theta}\right)=\sum_{k=1}^{N_{\mathrm{Train}}}&\Bigg{[}\lambda_{E}\lVert
E_{k}^{\mathrm{ref}}-E(S_{k},\bm{\theta})\rVert_{2}^{2}+\\\
&\quad\frac{\lambda_{F}}{3N_{\mathrm{at}}^{(k)}}\sum_{i=1}^{N_{\mathrm{at}}^{(k)}}\lVert\mathbf{F}_{i,k}^{\mathrm{ref}}-\mathbf{F}_{i}\left(S_{k},\bm{\theta}\right)\rVert_{2}^{2}\Bigg{]},\end{split}$
(2)

to optimize the respective parameters of trainable representation, fully-
connected neural network part, and the parameters which scale and shift the
output of the neural network. For more details on the exact architecture of
the network, on the atom-centered representation, and on the definition of
energy scale and shift parameters see Ref. 24.

In Eq. (2), $N_{\mathrm{at}}^{(k)}$ is the number of atoms in the respective
structure. The reference values for the energy and atomic force are denoted by
$E_{k}^{\mathrm{ref}}$ and $\mathbf{F}_{i,k}^{\mathrm{ref}}$, respectively.
Since atomic forces containing $3N_{\mathrm{at}}$ scalars provide much more
information about a structure and determine the dynamics of a system, the
parameters $\lambda_{E}$ and $\lambda_{F}$ were set to $1~{}\text{au}$ and
$12N_{\mathrm{at}}^{(k)}$ au Å2, respectively. The loss function in Eq. (2)
was minimised employing the Adam optimiser[77] with $32$ molecules per mini-
batch. The layer-wise learning rate was set to $0.03$ for the parameters of
the fully connected layers, $0.02$ for the trainable representation, 0.05 and
0.001 for the shift and scale parameters of atomic energies, respectively. All
learning rates were decayed linearly to zero. For more information, see Ref.
24. Besides the hyper-parameters defined above, we have used a cutoff radius
of 6.5 Å employed in the definition of the GM representation, see Ref. 26.

To run MD simulations with the GM-NN interatomic potentials, we interfaced
them with the ASE package (v. 3.21.0).[78] To assess the accuracy of GM-NN
during an MD simulation, its uncertainty has been computed by employing the
ensembling technique.[79, 80, 81, 82, 83, 84, 85] Note that this method is
often referred to as the query-by-committee (QbC) approach[86] and is used to
run active learning, different to the current work. Here, we trained a
committee of 3 models on the same split of the data, i.e. on the same training
data set, but using randomly initialised parameters and reported on the
respective uncertainty

$\sigma_{\mathrm{ens}}\left(S\right)=\sqrt{\frac{1}{N_{\mathrm{ens}}}\sum_{i=1}^{N_{\mathrm{ens}}}\left(y_{i}\left(S\right)-\bar{y}\left(S\right)\right)^{2}},$
(3)

where $N_{\mathrm{ens}}$ is the number of models in the committee, generally
we use $N_{\mathrm{ens}}=3$.
$\bar{y}\left(S\right)=1/N_{\mathrm{ens}}\sum_{i=1}^{N_{\mathrm{ens}}}y_{i}\left(S\right)$
is the mean of the property prediction, the atomic force component in this
specific case, over the committee. All models were trained within the
Tensorflow framework[87] on an NVIDIA Tesla V100-SXM-32GB GPU.

###  II.3 Stress Tensor

The molecular dynamics (MD) simulations at constant pressure (NPT) in Section
III.3 require virial stress tensors to compute the instantaneous pressure.[88]
Thus, we have to derive an expression for stress tensors specific to the atom-
centred neural networks (NNs) employed in this work. Recall the original
expression derived earlier for arbitrary pair potentials[88] and used more
recently for MLPs[89, 90]

$\mathbf{W}=-\frac{1}{2}\sum_{i=1}^{N_{\mathrm{at}}}\sum_{j=1}^{N_{\mathrm{at}}}\mathbf{r}_{ij}\otimes\mathbf{F}_{ij},$
(4)

where $\otimes$ denotes the tensor product, $\mathbf{F}_{ij}$ is the force on
atom $i$ caused by the interaction with atom $j$, and
$\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ is the atomic distance vector.
Here, $\mathbf{F}_{ij}$ is defined as

$\mathbf{F}_{ij}=-\bm{\nabla}_{\mathbf{r}_{j}}E_{i},$ (5)

where $E_{i}$ is the atomic energy computed by an atom-centered MLP and
$\mathbf{r}_{j}$ are the coordinates of atom $j$.

The atomic energy in the Gaussian moment neural network (GM-NN) framework is
defined by the atomic distance vectors, i.e.,
$E_{i}=E_{i}\left(\\{\mathbf{r}_{ij}\\}_{j\in r_{\mathrm{max}}}\right)$. Here,
$r_{\mathrm{max}}$ denotes the cutoff radius. Thus, employing the chain rule,
the expression in Eq. (5) can be simplified to

$\mathbf{F}_{ij}=\bm{\nabla}_{\mathbf{r}_{ij}}E\left(S,\bm{\theta}\right),$
(6)

where $E\left(S,\bm{\theta}\right)=\sum_{i=1}^{N_{\mathrm{at}}}E_{i}$ is the
total energy of an atomic arrangement. The resulting expression for the stress
tensor reads

$\mathbf{W}(S,\bm{\theta})=-\sum_{i=1}^{N_{\mathrm{at}}}\sum_{j\in
r_{\mathrm{max}}}\mathbf{r}_{ij}\otimes\bm{\nabla}_{\mathbf{r}_{ij}}E(S,\bm{\theta}).$
(7)

Here, we skip the pre-factor $1/2$ since the second sum in Equation (7) runs
only over the local neighbourhood of atom $i$. Therefore, we do not have to
account for the double-counting of interactions. Periodic boundary conditions
are taken into account when computing $\mathbf{r}_{ij}$. The atoms within the
cutoff radius can originate from the image atoms and the atoms within the
cell.

##  III  Results and Discussion

###  III.1  Interatomic Neural Network Potentials

In this section, we discuss the possibility of transferring the interatomic
potentials trained on cluster data to periodic systems based on simple error
estimates such as absolute error distribution in atomic forces as well as the
respective directional error. However, prior to that, the performance of the
GM-NN interatomic potentials on each data set in terms of their mean absolute
error (MAE) is studied.

Figure 2: Learning curves for periodic and local water data sets. The mean
absolute errors (MAEs) of total energies and atomic forces are plotted against
the training set size. Linear fits are displayed for clarity, and shaded areas
denote the 95 % confidence intervals for linear regression. The dashed black
lines represent the desired accuracies of 1 kcal/mol and 1 kcal/mol/Å,
respectively.

Fig. 2 shows the learning curves for models trained on periodic revPBE-D3-bulk
and local revPBE-D3-cluster, BLYP-D3-cluster, B3LYP-D3-cluster data sets. From
the figure, we see that the GM-NN interatomic potentials trained on
revPBE-D3-bulk reach the desired accuracy of 1 kcal/mol for total energies
already for a training set size of about 1000 reference structures. The
desired accuracy of 1 kcal/mol/Å for atomic forces is reached already after
training on 100 reference structures. A similar trend can be observed for GM-
NN potentials trained on local data sets. However, an MAE of 1 kcal/mol in
total energies is reached slightly earlier when training on about 800
configurations.

To test the transferability of the GM-NN interatomic potentials trained on
local data to periodic systems, we predict total energies and atomic forces
using the local GM-NN models for the periodic revPBE-D3-bulk test data and
compare the acquired results with the predictions of periodic GM-NN models.
Note that we perform a direct comparison only for the local models trained on
the revPBE-D3-cluster data set since for the other local data sets, similar
results are expected.

Figure 3: (Left) Correlation of the predicted potential energies with the
corresponding reference values for all structures in the revPBE-D3-bulk test
data. An arbitrary energy shift is introduced to the energies predicted by the
GM-NN model trained on local revPBE-D3-cluster data to allow for a comparison
with the models trained on the periodic revPBE-D3-bulk data. (Right)
Distribution of the absolute force error for GM-NN potentials trained on local
and periodic data. Here, $F_{ij,k}$ is the force component $j\in\\{1,2,3\\}$
acting on the atom $i\in\\{1,\dots,N_{\mathrm{at}}^{(k)}\\}$ in the structure
$k\in\\{1,\dots,N_{\mathrm{test}}\\}$. Periodic models were trained using 1800
reference structures; local models were trained on 1400 reference structures.

Fig. 3 (left) shows the correlation of the GM-NN potential energies with the
corresponding ab initio reference energies obtained for the structures in the
revPBE-D3-bulk test data split. As seen in the figure, the GM-NN models
trained on data containing only water clusters can reproduce the total energy
of a periodic bulk water system with 64 water molecules (192 atoms) up to an
arbitrary shift. The latter is introduced mainly by the differences in ab
initio methods employed in Section II.1. Note that we shifted the
revPBE-D3-cluster-based GM-NN potential energies in Fig. 3 (left) to allow for
a better comparison with revPBE-D3-bulk-based GM-NN predictions.

Fig. 3 (right) shows the distribution of the absolute error in the components
of the predicted atomic forces for the structures within the revPBE-D3-bulk
test data set, employing MLPs trained on the local and periodic data. From the
figure, one can see that the distribution is broader for the MLPs trained on
local data compared to the models trained on periodic data. However, that may
originate from the differences in levels of theory and other settings used to
generate the training data. Additionally, MLPs trained on periodic data can
explicitly model the interactions with images from neighbouring cells, while
local models are not trained to do so. Therefore, we study this phenomenon
more rigorously, e.g., by evaluating the error in directions of atomic forces
and by evaluating structural and dynamical properties from a molecular
dynamics (MD) simulation.

Figure 4: Samples of the angular error of the atomic force in degrees versus
the absolute atomic force for GM-NN models trained on local revPBE-D3-cluster
and periodic revPBE-D3-bulk data. The colours indicate the relative density in
log-space as estimated by a kernel density estimator. The gray lines
correspond to $F\mapsto\theta_{\mathrm{max}}(E,F)$ (in degrees) for different
values of $E$, see Eq. (8).

Fig. 4 shows a corresponding analysis of the error in the direction of the
predicted atomic forces. The plots show, for different atoms in the
revPBE-D3-bulk test set, the absolute reference force and the angle difference
between reference and predicted forces. Brighter colours indicate a higher
density of points in the log-log plot. It can be seen that the maximal angular
error decreases with increasing absolute values of the force vectors, for both
MLPs trained on local and those trained on periodic data. This can be
explained by the loss function trying to minimize the absolute force error.
For a given reference absolute force
$F=\|\mathbf{F}^{\mathrm{ref}}_{i,k}\|_{2}$ and given absolute force error
$E=\|\mathbf{F}^{\mathrm{ref}}_{i,k}-\mathbf{F}_{i}(S_{k},\mathbf{\theta})\|_{2}$,
the maximum possible angle is

$\theta_{\mathrm{max}}(E,F)=\begin{cases}\arcsin(E/F)&,\text{ if $E\leq F$}\\\
\pi&,\text{ if $E>F$.}\end{cases}$ (8)

The grey lines in the figure correspond to
$F\mapsto\theta_{\mathrm{max}}(E,F)$ for different values of $E$. It can be
seen that the maximum angles are well-aligned with these curves.

As seen from the figure, the direction of the force vector scatters slightly
stronger for the MLP trained on local data compared to the one trained on
periodic data. By Eq. (8), this is in line with the larger force errors shown
in Fig. 3. Since the directional errors are typically still small for the MLP
trained on water clusters, it can be expected that such MLPs apply to a broad
range of atomistic simulations of periodic liquid water.

###  III.2  Water Properties

The MAE of total energies and atomic forces and the respective error
distributions discussed earlier are abstract quality measures of machine-
learned interatomic potentials. In practice, the robustness and reliability of
the potential in a real-time application are more important. Therefore, in
this section, we assess the robustness and smoothness of the GM-NN interatomic
potentials based on local data sets by applying them to a periodic liquid
water system containing 216 water molecules (648 atoms) to run a molecular
dynamics (MD) simulation. This is because for an MD simulation, a smooth
energy surface is required to obtain an adequate estimate for structural
properties such as radial distribution functions and dynamical properties like
diffusion coefficients.

In this section, in order to analyse properties of liquid water, we run
molecular dynamics (MD) simulations in the canonical (NVT) statistical
ensemble carried out within the ASE simulation package[78] using a Langevin
thermostat at the temperature of 300 K but different atomic densities
depending on the level of theory employed. For revPBE-D3-bulk and
revPBE-D3-cluster, we have used the atomic density of liquid water at ambient
conditions of $0.998$ g cm-3. For BLYP-D3-cluster and B3LYP-D3-cluster, we
employed slightly higher values of 1.04 and 1.06 g cm-3, obtained by
performing several experiments in which the atomic density was varied between
$0.998$ g cm-3 and 1.06 g cm-3. Note that it is known that BLYP-D3-based MD
simulations produce slightly denser liquid water,[91] and the present study
shows that the same can be expected for B3LYP-D3-based simulations. All MD
runs were performed over 2.0 ns using a time step of 0.5 fs. The atomic
velocities were initialised with a Maxwell–Boltzmann distribution for a
temperature of 300 K.

Figure 5: Uncertainty distribution for atomic force components obtained during
molecular dynamics simulation on bulk liquid water obtained via the QbC
approach[86] for models trained on periodic revPBE-D3-bulk and local
revPBE-D3-cluster, BLYP-D3-cluster, and B3LYP-D3-cluster reference data sets.

Forces for molecular dynamics are generated by an ensemble of three machine-
learned interatomic potentials, see Eq. (3). The ensembling technique[79, 80,
81, 82, 83, 84, 85] provides us with an error estimate of the potential during
simulation. We have found the machine-learned potential to be very accurate
with uncertainty between models with a mean ranging from $0.17$ kcal/mol/Å to
$0.28$ kcal/mol/Å, see Fig. 5. Due to the cutoff radius, a model cannot
directly notice whether periodic boundary conditions are used or not. Hence,
the larger uncertainty of local models indicates that the bulk data is locally
different from cluster data used during training, requiring the local models
to extrapolate. These differences might arise because cluster data contain
fewer bulk-like atoms compared to periodic data and also by other aspects of
training data generation. Interestingly, we have observed only a minor
increase in uncertainty of a factor of a maximum of 1.6, leading to a still
well-behaving model with a mean uncertainty of $0.28$ kcal/mol/Å. These
observations suggest that the models trained on local data do not reach a
strongly extrapolative region, making them suitable for the following study of
the liquid water system.

Note that an ensemble of three models can be deemed too small, and larger
uncertainties could be expected. The respective values can be re-scaled to get
more accurate estimates; see Ref. 85. However, one should consider that the
uncertainty generally does not have to correlate with the actual error, i.e.,
a large uncertainty does not necessarily imply a large error. For example, in
Ref. 92, an upper bound for the correlation of uncertainty of atomistic NNs
with the actual error has been computed. Here, the estimated uncertainties
should only shed light on the relative stability of MLPs during MD
simulations. A more rigorous assessment of the quality of trained MLPs is
achieved by analyzing the estimated properties, i.e., the radial distribution
functions and diffusion coefficients.

####  III.2.1  Radial Distribution Function

In the following, our discussion on liquid water will be concerned with the
structure of bulk liquid water and the dynamics of the water molecules. We
shall pay particular attention to the three radial distribution functions
(RDFs) $g_{\alpha\beta}\left(r\right)$ and the self-diffusion coefficient $D$,
all of which are available from experiments. All these properties are obtained
from the aforementioned 2.0 ns long MD simulations in the canonical (NVT)
ensemble after the 20.0 ps equilibration period. We split the whole trajectory
into ten equally-sized parts to get statistics on the estimated RDF curves.

Figure 6: Oxygen–oxygen ($g_{\ce{OO}}$), oxygen–hydrogen ($g_{\ce{OH}}$),
hydrogen–hydrogen ($g_{\ce{HH}}$) radial distribution functions (RDF) of
liquid water at ambient conditions from experiment and from MD simulations
based on GM-NN models trained on periodic revPBE-D3-bulk, local
revPBE-D3-cluster, BLYP-D3-cluster, and B3LYP-D3-cluster data. Experimental
results are taken from Refs. 93, 94 for $g_{\ce{OO}}$, $g_{\ce{OH}}$, and
$g_{\ce{HH}}$, respectively.

As an initial orientation, we refer to Fig. 6 showing the oxygen–oxygen
($g_{\ce{OO}}$), oxygen–hydrogen ($g_{\ce{OH}}$), and hydrogen–hydrogen
($g_{\ce{HH}}$) RDFs for liquid water at ambient conditions obtained from
experiment[93, 94] and from MD simulations employing GM-NN interatomic
potentials. The experimental data for the oxygen–oxygen ($g_{\ce{OO}}$) RDF
originates from high-energy X-ray diffraction measurements performed at
temperatures close to 297 K.[93]

As shown in Fig. 6, the $g_{\ce{OO}}$ curve simulated employing the
revPBE-D3-bulk-based potential is in good agreement with the experimental
results for both the positions and intensities of the first two peaks, see
Table 1. The respective values obtained from MD simulation are very close to
the experimental ones and deviate from them by at most 0.03 Å in positions and
0.20 in intensities. One should note that employing the local variant
revPBE-D3-cluster leads to an even better agreement of $g_{\ce{OO}}$ with its
experimental counterpart, especially considering the higher values of the
trough between the first two peaks. The results are also consistent with the
ab initio study in Ref. 1.

The BLYP-D3-cluster and B3LYP-D3-cluster based GM-NN potentials produce fairly
good $g_{\ce{OO}}$ curves with an excellent minimum between the first two
peaks but an overestimated intensity of the first peak. Positions of all peaks
are in good agreement with the experiment. Note that such behaviour is known
for the BLYP-D3 level of theory,[91] and interestingly is also observed for
the hybrid B3LYP-D3 functional.

Table 1: Positions and intensities of the 1st maximum $(r_{1},g_{1})$, 1st minimum $(r_{2},g_{2})$, and 2nd maximum $(r_{3},g_{3})$ in $g_{\ce{OO}}$ with the respective standard deviations given in parentheses. Experimental data is from high-energy x-ray diffraction measurements.[93] Positions are given in Å. | Exp. X-ray | revPBE-D3-bulk | revPBE-D3-cluster | BLYP-D3-cluster | B3LYP-D3-cluster
---|---|---|---|---|---
$r_{1}$ | 2.80(1) | 2.79(1) | 2.82(1) | 2.79(1) | 2.78(1)
$g_{1}$ | 2.57(5) | 2.77(2) | 2.55(2) | 2.83(2) | 2.87(2)
$r_{2}$ | 3.45(4) | 3.43(2) | 3.54(3) | 3.41(3) | 3.37(2)
$g_{2}$ | 0.84(2) | 0.80(2) | 0.88(1) | 0.85(2) | 0.85(1)
$r_{3}$ | 4.5(1) | 4.47(3) | 4.51(3) | 4.47(3) | 4.42(3)
$g_{3}$ | 1.12(2) | 1.18(1) | 1.10(1) | 1.13(1) | 1.12(1)

All oxygen–oxygen RDFs, experimental and the ones sampled from MD simulations,
show a clear considerable disorder, since the $g_{\ce{OO}}$ curve is quite
close to unity over the range $3.0<r_{\ce{OO}}<4.0$ Å. The diffraction
measurements agree that the intensity of $g_{\ce{OO}}$ at its first minimum is
$0.84\pm 0.02$. All the GM-NN-based simulations resulted in only slightly
over- or under-structured liquids with the height of the first peak and the
first minimum in $g_{\ce{OO}}$ being 2.55–2.87 and 0.80–0.88, respectively.
Note that typical values for the revPBE functional range from 2.34 to 2.94 and
from 0.76 to 1.10 for $g_{1}$ and $g_{2}$, respectively,[91] depending on the
simulation algorithm and dispersion correction. For the BLYP-based
simulations, the respective values cover an even broader range, 2.40–3.17 and
0.52–1.0.[91] It shows that the approximation introduced by the GM-NN model
does not significantly influence the resulting structural properties sampled
during a long MD simulation.

For oxygen–hydrogen ($g_{\ce{OH}}$) and hydrogen–hydrogen ($g_{\ce{HH}}$)
curves, the intensities of all peaks shown in Fig. 6 are slightly
overestimated in comparison with the experimental results. It can be
attributed to the lack of quantum nuclear effects (QNE) in solutions and is
out of the scope of this paper. We should mention that, here, no direct
comparison for the first peaks of the $g_{\ce{OH}}$ and $g_{\ce{HH}}$ curves
could be made since those are not given by the experiment used as a reference
in this work.[94] However, it is known that they are, typically, very
overstructured in classical MD simulations, and the inclusion of QNE is
necessary to get correct estimates.[1] Note that the $g_{\ce{OH}}$ curves for
revPBE-D3-bulk and revPBE-D3-cluster agree with the results obtained for the
oxygen–hydrogen RDF in Ref. 95.

####  III.2.2  Diffusion Coefficient

Another physical quantity we are interested in is the self-diffusion
coefficient of the water molecules. It can be calculated from a linear fit of
the mean square displacement (MSD) of the molecular centre of mass and
performing a fit using

$\langle\Delta x^{2}\rangle=2Dd\Delta t,$ (9)

where $D$ is the diffusion coefficient, $d=3$ the dimensionality of the
system, and $\Delta t$ the time interval for which the respective MSD was
measured. In this work, we used time intervals in the range of 1–10 ps.
Because of the finite-size effect on diffusion, an extrapolation to infinite
system size is required to allow for comparison to the experiment.[96] To get
the extrapolated diffusion coefficient $D\left(\infty\right)$ we add a
correction to the the finite-system diffusion coefficient $D\left(L\right)$
and the resulting expression reads[96]

$D\left(\infty\right)=D\left(L\right)+\frac{\xi}{6\pi\beta\eta L}.$ (10)

Here, $L$ is the length of the simulation box, $\beta=1/k_{\mathrm{B}}T$,
$\xi$ is given by the geometry of the simulation box ($\xi=2.837297$ for a
cubic box) and $\eta=0.8925\times 10^{-3}$ Pa s is the experimental shear
viscosity.

Figure 7: Water molecule self-diffusion coefficients from molecular dynamics
simulations at 300 K based on revPBE-D3-bulk, revPBE-D3-cluster,
BLYP-D3-cluster, and B3LYP-D3-cluster GM-NN potentials. Error bars are
obtained as the standard deviation between 10 estimates. The black line
corresponds to the experimental value while the shaded area denotes its error
bar.[97]

Fig. 7 shows the system size corrected diffusion coefficients obtained from
the MD simulations based on GM-NN interatomic potentials along with the
respective experimental value.[97] Similar to Section III.2.1 to get some
uncertainty estimate for the respective values of diffusion coefficients, we
split the whole 2.0 ns long trajectory into ten parts. Error bars for the ML-
based diffusion coefficients are obtained as the standard deviation between 10
estimates.

The system size corrected diffusion coefficients of $\left(1.82\pm
0.06\right)\times 10^{-9}$ m2 s-1, $\left(2.15\pm 0.07\right)\times 10^{-9}$
m2 s-1, $\left(1.39\pm 0.06\right)\times 10^{-9}$ m2 s-1, and $\left(1.43\pm
0.04\right)\times 10^{-9}$ m2 s-1 obtained employing revPBE-D3-bulk,
revPBE-D3-cluster, BLYP-D3-cluster, and B3LYP-D3-cluster-based interatomic
potentials, respectively. All these values are in fair agreement with
experimental value of $\left(2.41\pm 0.15\right)\times 10^{-9}$ m2 s-1.[97]
Note that the empirical force fields usually produce a rather larger diffusion
coefficient of $>2.6\times 10^{-9}$ m2 s-1,[98] while the DFT-based AIMD
simulations usually give a relatively smaller value of $<2.0\times 10^{-9}$ m2
s-1.[99, 95] More specifically, the ab-initio MD simulations result in self-
diffusion coefficients of merely 0.71$\times 10^{-9}$ to 1.81$\times 10^{-9}$
m2 s-1 and 1.6$\times 10^{-9}$ to 3.4$\times 10^{-9}$ m2 s-1 for BLYP and
revPBE, respectively, depending on the simulation algorithm and the dispersion
correction,[91] and are close to the values predicted using our potentials.

Using the GM-NN potential trained on local revPBE-D3-cluster data produces
diffusion coefficients close to the one obtained in Ref. 1 of $\left(2.22\pm
0.05\right)\times 10^{-9}$ m2 s-1, where AIMD simulations were performed with
the revPBE-D3 functional. Moreover, using the local models, we are able to
produce results close to the one presented in Ref. 62 ($2.0\times 10^{-9}$ m2
s-1) where the fragment-based coupled-cluster approach was employed.

As a final remark, we would like to note that besides the usage of machine-
learned interatomic potentials, other approximations made in this study may
influence the observed quantities, such as the employed canonical (NVT)
statistical ensemble for the system and Langevin thermostat for adjusting the
system’s temperature. The agreement of the computational results obtained in
this work with their experimental counterparts is mainly achieved only due to
the accurate atom-centred machine-learned potentials and their ability to
transfer from cluster to periodic systems.

###  III.3  Estimating equilibrium density of bulk liquid water

The ability of MLPs trained on cluster data to reproduce the correct
equilibrium density of bulk liquid water at ambient conditions is a very
useful criterion to assess the generalization ability of the former. For this
purpose, we run molecular dynamics (MD) simulations in the NPT (constant
pressure) statistical ensemble. Specifically, the isobaric-isothermal form of
the Nosé–Hoover dynamics has been employed[100, 101] as implemented in
ASE.[78] All production simulations have been performed at ambient conditions,
i.e., $T=300$ K and $p=1$ bar, over 1 ns after 200 ps long equilibration. A
time step of 0.5 fs has been used to integrate the equations of motion, while
the characteristic time scales of the thermostat and barostat were set to 1 ps
each. The simulation box was allowed to change independently along the three
Cartesian axes, $x$, $y$, and $z$. However, the angle between axes has been
fixed to 90 degrees. Finally, the density of the bulk liquid water has been
computed from the time-averaged (1 ns) volume of the simulation box.

The MD simulations for the MLP trained on the periodic revPBE-D3-bulk data set
resulted in an underestimated density, $0.86$ g cm-3. By contrast, the MLP
trained on local revPBE-D3-cluster data produced bulk liquid water with a much
more realistic density of $1.02$ g cm-3. The well-established TIP4P water
model provides a density of $1.00$ g cm-3 for comparison.[64] The difference
between the estimated densities employing cluster and the periodic models
might arise from various aspects of training data generation. Moreover, on par
with the observations from Section III.2, BLYP-D3-cluster and
B3LYP-D3-cluster-based NN potentials produce bulk liquid water with an
overestimated density of $1.12$ and $1.10$ g cm-3, respectively.

We have observed in Section III.2, in agreement with previous studies,[91, 1]
that the revPBE-D3 level is well suited to describe bulk liquid water.
However, the primary purpose of this section is not to assess the reference
DFT but rather to study the limitations of the cluster models in terms of
their generalisation ability. From the presented results, the excellent
transferability of cluster models to periodic bulk systems can be deduced.
However, further studies on systems different from bulk liquid water are
necessary to get a deeper insight into this intriguing feature of MLPs.

Under reasonably general conditions, the only difference between the cluster
and bulk models are the local environments of surface atoms and the atoms at
the periodic boundaries. The former lack bonding partners and, thus,
experience increased strain energy. This effect leads to atomic arrangements
at the surface being different from bulk ones. Note that they are physically
reasonable and do not bias the model to a nonphysical behaviour. When running
simulations of periodic systems with the cluster models, these atomic
arrangements are not probable and, thus, are not visited. However, they might
contribute to a better performance of cluster models when computing stress
tensors. It could explain the better performance of the revPBE-D3-cluster-
based model than the revPBE-D3-bulk-based model.

Another aspect that could influence the performance of MLPs trained on cluster
data are long-range electrostatic interactions, which are absent in our MLPs.
However, no important contributions of the long-range interactions could be
found, which is in agreement with previous work,[26] where no significant
difference in the accuracy between MLPs with and without explicit long-range
electrostatics for a cutoff of 6.35 Å. We have employed a cutoff of 6.5 Å.

##  IV  Conclusions

We have studied the capability of atom-centred machine-learned potentials
(MLPs) to transfer the information gained from local training data to periodic
systems. Such a transfer is essential for the construction of training data
sets for periodic systems from ab initio methods not feasible otherwise, e.g.,
hybrid functionals such as B3LYP. Simular transferability is used for non-
conventional force fields other than MLPs, e.g., the recent fragment-based
coupled-cluster approach.[62] For more examples of non-conventional force
fields, see Ref. 56 and references therein. Here, the transfer of structure-
property predictions of atom-centred MLPs has been studied on an example of
liquid water, while more elaborate applications can be feasible, such as
molten alkali chloride salts.[30] Note that no explicit electrostatics for our
water models has been employed since it is irrelevant for large enough cutoff
radii.[26]

Applying the machine-learned potentials trained on local data, we could obtain
oxygen–oxygen, oxygen–hydrogen, and hydrogen–hydrogen radial distribution
functions (RDF). We have found a good agreement with the experiment for the
oxygen–oxygen RDF, while for the others, somewhat overestimated intensities
were observed due to the lack of quantum nuclear effects (QNE). Moreover, we
obtained values reasonably close to the experimental data and on par with the
recent ab initio molecular dynamics studies by studying the self-diffusion
coefficients. Finally, by running molecular dynamics simulations at constant
pressure, the equilibrium densities for our cluster models could be estimated.
For the revPBE-D3-cluster-based model, we have found an excellent agreement
with the experiment, while its periodic counterpart provided somewhat
underestimated density.

In summary, we want to state the central finding of this work. We have shown
that it is feasible to transfer machine-learned interatomic potentials trained
on local data to periodic systems. It reduces the computational cost required
for the creation of training data using ab-initio methods more suited for
cluster models, e.g. hybrid functionals such as B3LYP. Moreover, the
possibility to transfer a local machine-learned potential to periodic systems
opens the door to predictions of bulk physical quantities using traditional
wave-function-based approaches, like coupled-cluster approach. However,
further development to reach such goals is necessary, which will be a subject
of our future investigations.

###### Acknowledgements.

We thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
for supporting this work by funding EXC 2075 - 390740016 under Germany’s
Excellence Strategy. We acknowledge the support by the Stuttgart Center for
Simulation Science (SimTech). The authors acknowledge support by the state of
Baden-Württemberg through the bwHPC consortium for providing computer and GPU
time. The authors thank the International Max Planck Research School for
Intelligent Systems (IMPRS-IS) for supporting D. H., while V. Z. acknowledges
the financial support received in the form of a PhD scholarship from the
Studienstiftung des Deutschen Volkes (German National Academic Foundation).

## References

  * Marsalek and Markland [2017] O. Marsalek and T. E. Markland, “Quantum Dynamics and Spectroscopy of Ab Initio Liquid Water: The Interplay of Nuclear and Electronic Quantum Effects,” J. Phys. Chem. Lett. 8, 1545–1551 (2017).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, “Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces,” Phys. Rev. Lett. 98, 146401 (2007).
  * Behler [2011] J. Behler, “Atom-centered symmetry functions for constructing high-dimensional neural network potentials,” J. Chem. Phys. 134, 074106 (2011).
  * Bartók  _et al._ [2010] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, “Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons,” Phys. Rev. Lett. 104, 136403 (2010).
  * Rupp  _et al._ [2012] M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, “Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning,” Phys. Rev. Lett. 108, 058301 (2012).
  * Bartók, Kondor, and Csányi [2013] A. P. Bartók, R. Kondor, and G. Csányi, “On representing chemical environments,” Phys. Rev. B 87, 184115 (2013).
  * von Lilienfeld  _et al._ [2015] O. A. von Lilienfeld, R. Ramakrishnan, M. Rupp, and A. Knoll, “Fourier series of atomic radial distribution functions: A molecular fingerprint for machine learning models of quantum chemical properties,” Int. J. Quantum Chem. 115, 1084–1093 (2015).
  * Shapeev [2016] A. V. Shapeev, “Moment Tensor Potentials: A Class of Systematically Improvable Interatomic Potentials,” Multiscale Model. Simul. 14, 1153–1173 (2016).
  * Khorshidi and Peterson [2016] A. Khorshidi and A. A. Peterson, “Amp: A modular approach to machine learning in atomistic simulations,” Comput. Phys. Commun. 207, 310 – 324 (2016).
  * Schütt  _et al._ [2017] K. Schütt, P.-J. Kindermans, H. E. Sauceda Felix, S. Chmiela, A. Tkatchenko, and K.-R. Müller, “Schnet: A continuous-filter convolutional neural network for modeling quantum interactions,” in  _NeurIPS 30_ , edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017) pp. 991–1001.
  * Faber  _et al._ [2018] F. A. Faber, A. S. Christensen, B. Huang, and O. A. von Lilienfeld, “Alchemical and structural distribution based representation for universal quantum machine learning,” J. Chem. Phys. 148, 241717 (2018).
  * Kocer, Mason, and Erturk [2019] E. Kocer, J. K. Mason, and H. Erturk, “A novel approach to describe chemical environments in high-dimensional neural network potentials,” J. Chem. Phys. 150, 154102 (2019).
  * Zhang, Hu, and Jiang [2019] Y. Zhang, C. Hu, and B. Jiang, “Embedded Atom Neural Network Potentials: Efficient and Accurate Machine Learning with a Physically Inspired Representation,” J. Phys. Chem. Lett. 10, 4962–4967 (2019).
  * Christensen  _et al._ [2020] A. S. Christensen, L. A. Bratholm, F. A. Faber, and O. von Lilienfeld, “FCHL revisited: Faster and more accurate quantum machine learning,” J. Chem. Phys. 152, 44107 (2020).
  * Zaverkin and Kästner [2020] V. Zaverkin and J. Kästner, “Gaussian Moments as Physically Inspired Molecular Descriptors for Accurate and Scalable Machine Learning Potentials,” J. Chem. Theory Comput. 16, 5410–5421 (2020).
  * Schütt, Unke, and Gastegger [2021] K. T. Schütt, O. T. Unke, and M. Gastegger, “Equivariant message passing for the prediction of tensorial properties and molecular spectra,” ArXiv abs/2102.03150 (2021).
  * Artrith and Urban [2016] N. Artrith and A. Urban, “An implementation of artificial neural-network potentials for atomistic materials simulations: Performance for TiO2,” Comput. Mater. Sci. 114, 135–150 (2016).
  * Chmiela  _et al._ [2017] S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, “Machine learning of accurate energy-conserving molecular force fields,” Sci. Adv. 3 (2017).
  * Gubaev, Podryabinkin, and Shapeev [2018] K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, “Machine learning of molecular properties: Locality and active learning,” J. Chem. Phys. 148, 241727 (2018).
  * Lubbers, Smith, and Barros [2018] N. Lubbers, J. S. Smith, and K. Barros, “Hierarchical modeling of molecular energies using a deep neural network,” J. Chem. Phys. 148, 241715 (2018).
  * Yao  _et al._ [2018] K. Yao, J. E. Herr, D. W. Toth, R. Mckintyre, and J. Parkhill, “The TensorMol-0.1 model chemistry: a neural network augmented with long-range physics,” Chem. Sci. 9, 2261–2269 (2018).
  * Unke and Meuwly [2019] O. T. Unke and M. Meuwly, “PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments, and Partial Charges,” J. Chem. Theory Comput. 15, 3678–3693 (2019).
  * Cooper  _et al._ [2020a] A. M. Cooper, J. Kästner, A. Urban, and N. Artrith, “Efficient training of ANN potentials by including atomic forces via Taylor expansion and application to water and a transition-metal oxide,” npj Comput. Mater. 6 (2020a).
  * Zaverkin  _et al._ [2021] V. Zaverkin, D. Holzmüller, I. Steinwart, and J. Kästner, “Fast and Sample-Efficient Interatomic Neural Network Potentials for Molecules and Materials Based on Gaussian Moments,” J. Chem. Theory Comput. 17, 6658 – 6670 (2021).
  * Hornik [1991] K. Hornik, “Approximation capabilities of multilayer feedforward networks,” Neural Netw. 4, 251 – 257 (1991).
  * Morawietz  _et al._ [2016] T. Morawietz, A. Singraber, C. Dellago, and J. Behler, “How van der Waals interactions determine the unique properties of water,” PNAS 113, 8368–8373 (2016).
  * Brickel  _et al._ [2019] S. Brickel, A. K. Das, O. T. Unke, H. T. Turan, and M. Meuwly, “Reactive molecular dynamics for the [Cl–CH3–Br]- reaction in the gas phase and in solution: a comparative study using empirical and neural network force fields,” Electron. Struct. 1, 024002 (2019).
  * Käser, Unke, and Meuwly [2020] S. Käser, O. T. Unke, and M. Meuwly, “Isomerization and decomposition reactions of acetaldehyde relevant to atmospheric processes from dynamics simulations on neural network-based potential energy surfaces,” J. Chem. Phys. 152, 214304 (2020).
  * Molpeceres, Zaverkin, and Kästner [2020a] G. Molpeceres, V. Zaverkin, and J. Kästner, “Neural-network assisted study of nitrogen atom dynamics on amorphous solid water – i. adsorption and desorption,” Mon. Not. R. Astron. Soc. 499, 1373–1384 (2020a).
  * Tovey  _et al._ [2020] S. Tovey, A. Narayanan Krishnamoorthy, G. Sivaraman, J. Guo, C. Benmore, A. Heuer, and C. Holm, “DFT Accurate Interatomic Potential for Molten NaCl from Machine Learning,” J. Phys. Chem. C 124, 25760–25768 (2020).
  * Korotaev and Shapeev [2020] P. Korotaev and A. Shapeev, “Lattice dynamics of ${\mathrm{Yb}}_{x}{\mathrm{Co}}_{4}{\mathrm{Sb}}_{12}$ skutterudite by machine-learning interatomic potentials: Effect of filler concentration and disorder,” Phys. Rev. B 102, 184305 (2020).
  * Calegari Andrade  _et al._ [2020] M. F. Calegari Andrade, H.-Y. Ko, L. Zhang, R. Car, and A. Selloni, “Free energy of proton transfer at the water–TiO2 interface from ab initio deep potential molecular dynamics,” Chem. Sci. 11, 2335–2341 (2020).
  * Molpeceres, G.  _et al._ [2021] Molpeceres, G., Zaverkin, V., Watanabe, N., and Kästner, J., “Binding energies and sticking coefficients of H2 on crystalline and amorphous CO ice,” A&A 648, A84 (2021).
  * Zhang  _et al._ [2021] L. Zhang, H. Wang, R. Car, and W. E, “Phase Diagram of a Deep Potential Water Model,” Phys. Rev. Lett. 126, 236001 (2021).
  * Blank  _et al._ [1995] T. B. Blank, S. D. Brown, A. W. Calhoun, and D. J. Doren, “Neural network models of potential energy surfaces,” J. Chem. Phys. 103, 4129–4137 (1995).
  * Lorenz, Groß, and Scheffler [2004] S. Lorenz, A. Groß, and M. Scheffler, “Representing high-dimensional potential-energy surfaces for reactions at surfaces by neural networks,” Chem. Phys. Lett. 395, 210–215 (2004).
  * Becke [1993] A. D. Becke, “Density-functional thermochemistry. III. The role of exact exchange,” J. Chem. Phys. 98, 5648–5652 (1993).
  * Stephens  _et al._ [1994] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, “Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields,” J. Phys. Chem. 98, 11623–11627 (1994).
  * Ghasemi  _et al._ [2015] S. A. Ghasemi, A. Hofstetter, S. Saha, and S. Goedecker, “Interatomic potentials for ionic systems with density functional accuracy based on charge densities obtained by a neural network,” Phys. Rev. B 92, 045131 (2015).
  * Faraji  _et al._ [2017] S. Faraji, S. A. Ghasemi, S. Rostami, R. Rasoulkhani, B. Schaefer, S. Goedecker, and M. Amsler, “High accuracy and transferability of a neural network potential through charge equilibration for calcium fluoride,” Phys. Rev. B 95, 104105 (2017).
  * Rappe and Goddard [1991] A. K. Rappe and W. A. Goddard, “Charge equilibration for molecular dynamics simulations,” J. Phys. Chem. 95, 3358–3363 (1991).
  * Wilmer, Kim, and Snurr [2012] C. E. Wilmer, K. C. Kim, and R. Q. Snurr, “An Extended Charge Equilibration Method,” J. Phys. Chem. Lett. 3, 2506–2511 (2012).
  * Cheng  _et al._ [2014] Y.-T. Cheng, T.-R. Shan, T. Liang, R. K. Behera, S. R. Phillpot, and S. B. Sinnott, “A charge optimized many-body (comb) potential for titanium and titania,” J. Condens. Matter Phys. 26, 315007 (2014).
  * Gastegger  _et al._ [2016] M. Gastegger, C. Kauffmann, J. Behler, and P. Marquetand, “Comparing the accuracy of high-dimensional neural network potentials and the systematic molecular fragmentation method: A benchmark study for all-trans alkanes,” J. Chem. Phys. 144, 194110 (2016).
  * Monserrat  _et al._ [2020] B. Monserrat, J. G. Brandenburg, E. A. Engel, and B. Cheng, “Liquid water contains the building blocks of diverse ice phases,” Nat. Commun. 11, 5757 (2020).
  * Rowe  _et al._ [2020] P. Rowe, V. L. Deringer, P. Gasparotto, G. Csányi, and A. Michaelides, “An accurate and transferable machine learning potential for carbon,” J. Chem. Phys. 153, 34702 (2020).
  * Westermayr and Marquetand [2020] J. Westermayr and P. Marquetand, “Deep learning for UV absorption spectra with SchNarc: First steps toward transferability in chemical compound space,” J. Chem. Phys. 153, 154112 (2020).
  * Schran, Brieuc, and Marx [2021] C. Schran, F. Brieuc, and D. Marx, “Transferability of machine learning potentials: Protonated water neural network potential applied to the protonated water hexamer,” J. Chem. Phys. 154, 51101 (2021).
  * Hajibabaei  _et al._ [2021] A. Hajibabaei, M. Ha, S. Pourasad, J. Kim, and K. S. Kim, “Machine Learning of First-Principles Force-Fields for Alkane and Polyene Hydrocarbons,” J. Phys. Chem. A 125, 9414–9420 (2021).
  * Collins and Deev [2006] M. A. Collins and V. A. Deev, “Accuracy and efficiency of electronic energies from systematic molecular fragmentation,” J. Chem. Phys. 125, 104104 (2006).
  * Netzloff and Collins [2007] H. M. Netzloff and M. A. Collins, “Ab initio energies of nonconducting crystals by systematic fragmentation,” J. Chem. Phys. 127, 134113 (2007).
  * Collins [2012] M. A. Collins, “Systematic fragmentation of large molecules by annihilation,” Phys. Chem. Chem. Phys. 14, 7744–7751 (2012).
  * Gordon  _et al._ [2012] M. S. Gordon, D. G. Fedorov, S. R. Pruitt, and L. V. Slipchenko, “Fragmentation Methods: A Route to Accurate Calculations on Large Systems,” Chem. Rev. 112, 632–672 (2012).
  * Collins and Bettens [2015] M. A. Collins and R. P. A. Bettens, “Energy-Based Molecular Fragmentation Methods,” Chem. Rev. 115, 5607–5642 (2015).
  * Tröster  _et al._ [2013] P. Tröster, K. Lorenzen, M. Schwörer, and P. Tavan, “Polarizable Water Models from Mixed Computational and Empirical Optimization,” J. Phys. Chem. B 117, 9486–9500 (2013).
  * Koner  _et al._ [2020] D. Koner, S. M. Salehi, P. Mondal, and M. Meuwly, “Non-conventional force fields for applications in spectroscopy and chemical reaction dynamics,” J. Chem. Phys. 153, 10901 (2020).
  * Huang, Braams, and Bowman [2006] X. Huang, B. J. Braams, and J. M. Bowman, “Ab Initio Potential Energy and Dipole Moment Surfaces of (H2O)2,” J. Phys. Chem. A 110, 445–451 (2006).
  * Bukowski  _et al._ [2007] R. Bukowski, K. Szalewicz, G. C. Groenenboom, and A. van der Avoird, “Predictions of the Properties of Water from First Principles,” Science 315, 1249–1252 (2007).
  * Babin, Medders, and Paesani [2012] V. Babin, G. R. Medders, and F. Paesani, “Toward a Universal Water Model: First Principles Simulations from the Dimer to the Liquid Phase,” J. Phys. Chem. Lett. 3, 3765–3769 (2012).
  * Babin, Leforestier, and Paesani [2013] V. Babin, C. Leforestier, and F. Paesani, “Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient,” J. Chem. Theory Comput. 9, 5395–5403 (2013).
  * He  _et al._ [2014] X. He, T. Zhu, X. Wang, J. Liu, and J. Z. H. Zhang, “Fragment Quantum Mechanical Calculation of Proteins and Its Applications,” Acc. Chem. Res. 47, 2748–2757 (2014).
  * Liu  _et al._ [2018] J. Liu, X. He, J. Z. H. Zhang, and L.-W. Qi, “Hydrogen-bond structure dynamics in bulk water: insights from ab initio simulations with coupled cluster theory,” Chem. Sci. 9, 2065–2073 (2018).
  * Liu, Zhang, and He [2018] J. Liu, J. Z. H. Zhang, and X. He, “Probing the Ion-Specific Effects at the Water/Air Interface and Water-Mediated Ion Pairing in Sodium Halide Solution with Ab Initio Molecular Dynamics,” J. Phys. Chem. B 122, 10202–10209 (2018).
  * Jorgensen and Jenson [1998] W. L. Jorgensen and C. Jenson, “Temperature dependence of TIP3P, SPC, and TIP4P water from NPT Monte Carlo simulations: Seeking temperatures of maximum density,” J. Comput. Chem. 19, 1179–1186 (1998).
  * Blöchl [1994] P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, 17953–17979 (1994).
  * Zhang and Yang [1998] Y. Zhang and W. Yang, “Comment on “Generalized Gradient Approximation Made Simple”,” Phys. Rev. Lett. 80, 890–890 (1998).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu,” J. Chem. Phys. 132, 154104 (2010).
  * Cooper  _et al._ [2020b] A. Cooper, J. Kästner, A. Urban, and N. Artrith, “Efficient Training of ANN Potentials by Including Atomic Forces via Taylor Expansion and Application to Water and a Transition-Metal Oxide,” (2020b), 10.24435/materialscloud:2020.0037/v1 .
  * Molpeceres, Zaverkin, and Kästner [2020b] G. Molpeceres, V. Zaverkin, and J. Kästner, “N-ASW: Molecular Dynamics Data [Data set],” (2020b), http://doi.org/10.5281/zenodo.4013889 .
  * Weigend and Ahlrichs [2005] F. Weigend and R. Ahlrichs, “Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy,” Phys. Chem. Chem. Phys. 7, 3297–3305 (2005).
  * Weigend [2006] F. Weigend, “Accurate Coulomb-fitting basis sets for H to Rn,” Phys. Chem. Chem. Phys. 8, 1057–1065 (2006).
  * Furche  _et al._ [2014] F. Furche, R. Ahlrichs, C. Hättig, W. Klopper, M. Sierka, and F. Weigend, “Turbomole,” WIREs Comput. Mol. Sci. 4, 91–100 (2014).
  * Metz  _et al._ [2014] S. Metz, J. Kästner, A. A. Sokol, T. W. Keal, and P. Sherwood, “ChemShell—a modular software package for QM/MM simulations,” WIREs Comput. Mol. Sci. 4, 101–110 (2014).
  * Sherwood  _et al._ [2003] P. Sherwood, A. H. de Vries, M. F. Guest, G. Schreckenbach, C. A. Catlow, S. A. French, A. A. Sokol, S. T. Bromley, W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, J. Kendrick, S. C. Rogers, J. Casci, M. Watson, F. King, E. Karlsen, M. Sjøvoll, A. Fahmi, A. Schäfer, and C. Lennartz, “QUASI: A general purpose implementation of the QM/MM approach and its application to problems in catalysis,” J. Mol. Struc.-THEOCHEM 632, 1 – 28 (2003).
  * Lee, Yang, and Parr [1988] C. Lee, W. Yang, and R. G. Parr, “Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density,” Phys. Rev. B 37, 785–789 (1988).
  * Zaverkin  _et al._ [2022] V. Zaverkin, D. Holzmüller, R. Schuldt, and J. Kästner, “Predicting Properties of Periodic Systems from Cluster Data: A Case Study of Liquid Water [Data set],” (2022), https://doi.org/10.5281/zenodo.5817133 .
  * Kingma and Ba [2015] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” CoRR abs/1412.6980 (2015).
  * Larsen  _et al._ [2017] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, “The atomic simulation environment—a Python library for working with atoms,” J. Phys. Condens. Matter 29, 273002 (2017).
  * Gastegger, Behler, and Marquetand [2017] M. Gastegger, J. Behler, and P. Marquetand, “Machine learning molecular dynamics for the simulation of infrared spectra,” Chem. Sci. 8, 6924–6935 (2017).
  * Smith  _et al._ [2018] J. S. Smith, B. Nebgen, N. Lubbers, O. Isayev, and A. E. Roitberg, “Less is more: Sampling chemical space with active learning,” J. Chem. Phys. 148, 241733 (2018).
  * Zhang  _et al._ [2019] L. Zhang, D.-Y. Lin, H. Wang, R. Car, and W. E, “Active learning of uniformly accurate interatomic potentials for materials simulation,” Phys. Rev. Materials 3, 023804 (2019).
  * Janet  _et al._ [2019] J. P. Janet, C. Duan, T. Yang, A. Nandy, and H. J. Kulik, “A quantitative uncertainty metric controls error in neural network-driven chemical discovery,” Chem. Sci. 10, 7913–7922 (2019).
  * Schran, Behler, and Marx [2020] C. Schran, J. Behler, and D. Marx, “Automated Fitting of Neural Network Potentials at Coupled Cluster Accuracy: Protonated Water Clusters as Testing Ground,” J. Chem. Theory Comput. 16, 88–99 (2020).
  * Schran, Brezina, and Marsalek [2020] C. Schran, K. Brezina, and O. Marsalek, “Committee neural network potentials control generalization errors and enable active learning,” J. Chem. Phys. 153, 104105 (2020).
  * Imbalzano  _et al._ [2021] G. Imbalzano, Y. Zhuang, V. Kapil, K. Rossi, E. A. Engel, F. Grasselli, and M. Ceriotti, “Uncertainty estimation for molecular dynamics and sampling,” J. Chem. Phys. 154, 74102 (2021).
  * Settles [2009] B. Settles, “Active Learning Literature Survey,” Computer Sciences Technical Report 1648 (University of Wisconsin–Madison, 2009).
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “ TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems,”  (2015), software available from tensorflow.org.
  * Thompson, Plimpton, and Mattson [2009] A. P. Thompson, S. J. Plimpton, and W. Mattson, “General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions,” J. Chem. Phys. 131, 154107 (2009).
  * Zhang  _et al._ [2018] L. Zhang, J. Han, H. Wang, R. Car, and W. E, “Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics,” Phys. Rev. Lett. 120, 143001 (2018).
  * Shao  _et al._ [2020] Y. Shao, M. Hellström, P. D. Mitev, L. Knijff, and C. Zhang, “PiNN: A Python Library for Building Atomic Neural Networks of Molecules and Materials,” J. Chem. Inf. Model. 60, 1184–1193 (2020).
  * Gillan, Alfè, and Michaelides [2016] M. J. Gillan, D. Alfè, and A. Michaelides, “Perspective: How good is DFT for water?” J. Chem. Phys. 144, 130901 (2016).
  * Zaverkin and Kästner [2021] V. Zaverkin and J. Kästner, “Exploration of transferable and uniformly accurate neural network interatomic potentials using optimal experimental design,” Mach. Learn.: Sci. Technol. 2, 035009 (2021).
  * Skinner  _et al._ [2013] L. B. Skinner, C. Huang, D. Schlesinger, L. G. M. Pettersson, A. Nilsson, and C. J. Benmore, “Benchmark oxygen-oxygen pair-distribution function of ambient water from x-ray diffraction measurements with a wide Q-range,” J. Chem. Phys. 138, 74506 (2013).
  * Soper [2013] A. K. Soper, “The Radial Distribution Functions of Water as Derived from Radiation Total Scattering Experiments: Is There Anything We Can Say for Sure?” ISRN Phys. Chem. 2013, 279463 (2013).
  * Ruiz Pestana  _et al._ [2017] L. Ruiz Pestana, N. Mardirossian, M. Head-Gordon, and T. Head-Gordon, “Ab initio molecular dynamics simulations of liquid water using high quality meta-GGA functionals,” Chem. Sci. 8, 3554–3565 (2017).
  * Yeh and Hummer [2004] I.-C. Yeh and G. Hummer, “System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions,” J. Phys. Chem. B 108, 15873–15879 (2004).
  * Holz, Heil, and Sacco [2000] M. Holz, S. R. Heil, and A. Sacco, “Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1H NMR PFG measurements,” Phys. Chem. Chem. Phys. 2, 4740–4742 (2000).
  * Price and Brooks [2004] D. J. Price and C. L. Brooks, “A modified TIP3P water potential for simulation with Ewald summation,” J. Chem. Phys. 121, 10096–10103 (2004).
  * DiStasio  _et al._ [2014] R. A. DiStasio, B. Santra, Z. Li, X. Wu, and R. Car, “The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water,” J. Chem. Phys. 141, 84502 (2014).
  * Melchionna, Ciccotti, and Holian [1993] S. Melchionna, G. Ciccotti, and B. L. Holian, “Hoover NPT dynamics for systems varying in shape and size,” Mol. Phys. 78, 533–544 (1993).
  * Melchionna [2000] S. Melchionna, “Constrained systems and statistical distribution,” Phys. Rev. E 61, 6165–6170 (2000).


# Rabi oscillations at the exceptional point in anti-parity-time symmetric
diffusive systems

Gabriel González Cátedra CONACYT, Universidad Autónoma de San Luis Potosí, San
Luis Potosí, 78000 MEXICO. E-mail<EMAIL_ADDRESS>
###### Abstract

The motivation for this theoretical paper comes from recent experiments of a
heat transfer system of two thermally coupled rings rotating in opposite
directions with equal angular velocities that present anti-parity-time (APT)
symmetry. The theoretical model predicted a rest-to-motion temperature
distribution phase transition during the symmetry breaking for a particular
rotation speed. In this work we show that the system exhibits a parity-time
($\mathcal{PT}$) phase transition at the exceptional point in which
eigenvalues and eigenvectors of the corresponding non-Hermitian Hamiltonian
coalesce. We analytically solve the heat diffusive system at the exceptional
point and show that one can pass through the phase transition that separates
the unbroken and broken phases by changing the radii of the rings. In the case
of unbroken $\mathcal{PT}$ symmetry the temperature profiles exhibit damped
Rabi oscillations at the exceptional point. Our results unveils the behavior
of the system at the exceptional point in heat diffusive systems.

A closed or conservative system evolves according to a Hermitian Hamiltonian
in contrast with open or non conservative systems which are described by non-
Hermitian Hamiltonians. There are a special class of non-Hermitian systems in
which the energy exchange between the system and the environment is balanced.
The entire balanced system exhibits a symmetry called $\mathcal{PT}$ symmetry
where the symbol $\mathcal{P}$ stands for parity and interchanges the gain and
loss components of the total system and $\mathcal{T}$ represents the operation
of time reversal and has the effect of turning a system with loss into a
system with gain and viceversa.bender ; bender1 ; bender2
Non-Hermitian $\mathcal{PT}$ symmetric systems can exhibit a rich and
unexpected behavior and have broad applications in classical and quantum
physics.bender3 ; bender4 ; bender5 ; bender6  $\mathcal{PT}$ symmetric
systems have been intensively studied in optics in which many intriguing
phenomena haven been experimentally confirmed and has led to the development
of new ways of controlling light propagation.guo ; ruter ; chang ; wimmer
Recently, anti-$\mathcal{PT}$ (APT) symmetric systems have attracted a lot of
attention because they exhibit noteworthy effects different from the
$\mathcal{PT}$ counterpart. An APT symmetric Hamiltonian can be defined in
terms of a $\mathcal{PT}$ symmetric Hamiltonian by $H^{(APT)}=\pm iH^{(PT)}$,
but physically it is really difficult to implement it in the laboratory since
it requires the coupling between the two subsystems to be a purely imaginary
value, in contrast with the $\mathcal{PT}$ systems which requires a real
coupling. Anti-PT symmetry has been demonstrated by using dissipatively
coupled atomic beams,peng  cold atoms,jiang  electrical circuits,choi  and
optical devices.zhang ; lic ; zhao  These breakthroughs have initiated the
field of exploring unique APT effects. More recently, Li et al reported the
experimental realization of an APT symmetric diffusive system in Ref.X . The
system investigated in Ref.X  is depicted in Fig. 1 and consists of two
identical solid rings with inner and outer radius given by $R$ and $R+\delta
R$, respectively. The thickness is $b$. The upper ring is rotating with
angular velocity $\omega_{1}$, while the lower ring is rotating with angular
velocity $\omega_{2}=-\omega_{1}$. There is an interface of thickness $d$ and
thermal conductivity $k_{i}$ between the two rings. The temperature
distribution along the inner edges of the upper and lower rings is given by
the following diffusion coupled partial differential equations

$\displaystyle\frac{\partial T_{1}}{\partial
t}=D\frac{\partial^{2}T_{1}}{\partial x^{2}}-v\frac{\partial T_{1}}{\partial
x}+h_{c}(T_{2}-T_{1})$ $\displaystyle\frac{\partial T_{2}}{\partial
t}=D\frac{\partial^{2}T_{2}}{\partial x^{2}}+v\frac{\partial T_{2}}{\partial
x}+h_{c}(T_{1}-T_{2})$ (1)

|
---|---
(a) | (b)
Figure 1: The figure (a) shows two identical rotating rings with equal but
opposite angular velocities joined together by a stationary intermediate layer
and (b) the imaginary and real parts of the eigenfrequencies as a function of
the tangential velocity where the dotted line represents the exceptional point
$v_{EP}=h_{c}/\kappa$.

where $x$ is the coordinate along each edge, $D=k/\rho c$ is the diffusivity,
$v$ is the tangential velocity in the inner edge of the rings,$h_{c}=h/\rho
cb$ is the rate of heat exchange coupling, $\rho$ is the density, $c$ is the
heat capacity and $h=k_{i}/d$ is a coefficient that represent the heat
exchange between the two rings. Using plane wave solutions, i.e.
$T_{i}=A_{i}e^{i(\kappa x-\omega t)}$, the system given in Eq. (1) can be cast
into an APT symmetric Hamiltonian given by

$H^{(APT)}=\begin{pmatrix}-i(\kappa^{2}D+h_{c})+\kappa v&ih_{c}\\\
ih_{c}&-i(\kappa^{2}D+h_{c})-\kappa v\end{pmatrix}$ (2)

where $\kappa$ is the wave number and $\omega$ are the eigenvalues of the APT
Hamiltonian which are given by

$\omega_{\pm}=-i\left[(\kappa^{2}D+h_{c})\pm\sqrt{h_{c}^{2}-\kappa^{2}v^{2}}\right].$
(3)

The exceptional point where the two eigenvectors coalesce is when
$v_{EP}^{2}=h_{c}^{2}/\kappa^{2}$, i.e. $\omega_{+}=\omega_{-}$. The sudden
collapse of the eigenvectors and eigenvalues at the exceptional point leads to
an abrupt reduction in dimensionality. Many of the interesting properties of
non Hermitian systems are found at or close to the exceptional point which
have led to many novel and exotic phenomena. Exceptional points are currently
the subject of many interesting and counter-intuitive phenomena associated
with them such as topological mode switching,topo1 ; topo2  reflection and
transmission,refle1 ; refle2 ; refle3  instrinsic single-mode lasingintri1 ;
intri2  and coherent perfect absorption.coher
In this work we study the APT symmetric diffusive system given by Eq. (1) when
$v=v_{EP}$ and show that the system behaves as a pair of coupled linear
oscillators one with gain and the other one with loss. The noteworthy feature
of the exceptional point $v_{EP}$ is that it exhibits damped Rabi oscillations
in the unbroken $\mathcal{PT}$ phase transition that depends on the radii of
the rotating rings. We obtain the analytical temperature distribution of each
ring at the exceptional point and obtain the conditions that have to be
fulfill in order for the system to be in equilibrium. We start our
investigation by making the following variables change in Eq. (1):
$\tau=h_{c}t$, $z=\sqrt{h_{c}(\lambda-1)/D}x$ where $\lambda>1$ is an
auxiliary constant to be determined and $\Delta T=T_{i}-T_{0}$ where $T_{0}$
is a reference temperature. Rewriting Eq. (1) in terms of the new variables we
have

$\displaystyle h_{c}\frac{\partial\Delta
T_{1}}{\partial\tau}=h_{c}(\lambda-1)\frac{\partial^{2}\Delta T_{1}}{\partial
z^{2}}-v_{EP}\sqrt{\frac{h_{c}(\lambda-1)}{D}}\frac{\partial\Delta
T_{1}}{\partial z}+h_{c}(\Delta T_{2}-\Delta T_{1})$ $\displaystyle
h_{c}\frac{\partial\Delta
T_{2}}{\partial\tau}=h_{c}(\lambda-1)\frac{\partial^{2}\Delta T_{2}}{\partial
x^{2}}+v_{EP}\sqrt{\frac{h_{c}(\lambda-1)}{D}}\frac{\partial\Delta
T_{2}}{\partial z}+h_{c}(\Delta T_{1}-\Delta T_{2})$ (4)

Looking for solutions of the form $\Delta T_{i}=e^{-\lambda\tau}f_{i}(z)$ in
Eq. (4) we end up with the following system of coupled ordinary differential
equations

$\displaystyle\frac{d^{2}f_{1}}{dz^{2}}-\frac{v_{EP}}{\sqrt{Dh_{c}(\lambda-1)}}\frac{df_{1}}{dz}+f_{1}+\frac{1}{\lambda-1}f_{2}=0$
$\displaystyle\frac{d^{2}f_{2}}{dz^{2}}+\frac{v_{EP}}{\sqrt{Dh_{c}(\lambda-1)}}\frac{df_{2}}{dz}+f_{2}+\frac{1}{\lambda-1}f_{1}=0.$
(5)

Inspection of Eqs.(5) reveals that they are invariant under combined parity,
i.e. $f_{1}\leftrightarrow f_{2}$, and time reversal $t\rightarrow-t$
transformation. To solve the system of equations analytically we first
differentiate one of the equations and then use the other equation to
eliminate $f_{2}$ in order to get the following fourth order differential
equation

$\bigg{(}\frac{d^{4}}{dz^{4}}+\left(2-\frac{\epsilon
v^{2}_{EP}}{Dh_{c}}\right)\frac{d^{2}}{dz^{2}}+(1-\epsilon^{2})\bigg{)}f_{1}(z)=0.$
(6)

where $\epsilon=1/(\lambda-1)$. By assuming a solution of the form
$f_{1}(z)\propto\cosh(\chi z)$ for Eq. (6) we get the following condition over
$\chi$:

$\chi^{4}+(2-a^{2})\chi^{2}+(1-\epsilon^{2})=0,$ (7)

where $a^{2}=\epsilon v_{EP}^{2}/Dh_{c}$. The solution of Eq. (7) is given by

$\chi^{2}=\frac{1}{2}\left(a^{2}-2\pm\sqrt{a^{4}-4a^{2}+4\epsilon^{2}}\right).$
(8)

In order to have an oscillatory behavior we must demand that $\chi^{2}<0$,
which implies that

  * (i)

$a^{4}-4a^{2}+4\epsilon^{2}>0$

  * (ii)

$a^{2}-2+\sqrt{a^{4}-4a^{2}+4\epsilon^{2}}<0$.

Condition (ii) gives $\epsilon<1$ and condition (i) gives

$a<a_{crit}=2\left(1-\sqrt{1-\epsilon^{2}}\right).$ (9)

If $\epsilon<1$ and $a<a_{crit}$ we get the following oscillatory solution for
$f_{1}$

$f_{1}(z)=A_{1}\cos\left(\chi_{1}z\right)+B_{1}\cos\left(\chi_{2}z\right),$
(10)

where $\chi_{1,2}=\sqrt{|\chi_{\pm}^{2}|}$ and $A_{1}$ and $B_{1}$ are
constants to be determined. In order to obtain the value of $\epsilon$ we must
consider the periodicty of $f_{1}(z)$, i.e. $f(0)=f(2\pi
R\sqrt{h_{c}/D\epsilon})$, which gives us the following conditions

$\chi_{1,2}R\sqrt{\frac{h_{c}}{D\epsilon}}=n,$ (11)

where $n=\pm 1,\pm 2,\ldots$. Solving Eq. (11) we get the following value for
$\epsilon$

$\epsilon=\frac{h_{c}R^{2}}{Dn^{2}}.$ (12)

Using the fact that $\epsilon=1/(\lambda-1)$ we get the following value for
$\lambda$

$\lambda=1+\frac{Dn^{2}}{h_{c}R^{2}}$ (13)

Equation (13) is in agreement with Eq. (3) when
$v_{EP}^{2}=h_{c}^{2}/\kappa^{2}$ and $k=n/R$. Interestingly, Eq. (13) is
valid only when conditions (i) and (ii) are fulfilled.
Once we know the value of $\epsilon$ we can substitute in $a^{2}=\epsilon
v_{EP}^{2}/Dh_{c}$ in order to get $a=\epsilon$, substituting this value into
Eq. (9) we have the conditions that have to be satisfied in order to have
unbroken-$\mathcal{PT}$ symmetry at the exceptional point which are
$0<\epsilon<1$ and $0<\epsilon<2(1-\sqrt{1-\epsilon^{2}})$, which gives us the
following solution

$\frac{4}{5}<\epsilon<1$ (14)

Equation (14) is the main result of this study which states that two phase
transitions take place at the exceptional point and depends only on the radii
of the rotating rings.
Substituting $a=\epsilon$ in Eq. (8) we get $\chi_{+}^{2}=\epsilon^{2}-1$ and
$\chi_{-}^{2}=-1$, therefore
$f_{1}(z)=A_{1}\cos(\sqrt{1-\epsilon^{2}}z)+B_{1}\cos(z)$ which means we have
to choose $A_{1}=0$ in order to fulfill the periodicity condition.
Substituting $f_{1}$ into Eq. (5) we obtain the following ordinary
differential equation for $f_{2}$:

$\frac{d^{2}f_{2}}{dz^{2}}+\epsilon\frac{df_{2}}{dz}+f_{2}=-\epsilon
B_{1}\cos(z).$ (15)

The general solution for Eq. (15) is given by

$f_{2}(z)=A_{2}e^{-z/2\epsilon}\cos\left(\sqrt{1-(\epsilon/2)^{2}}z+\phi\right)-B_{1}\sin(z),$
(16)

where $A_{2}$, $B_{1}$, $\phi$ and $n$ are constants to be determined by the
initial conditions.
If we impose the following initial conditions over the temperature profiles in
the rings

$T_{1}(x,0)=T_{2}(x,0)=T_{0}+A\cos(x/R)$ (17)

we need to choose $n=1$ and $B_{1}=A$ in order to get

$f_{1}(x)=A\cos(x/R)$ (18)

and

$f_{2}(x)=Ae^{-Dx/2R^{3}h_{c}}\sec(\phi)\cos\left(\alpha\frac{x}{R}+\phi\right)-A\sin(x/R)$
(19)

where $\alpha=\sqrt{1-(h_{c}R^{2}/2D)^{2}}$ and

$\phi=\arctan\left[\cot(2\pi\alpha)-\csc(2\pi\alpha)e^{\pi
D/h_{c}R^{2}}\right].$ (20)

The solution given in Eq. (18) means that $\Delta T_{1}(x,0)=f_{1}(x)$,
therefore the temperature distribution in the first ring will not change in
position but will only decay on time, in contrast with the solution given in
Eq. (19) which is different from the initial condition, therefore the
temperature profile will change in position and decay on time. In Fig. (2) we
show the Rabi oscillations as a function of position for different values of
$\epsilon$.

|
---|---
(a) | (b)
Figure 2: Numerical solution of the coupled equations given in Eq. (5). The
graph shows the Rabi oscillations for (a) $\epsilon=0.86$ and (b)
$\epsilon=0.96$, respectively. The Rabi oscillations are harder to see near
the upper end of the unbroken $\mathcal{PT}$ symmetry.

If we impose the following new conditions over the temperature profiles in the
rings

$T_{1}(x,0)=T_{0}+A\cos(x/R)\quad\mbox{and}\quad T_{2}(x,0)=T_{0}+A\sin(x/R)$
(21)

we need to choose $n=1$, $B_{1}=-A$, $A_{2}=0$, which means that $\Delta
T_{i}(x,0)=f_{i}(x)$, therefore both temperature distributions will remain
invariant and will only decay on time.

|
---|---
(a) | (b)
Figure 3: The graph shows the temperature profiles for both rings for the
following values $D=100mm^{2}/s$, $\rho=1000Kg/m^{3}$, $c=1000J/Kg^{\circ}K$,
$k_{i}=1W/m^{\circ}K$, $a=100mm$, $b=5mm$, $d=1mm$ and for the radii: (a)
$R=21mm$ and $v_{EP}=4.2mm/s$ and (b) $R=22mm$ and $v_{EP}=4.4mm/s$.

Using the same experimental values given in Ref. LABEL:Li, i.e.
$D=100mm^{2}/s$, $\rho=1000Kg/m^{3}$, $c=1000J/Kg^{\circ}K$,
$k_{i}=1W/m^{\circ}K$, $a=100mm$, $b=5mm$ and $d=1mm$, we find that Rabi
oscillations take place for the fundamental wave if the inner ring radius is
between $20mm<R<22mm$ and the rings are rotating with equal but opposite
velocities given by $v_{EP}=h_{c}R$. In Fig. (3) we show the temperature
fields for the unbroken-$\mathcal{PT}$ regime where damped Rabi oscillations
occur in which the maximum and minimum temperatures are 90∘ out of phase.
Let us now consider the case when the rings are rotating with different
velocities close to the exceptional point, specifically we will like to solve
the following system

$\displaystyle\frac{\partial T_{1}}{\partial
t}=D\frac{\partial^{2}T_{1}}{\partial x^{2}}-\left(v_{EP}+\delta
v\right)\frac{\partial T_{1}}{\partial x}+h_{c}(T_{2}-T_{1})$
$\displaystyle\frac{\partial T_{2}}{\partial
t}=D\frac{\partial^{2}T_{2}}{\partial x^{2}}+\left(v_{EP}-\delta
v\right)\frac{\partial T_{2}}{\partial x}+h_{c}(T_{1}-T_{2})$ (22)

where $\delta v<<v_{EP}$. At first it seems that the system given in Eq. (22)
is not APT symmetric, however if we make the following transformation
$\xi=x-\delta vt$ we obtain a system of equations identical to the one given
in Eq. (1) replacing $x\rightarrow\xi$ and $v\rightarrow v_{EP}$. The solution
for Eq. (22) is given by

$\displaystyle\Delta T_{1}=e^{-\lambda h_{c}t}\cos((x-\delta vt)/R)$ (23)
$\displaystyle\Delta T_{2}=-e^{-\lambda h_{c}t}\sin((x-\delta vt)/R)$ (24)

which means that the temperature profiles are moving. This result shows that
we can have a rest-to-motion temperature profile without having equal opposite
rotating velocities.
In summary, we have predicted the existence of Rabi oscillations at the
exceptional point in the diffusive system proposed by Li et al.. We showed
that at the exceptional point the system exhibits two $\mathcal{PT}$ phase
transitions which takes place at critical values for the radii of the rotating
rings. Specifically, if the rings are rotating in opposite directions with
equal tangential velocity given by $v_{EP}=h_{c}/|\kappa|$ and the ring radius
lies between $\sqrt{4D/5h_{c}}<R<\sqrt{D/h_{c}}$ for the fundamental wave,
i.e. $\kappa=\pm 1/R$, the temperature fields exhibit damped Rabi
oscillations. We have shown also that it is not essential to have identical
opposite rotating velocities for the rings in order to have a rest-to-motion
temperature transition, we can do this also by increasing/decreasing the
upper/lower ring velocity away from the exceptional point in order to obtain
traveling wave solutions with positive/negative velocity. Our work reveals the
rich structure of exceptional points in anti-parity-time symmetric diffusive
systems.

I would like to acknowledge support by the program Cátedras Conacyt through
project 1757 and from project A1-S-43579 of SEP-CONACYT Ciencia Básica and
Laboratorio Nacional de Ciencia y Tecnología de Terahertz.

## References

  * (1) C. M. Bender and S. Boettcher, Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry, Phys. Rev. Lett., 80, 5243-5246 (1998).
  * (2) C.M. Bender, S. Böttcher, S. and P. N. Meisinger, PT-symmetric quantum mechanics, J. Math. Phys. 40, 2201–2229 (1999).
  * (3) C. M. Bender, Making sense of non-Hermitian Hamiltonians, Rep. Prog. Phys., 70, 947 (2007).
  * (4) C. M. Bender, D. C. Brody and H. F. Jones, Complex Extension of Quantum Mechanics, Phys. Rev. Lett. 89, 270401 (2002).
  * (5) A. Mostafazadeh, Pseudo-Hermiticity versus PT symmetry: The necessary condition for the reality of the spectrum of a non-Hermitian Hamiltonian, J. of Math. Phys., 43, 205 (2002).
  * (6) C. M. Bender, D. C. Brody and H. F. Jones, Must a Hamiltonian be Hermitian?, Am. J. Phys., 71 (11), 1095-1102 (2003).
  * (7) C. M. Bender, B. K. Berntson, D. Parker and E. Samuel, Observation of PT phase transition in a simple mechanical system, Am. J. Phys., 81 (3), 173-179 (2013).
  * (8) A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. VolatierRavat, V. Aimez, G. A. Siviloglou and D. N. Christodoulides, Observation of PT-Symmetry Breaking in Complex Optical Potentials, Phys. Rev. Lett., 103, No. 093902 (2009).
  * (9) C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev and D. Kip, Observation of Parity-Time Symmetry in Optics, Nat. Phys., 6, 192-195 (2010).
  * (10) L. Chang, X. Jiang,S. Hua, C. Yang, J. Wen, L. Jiang, G. Li, G. Wang and M. Xiao, Parity-Time Symmetry and Variable Optical Isolation in Active-Passive-Coupled Microresonators, Nat. Photonics, 8, 524-529 (2014).
  * (11) M. Wimmer, A. Regensburger, M. A. Miri, C. Bersch, D. N. Christodoulides and U. Peschel, Observation of Optical Solitons in PT-Symmetric Lattices, Nat. Commun., 6, 7782 (2015).
  * (12) P. Peng, W. Cao, C. Shen, W. Qu, J. Wen, L. Jiang and Y. Xiao, Anti-Parity-Time Symmetry with Flying Atoms, Nat. Phys. 12 1139-1145 (2016).
  * (13) Y. Jiang, Y. Mei, Y. Zuo, Y. Zhai, J. Li, J. Wen and S. Du, Anti-Parity-Time Symmetric Optical Four-Wave Mixing in Cold Atoms, Phys. Rev. Lett. 123, 193604 (2019).
  * (14) Y. Choi, C. Hahn, J. W. Yoon and S.H. Song, Observation of an Anti-PT-Symmetric Exceptional Point and Energy-Difference Conserving Dynamics in Electrical Circuit Resonators, Nat. Commun. 9, 2182 (2018).
  * (15) X. L. Zhang, T. Jiang and C.T. Chan, Dynamically Encircling an Exceptional Point in Anti-Parity-Time Symmetric Systems: Asymmetric Mode Switching for Symmetry-Broken Modes, Light: Sci. Appl. 8, 88 (2019)
  * (16) C. F. Li and G. C. Guo, Experimental Simulation of Anti-Parity-Time Symmetric Lorentz Dynamics, Optica 6, 67 (2019).
  * (17) J. Zhao, Y. Liu, L. Wu, C. K. Duan, Y. X. Liu, and J. Du, Observation of Anti-PT-Symmetry Phase Transition in the Magnon-Cavity-Magnon Coupled System, Phys. Rev. Appl. 13, No. 014053 (2020).
  * (18) Ying Li, Yu-Gui Peng, Lei Han, Mohammad-Ali Miri, Wei Li, Meng Xiao, Xue-Feng Zhu, Jianlin Zhao, Andrea Alu, Shanhui Fan and Cheng-Wei Qiu, Anti–parity-time symmetry in diffusive systems, Science 364, 6436, 170-173 (2019)
  * (19) J. Doppler, A. A. Mailybaev, J. Böhm, A. G. U. Kuhl and, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Dynamically encircling an exceptional point for asymmetric mode switching, Nature 537, 76–79 (2016).
  * (20) H. Xu, D. Mason, L. Jiang, and J. G. E. Harris, Topological energy transfer in an optomechanical system with exceptional points, Nature 537, 80–83 (2016).
  * (21) Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, and D. N. Christodoulides, Unidirectional invisibility induced by PT-symmetric periodic structures, Phys. Rev. Lett. 106, 213901 (2011).
  * (22) B. Peng, Ş. K. Özdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, Parity-time-symmetric whispering-gallery microcavities, Nat. Phys. 10 (2014).
  * (23) L. Feng, Y.-L. Xu, W. S. Fegadolli, M.-H. Lu, J. E. Oliveira, V. R. Almeida, Y.-F. Chen, and A. Scherer, Experimental demonstration of a unidirectional reflectionless parity-time metamaterial at optical frequencies, Nat. Mater. 12, 108–113 (2013).
  * (24) L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, and X. Zhang, Single-mode laser by parity-time symmetry breaking, Science 346, 972–975 (2014).
  * (25) H. Hodaei, M.-A. Miri, M. Heinrich, D. N. Christodoulides, and M. Khajavikhan, Parity-time–symmetric microring lasers, Science 346, 975–978 (2014).
  * (26) Y. Sun, W. Tan, H.Q. Li, J. Li and H. Chen, Experimental demonstration of a coherent perfect absorber with PT phase transition, Phys. Rev. Lett. 121 073901 (2018)