y)+0.4\cos(6\pi x)\cos(6\pi y)$ $\displaystyle\quad+0.3\cos(10\pi x)\cos(4\pi
y)+0.02\sin(2\pi y)+0.02\sin(2\pi x)$ (6.46) $\displaystyle b(0,x,y)$
$\displaystyle=\sin(2\pi y)-1,$

see Figure 4 for illustrations, as well as the following bathymetry and
rotation fields

(6.47) $\displaystyle h(x,y)$ $\displaystyle=\cos(2\pi x)+0.5\cos(4\pi
x)+\frac{1}{2}\cos(6\pi x)$ (6.48) $\displaystyle f(x,y)$
$\displaystyle=0.4\cos(4\pi x)\cos(4\pi y).$

For time stepping, we used $\Delta t=0.0005$ in all cases to facilitate
comparisons. This choice satisfies the CFL condition for the $\alpha=0$ case.

We computed each solution for $5000$ time steps111Go to
https://youtu.be/a2a4xzft3Pg for a video of the full simulation.. Figures 5, 6
and 7 show $\alpha$-TQG solution snapshots of buoyancy, potential vorticity
and velocity magnitude respectively, at the $1280$’th and $2600$’th time
steps, for $\alpha$ values $0,\frac{1}{128^{2}},\frac{1}{64^{2}}$ and
$\frac{1}{16^{2}}$. The interpretation of the regularisation parameter
$\alpha$ is that its square root value corresponds to the fraction of the
domain’s length scale that get regularised. At the $1280$’th time step, the
flows are in early spin-up phase. At the $2600$’th time step, although the
flows are still in spin-up phase, much more flow features have developed. The
sub-figures illustrate via comparisons, how $\alpha$ controls the development
of small scale features and instabilities. Increasing $\alpha$ leads to more
regularisation at larger scales.

In view of Proposition 4.1, we investigated numerically the convergence of
$\alpha$-TQG to TQG using our numerical setup. Consider the relative error
between TQG buoyancy and $\alpha$-TQG buoyancy,

(6.49) $\displaystyle e_{b}(t,\alpha)$
$\displaystyle:=\frac{\|b(t)-b^{\alpha}(t)\|_{H^{1}}}{\|b(t)\|_{H^{1}}}$

and the relative error between TQG potential vorticity and $\alpha$-TQG
potential vorticity,

(6.50) $\displaystyle e_{\omega}(t,\alpha)$
$\displaystyle:=\frac{\|\omega(t)-\omega^{\alpha}(t)\|_{L^{2}}}{\|\omega(t)\|_{L^{2}}}.$

The norms were chosen in view of Proposition 4.1. Figures 8(a) and 8(b) show
plots of $e_{b}(t,\alpha)$ and $e_{\omega}(t,\alpha)$ as functions of time
only, for fixed $\alpha$ values
$\frac{1}{16^{2}},\frac{1}{32^{2}},\frac{1}{64^{2}},\frac{1}{128^{2}}$ and
$\frac{1}{220^{2}}$, from the initial time up to and including the $2800$’th
time step. We observe that starting from $0$, the relative errors $e_{b}$ and
$e_{\omega}$ initially increase over time but plateau at around and beyond the
$2000$’th time step. Up to the $1400$’th time step, the plotted relative
errors remain less than $1.0$, and are arranged in the ascending order of
$\alpha$, i.e. smaller $\alpha$’s give smaller relative errors.

We note that, if we compare and contrast the relative error results with the
solution snapshots of buoyancy and potential vorticity, particularly at the
$2600$’th time step (shown in Figures 5(b) and 6(b)), we observe
”discrepancies” – for lack of a better term – between the results. For
example, in Figure 5(b), if we compare the $\alpha=\frac{1}{128^{2}}$ snapshot
with the $\alpha=\frac{1}{16^{2}}$ one, the former show small scale features
that are much closer to those that exist in the reference TQG $\alpha=0$
solution. However, according to the relative error measurements, the two
regularised solutions are more or less equivalent in their differences to the
reference TQG solution at time step $2600$.

If we fix the value of the time parameter in $e_{b}(t,\alpha)$ and
$e_{\omega}(t,\alpha)$, and vary $\alpha$, we can estimate the convergence
rate of $\alpha$-TQG to TQG for our numerical setup. Proposition 4.1 predicts
a convergence rate of $1$ in $\alpha$, using the $H^{1}$ norm on the buoyancy
differences and the $L^{2}$ norm on the potential vorticity differences. Note
that, in (4.1), the left hand side evaluates a supremum over a given compact
time interval. We see in Figures 8(a) and 8(b) that the relative errors are
monotonic up to around the $1500$’th time step mark. So for a given time point
$T$ between the initial time and the $1500$’th time step, we assume we can
evaluate $e_{b}$ and $e_{\omega}$ at $T$ to estimate the supremum over
$[0,T]$.

Figures 9(a) and 9(b) show plots – in log-log scale – of $e_{b}$ and
$e_{\omega}$ as functions of $\alpha$ only at the $600$’th, $800$’th,
$1000$’th and $1200$’th time steps, for all the chosen $\alpha$ values. These
numbers of time steps correspond to $T=0.3$, $T=0.4$, $T=0.5$ and $T=0.6$
respectively. We also plotted in Figures 9(a) and 9(b) linear functions of
$\alpha$ to provide reference order 1 slopes. Comparing the results to the
reference, we see that the theoretical convergence rate of 1 is attained for
$T=0.3$ ($600$’th time step), $T=0.4$ ($800$’th time step) and $T=0.5$
($1000$’th time step). However for $T=0.6$ ($1200$’th time step), the
convergence rates are at best order $1/2$.

Figure 4. Initial conditions – buoyancy (6.46) on the left, and potential
vorticity (6.45) on the right. The different colours of the PV plot correspond
to positive and negative values of PV, which can be interpreted as clockwise
and anticlockwise eddies respectively.

(a)

(b)

Figure 5. Comparisons of solution snapshots of the $\alpha$-TQG buoyancy field
that correspond to four different values of $\alpha$, at two different points
in time. Potential vorticity and velocity magnitude snapshots of the same
solutions are shown in Figure 6 and Figure 7 respectively. Shown in each sub-
figure are the results corresponding to $\alpha=0$ (top left),
$\alpha=\frac{1}{128^{2}}$ (top right), $\alpha=\frac{1}{64^{2}}$ (bottom
right) and $\alpha=\frac{1}{16^{2}}$ (bottom left). When $\alpha=0$, the
solution is of the TQG system. Sub-figure (A) shows the solutions at the
$1280$’th time step, or equivalently when $t=0.64$. Sub-figure (B) shows the
solutions at the $2600$’th time step, or equivalently when $t=1.30$. The flows
in both sub-figures are at different stages of spin-up, with the flow at
$t=1.30$ showing more developed features. The $\alpha$ parameter is
interpreted as the fraction of the domain’s length-scale squared value at
which regularisation is applied. Due to regularisation, the flows do develop
differently. Nevertheless, we observe that as $\alpha$ gets smaller, the
$\alpha$-TQG flow features converge to that of the TQG flow.

(a)

(b)

Figure 6. Comparisons of solution snapshots of the $\alpha$-TQG potential
vorticity field that correspond to four different values of $\alpha$, at two
different points in time. See the caption of Figure 5 for explanations of the
arrangements of these plots. Buoyancy and velocity magnitude snapshots of the
same solutions are shown in Figure 5 and Figure 7 respectively. As in Figure
4, the different colours of the PV field correspond to positive and negative
values of PV, which can be interpreted as clockwise and anticlockwise eddies
respectively. In sub-figure (B), we observe more clearly the regularisation
effects of $\alpha$. As $\alpha$ increases, the flows develop in different
ways – features at smaller scales no longer develop, and larger features
evolve differently as a result. Nevertheless, we observe that as $\alpha$ gets
smaller, the $\alpha$-TQG flow features converge to that of the TQG flow.

(a)

(b)

Figure 7. Comparisons of solution snapshots of the $\alpha$-TQG velocity
magnitudes that correspond to four different values of $\alpha$, at two
different points in time. See the caption of Figure 5 for explanations of the
arrangements of these plots. Buoyancy and potential vorticity snapshots of the
same solutions are shown in Figure 5 and Figure 6 respectively. Using the same
scale for colouring, we observe in both sub-figures the strength of the
colours weaken as $\alpha$ increases. This indicates smoothing of large
velocity magnitudes. Additionally, in sub-figure (B), we observe how features
at smaller scales get smoothed out. In particular, in the
$\alpha=\frac{1}{16^{2}}$ plot, we see that essentially only large scale
vorticies remain. Hence, considering the velocity field as a part of the
system’s nonlinear advection operator, these figures help us to better
visualise how the flow features developed in Figure 5 and Figure 6.

(a)

(b)

Figure 8. Sub-figures (A) and (B) show, respectively, plots of the relative
error functions $e_{b}(t,\alpha)$ and $e_{\omega}(t,\alpha)$ as functions of
time only, for the fixed $\alpha$ values $\frac{1}{16^{2}}$,
$\frac{1}{32^{2}}$, $\frac{1}{64^{2}}$, $\frac{1}{128^{2}}$ and
$\frac{1}{220^{2}}$. See equations (6.49) and (6.50) for the definitions of
$e_{b}$ and $e_{\omega}$ respectively. The plots are shown from $t=0$ up to
and including $t=1.4$, which is equivalent to $2800$ time steps using $\Delta
t=0.0005$. We observe that, starting from $0$, all plots of $e_{b}$ and
$e_{\omega}$ increase initially, and plateau at around the $2000$’th time
step. Further, up to the $1400$’th time step the plotted relative errors
remain less than $1.0$ and are arranged in the ascending order of $\alpha$.
Relating these relative error results at the $2600$’th time step to the
solution snapshots at the same time point (shown in figures 5(b) and 6(b)), we
see that although the snapshots show convergence of flow features, the
relative error values suggest all the $\alpha$-TQG solutions are more or less
equally far away from the TQG flow.

(a)

(b)

Figure 9.  Sub-figures (A) and (B) show, respectively, plots in log-log scale
of the relative error functions $e_{b}(t,\alpha)$, and $e_{\omega}(t,\alpha)$
as functions of $\alpha$ only, at the fixed time values $t=0.3$ (equivalently
at the $600$’th time step), $t=0.4$ (equivalently at the $800$’th time step),
$t=0.5$ (equivalently at the $1000$’th time step) and $t=0.6$ (equivalently at
the $1200$’th time step). In view of Proposition 4.1, we assumed that, based
on the initial monotonicity of the relative errors (see Figure 8), the
supremum over $[0,T]$ in (4.1) can be estimated by evaluating $e_{b}$ and
$e_{\omega}$ at $T$. Proposition 4.1 predicts that up to a certain time, the
rate of convergence of $e_{b}$ and $e_{\omega}$ to $0$ should be no less than
order 1 in $\alpha$. In both sub-figures we have plotted, as the reference for
comparison, linear functions of $\alpha$. Thus, comparing slopes to the
reference, we see that numerically we get order $1$ convergence for $t=0.3$,
$t=0.4$ and $t=0.5$. However, for $t=0.6$, the rate of convergence is no more
than $1/2$.

##  7\. Conclusion and outlook

In conclusion, we have formulated the thermal quasi-geostrophic (TQG) model
equations in the regime of approximations relevant to GFD and we have explored
their analytical and numerical properties. In respect to the analytical
properties, we have shown that the TQG model and its $\alpha$-regularized
version, the $\alpha$-TQG model, admit unique local strong solutions that are
stable in a larger space with weaker norm, that both models have a maximum
time of existence, solutions of the latter model converge to solutions of the
former model on any time interval in which a solution to the former lives, and
we have also determined conditions under which the solution of the
$\alpha$-TQG will blow up.

With respect to numerics, we described our discretisation methods for
approximating TQG and $\alpha$-TQG solutions and showed that the FEM semi-
discrete scheme conserves numerical energy in the TQG case. Using example
simulation results, we numerically verified that $\alpha$-TQG solutions
converge to that of the TQG system and attains the order $1$ in $\alpha$
convergence rate predicted by Proposition 4.1. Additionally, we derived a
dispersion relation for the linearised $\alpha$-TQG system and showed how
$\alpha$-regularisation could be used to control the development of high
wavenumber instabilites. Given that we have convergence of $\alpha$-TQG
solutions to TQG, the linear stability results can be viewed as
generalisations of those shown in [17] for the TQG system.

We now end this section with some open problems.

  * •

Can one construct a global-in-time strong solution (or a global weak solution)
of either the TQG model equations (1.1)–(1.5), or the $\alpha$-TQG model model
equations (2.12)–(2.14)?

  * •

Can one give a Beale–Kato–Majda condition for the blowup of a strong solution
to the TQG in terms of a single unknown? That is, is there a BKM condition in
terms of either $b$ or $\omega$ (or $\mathbf{u}$) that does not require a
combination of both variables?

  * •

As discussed in Section 1.3, an application of the Stochastic Advection by Lie
Transport (SALT) approach to the Hamiltonian formulation of the TQG equations
results in a stochastic Hamiltonian formulation of TQG equations in
(1.6)-(1.8). This stochastic version of TQG preserves an infinite family of
integral conserved quantities, as is shown in [17]. The SALT version of TQG
represents an outstanding challenge for uncertainty quantification and data
assimilation which will surely spur us on to further investigations of these
problems.

  * •

As also discussed in Section 1.3, the parallels between TQG and the Rayleigh-
Bénard equations should also attract our attention for the consideration of a
SALT version of the deterministic Rayleigh-Bénard convection equations.

## Acknowledgements

This work has been partially supported by European Research Council (ERC)
Synergy grant STUOD-DLV-856408. We would like to thank our friends and
colleagues for their encouraging comments, especially W. Bauer, F. J. Beron-
Vera, B. Chapron, L. Cope, C. J. Cotter, E. Dinvay, O. Lang, E. Memin, A.
Radomska - Botelho Moniz, S. Takao.

##  8\. Appendix

We present a result in this section, which along with its proof, can be found
in [20]. Our construction of local solutions relies on this result and thus,
we present it here for the sake of completeness.
We consider the abstract equation

(8.1) $\displaystyle\frac{\partial}{\partial t}u+\mathcal{A}(t,u)=0,\qquad
t\geq 0,\qquad u\big{|}_{t=0}=u_{0},$

where $\mathcal{A}$ is a nonlinear operator.

######  Theorem 8.1.

Let $\\{V,H,X\\}$ be real separable Banach spaces such that

  * •

the embeddings $V\hookrightarrow H\hookrightarrow X$ are continuous and dense;

  * •

$H$ is a Hilbert space;

  * •

There is a continuous, nondegenerate bilinear form $(\,,\,)$ on ${X\times V}$
such that $(u,v)=\langle u,v\rangle_{H}$ for $u\in H$ and $v\in V$.

Let $\mathcal{A}$ be a (sequentially) weakly continuous map on
$[0,T_{*}]\times H$ into $X$ such that

(8.2)
$\displaystyle(\mathcal{A}(t,v),v)\geq-\rho(\|v\|_{H}^{2})\qquad\text{for}\quad
t\in[0,T_{*}],\quad v\in V,$

where $\rho(r)\geq 0$ is a monotone increasing function of $r\geq 0$. Then for
any $u_{0}\in H$, there exists $T>0$, $T\leq T_{*}$, and a solution $u$ of
(8.1) in the class

(8.3) $\displaystyle u\in C_{w}([0,T];H)\cap C_{w}^{1}([0,T];X).$

Moreover,

$\displaystyle\|u(t)\|_{H}^{2}\leq r(t),\qquad t\in[0,T],$

where $r$ is a monotone increasing function on $[0,T]$; $T$ and $r$ can be
chosen so as to depend only on $\rho$ and $\|u_{0}\|_{H}$.

######  Remark 8.2.

(a) If $\mathcal{A}$ maps $[0,T_{*}]\times V$ into $H$, we may replace (8.2)
with

(8.4)
$\displaystyle\langle\mathcal{A}(t,v),v\rangle_{H}\geq-\rho(\|v\|_{H}^{2})\qquad\text{for}\quad
t\in[0,T_{*}],\quad v\in V.$

(b) $T$ and $r$ can be determined by solving the scalar differential equation

(8.5) $\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}r=2\rho(r),\qquad
r(0)=\|u_{0}\|_{H}^{2}.$

$T$ may be any value such that $r$ exists on $[0,T]$. If the solution to (8.5)
is not unique, $r$ should be the maximal soluton.
(c) $u(t)\rightarrow u_{0}$ holds strongly in $H$ as $t\rightarrow 0$. In
other words, the solution (8.3) is strongly continuous at the initial time
$t=0$.

## References

  * [1] Anderson, D.L., McCreary Jr, J.P.: On the role of the indian ocean in a coupled ocean–atmosphere model of el niño and the southern oscillation.  Journal of the atmospheric sciences 42(22), 2439–2442 (1985)
  * [2] Beale, J.T., Kato, T., Majda, A.: Remarks on the breakdown of smooth solutions for the $3$-D Euler equations.  Comm. Math. Phys. 94(1), 61–66 (1984)
  * [3] Beier, E.: A numerical investigation of the annual variability in the gulf of california.  Journal of physical oceanography 27(5), 615–632 (1997)
  * [4] Bernsen, E., Bokhove, O., van der Vegt, J.J.: A (Dis)continuous finite element model for generalized 2D vorticity dynamics.  Journal of Computational Physics 211(2), 719–747 (2006).  DOI 10.1016/j.jcp.2005.06.008.  URL https://doi.org/10.1016/j.jcp.2005.06.008
  * [5] Beron-Vera, F.: Multilayer shallow-water model with stratification and shear.  Revista Mexicana de Física 67, 351–364 (2021).
  * [6] Beron-Vera, F.: Nonlinear saturation of thermal instabilities.  Physics of Fluids 33(3), 036608 (2021).
  * [7] Brenner, S.C., Scott, L.R.: The Mathematical Theory of Finite Element Methods.  Springer New York (2008).  DOI 10.1007/978-0-387-75934-0.  URL https://doi.org/10.1007/978-0-387-75934-0
  * [8] Brézis, H., Gallouet, T.: Nonlinear Schrödinger evolution equations.  Nonlinear Anal. 4(4), 677–681 (1980)
  * [9] Cao, Y., Jolly, M.S., Titi, E.S., Whitehead, J.P.: Algebraic bounds on the rayleigh–bénard attractor.  Nonlinearity 34(1), 509 (2021)
  * [10] Dinvay, E.: Well-posedness for a Whitham-Boussinesq system with surface tension.  Math. Phys. Anal. Geom. 23(2), Paper No. 23, 27 (2020)
  * [11] Foias, C., Holm, D.D., Titi, E.S.: The Navier-Stokes-alpha model of fluid turbulence.  pp. 505–519 (2001).  Advances in nonlinear mathematics and science
  * [12] Foias, C., Holm, D.D., Titi, E.S.: The three dimensional viscous Camassa-Holm equations, and their relation to the Navier-Stokes equations and turbulence theory.  J. Dynam. Differential Equations 14(1), 1–35 (2002)
  * [13] Gibson, T.H., McRae, A.T., Cotter, C.J., Mitchell, L., Ham, D.A.: Compatible Finite Element Methods for Geophysical Flows.  Springer International Publishing (2019).  DOI 10.1007/978-3-030-23957-2.  URL https://doi.org/10.1007/978-3-030-23957-2
  * [14] Hesthaven, J.S., Warburton, T.: Nodal Discontinuous Galerkin Methods.  Springer New York (2008).  DOI 10.1007/978-0-387-72067-8.  URL https://doi.org/10.1007/978-0-387-72067-8
  * [15] Holm, D.D.: Variational principles for stochastic fluid dynamics.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471(2176), 20140,963 (2015)
  * [16] Holm, D.D., Luesink, E.: Stochastic wave-current interaction in thermal shallow water dynamics.  Journal of Nonlinear Science 31(2), 1–56 (2021).  DOI 10.1007/s00332-021-09682-9.  URL https://doi.org/10.1007/s00332-021-09682-9
  * [17] Holm, D.D., Luesink, E., Pan, W.: Stochastic mesoscale circulation dynamics in the thermal ocean.  Physics of Fluids 33(4), 046,603 (2021).  DOI 10.1063/5.0040026.  URL https://doi.org/10.1063/5.0040026
  * [18] Holm, D.D., Marsden, J.E., Ratiu, T.S.: The euler–poincaré equations and semidirect products with applications to continuum theories.  Advances in Mathematics 137(1), 1–81 (1998)
  * [19] Kato, T.: Liapunov functions and monotonicity in the navier-stokes equation.  In: Functional-analytic methods for partial differential equations, pp. 53–63. Springer (1990)
  * [20] Kato, T., Lai, C.Y.: Nonlinear evolution equations and the Euler flow.  J. Funct. Anal. 56(1), 15–28 (1984)
  * [21] Kato, T., Ponce, G.: Commutator estimates and the Euler and Navier-Stokes equations.  Comm. Pure Appl. Math. 41(7), 891–907 (1988)
  * [22] Klainerman, S., Majda, A.: Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids.  Commun. Pure Appl. Math. 34(4), 481–524 (1981)
  * [23] Marsden, J.E., Shkoller, S.: The anisotropic Lagrangian averaged Euler and Navier-Stokes equations.  Arch. Ration. Mech. Anal. 166(1), 27–46 (2003)
  * [24] McCreary Jr, J.P., Zhang, S., Shetye, S.R.: Coastal circulations driven by river outflow in a variable-density 1$1/2$-layer model.  Journal of Geophysical Research: Oceans 102(C7), 15,535–15,554 (1997)
  * [25] McWilliams, J.C.: A survey of submesoscale currents.  Geoscience Letters 6(1), 1–15 (2019)
  * [26] O’Brien, J.J., Reid, R.O.: The non-linear response of a two-layer, baroclinic ocean to a stationary, axially-symmetric hurricane: Part i. upwelling induced by momentum transfer.  Journal of Atmospheric Sciences 24(2), 197–207 (1967)
  * [27] Ripa, P.: Conservation laws for primitive equations models with inhomogeneous layers.  Geophysical & Astrophysical Fluid Dynamics 70(1-4), 85–111 (1993)
  * [28] Ripa, P.: Low frequency approximation of a vertically averaged ocean model with thermodynamics.  Revista Mexicana de Física 42(1), 117–135 (1996).
  * [29] Ripa, P.: On improving a one-layer ocean model with thermodynamics.  Journal of Fluid Mechanics 303, 169–201 (1995)
  * [30] Ripa, P.: On the validity of layered models of ocean dynamics and thermodynamics with reduced vertical resolution.  Dynamics of atmospheres and oceans 29(1), 1–40 (1999)
  * [31] Schopf, P.S., Cane, M.A.: On equatorial dynamics, mixed layer physics and sea surface temperature.  Journal of physical oceanography 13(6), 917–935 (1983)
  * [32] Volkov, D.L., Kubryakov, A.A., Lumpkin, R.: Formation and variability of the lofoten basin vortex in a high-resolution ocean model.  Deep Sea Research Part I: Oceanographic Research Papers 105, 142–157 (2015)
  * [33] Wang, C., Zhang, Z.: Global well-posedness for the 2-D Boussinesq system with the temperature-dependent viscosity and thermal diffusivity.  Adv. Math. 228(1), 43–62 (2011)


††thanks: Work done while at Indian Institute of Science Bangalore††thanks:
Corresponding Author and Work done while at Shell Technology Centre
Bangalore.††thanks: Work done while at Indian Institute of Science
Bangalore††thanks: Corresponding Author and Work done while at Shell
Technology Centre Bangalore.

# Potential energy surface prediction of Alumina polymorphs using graph neural
network

Soumya Sanyal  University of Southern California, USA  Arun Kumar Sagotra
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>Narendra Kumar
Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore  Sharad
Rathi  Vellore Institute of Technology, Vellore  Mohana Krishna  Vellore
Institute of Technology, Vellore  Nagesh Somayajula  Centura, Colarado, USA
Duraivelan Palanisamy  Shell Technology Centre Bangalore, Bangalore  Ram R.
Ratnakar  Shell Technology Centre Houston, USA  Suchismita Sanyal  Shell
Technology Centre Bangalore, Bangalore  Partha Talukdar  Google Research,
Bangalore, India  Umesh Waghmare  Jawaharlal Nehru Centre for Advanced
Scientific Research, Bangalore  Janakiraman Balachandran
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
The process of design and discovery of new materials can be significantly
expedited and simplified if we can learn effectively from available data. Deep
learning (DL) approaches have recently received a lot of interest for their
ability to speed up the design of novel materials by predicting material
properties with precision close to experiments and ab-initio calculations. The
application of deep learning to predict materials properties measured by
experiments are valuable yet challenging due to the limited amount of
experimental data. Most of the existing approaches to predict properties from
computational data have also been directed towards specific material
properties. In this work, we extend this approach, by proposing Landscape
Crystal Graph Convolution Network(LCGCN), an accurate and transferable deep
learning framework based on graph convolutional networks. LCGCN directly
learns the potential energy surface (PES) from atomic configurations. This
approach can enable transferable models that can predict different material
properties. We apply this framework to bulk crystals (i.e. Al2O3), and test it
by calculating potential energy surfaces at different temperatures and across
different phases of crystal.

##  I Introduction

The development of next generation technologies in energy, chemicals,
electronics and transportation industries are closely tied to the development
of new materials with the desired properties [1]. For example in case of
energy, the need of the hour is the development of inexpensive and
environmentally safe materials that can generate and store electrons with very
high density and efficiency while operating at high power. In case of
chemicals, there is a critical need to develop new materials with better
catalytic performance [2], as well as those which can make chemicals from
alternative feedstock such as recycled plastics.

Historically, discovery of materials has been driven by empiricism, through
the knowledge and intuition of researchers working on these problems. Over the
past few decades, they are being augmented by the insights and knowledge
derived from physics based models, in particular, atomistic models and
simulations that can predict the various properties of interest for a wide
range of materials [3, 4]. However, these physics based models are
computationally expensive and as a result, only aid us to explore a small
subset of the materials search space. Further, within this approach, there is
no formal way to transfer information and knowledge from one system to
another, and from one application to another. Data driven models based on
statistics and machine learning can help in speeding up the materials
discovery process by providing a framework where algorithms can learn from
existing data obtained from models and experiments on a subset of materials
and predict the properties and functionalities of new materials. In order to
facilitate such developments, various research groups from across the world
have also provided free public access to the datasets of material structures
[5, 6, 7, 8, 9] and properties obtained predominantly from atomistic
simulations.

These datasets and other datasets obtained from experiments have been used in
the recent literature to develop data driven models to predict material
properties [10, 11, 12, 13]. Most of these models are based on hand crafted
descriptors that are suitable only for certain subset of materials and for
targeted properties of interest. However, a transferable model needs to
operate on raw data to remove these biases introduced through hand crafted
descriptors. In case of materials and molecules, such raw input data are the
atoms and their position coordinates. This basic materials information can be
naturally expressed in the form of graphs. Graph Neural Networks (GNNs) are a
class of deep learning models that have been successfully used to model graph-
structured data [14, 15, 16]. GNNs can learn to effectively represent a
molecular/crystal graph from raw data of atom/bond features and use neural
networks to perform regression and classification tasks [17, 18, 19, 20].
Recently, Xie $et\,\,al.$ [13] developed crystal graph convolution neural
network (CGCNN) that was able to take this basic structural and chemical
information of a crystal as graph to perform graph convolutions on these
graphs. They were able to obtain excellent results to predict properties of a
wide range of crystals, in their ground-state equilibrium structure. Later,
Sanyal $et\,\,al.$ [21] integrated CGCNN with multi-task learning to further
improve the prediction accuracy.

##  II Outline of the Work

Although CGCNN has been used to predict equilibrium properties of materials,
many thermodynamic, kinetic and transport properties at non-zero temperatures
are determined only after incorporating the phononic contributions that
distort the crystal and its atomic components from the equilibrium position.
These distortion and the resultant phononic contributions depend upon the
potential energy surface (PES) associated with it. PES is a highly non-linear
complex function of the chemical composition, symmetry and atomic arrangement
of the crystals. Even a small change in these parameters can significantly
change the PES. The PES is the input to various statistical approaches such as
Monte Carlo and molecular dynamics simulations to predict the phononic degrees
of freedom and their influence on various material properties.

Conventionally PES are fitted to explicit function forms from ab-initio
simulation results. However, over the past few years gaussian process based
machine learning has been applied to develop interatomic potentials for
various materials ranging from elemental (amorphous carbon, graphene, boron)
[22, 23, 24] to multi-elemental (methane) [25]. Although very successful,
development of Gaussian Approximation Potential (GAP) [26] requires careful
design of experiments to simulate small systems whose inputs are needed to fit
the interatomic potentials. Further, there are also challenges in transfering
the potential to a new system.

These shortcomings can potentially be overcome by GNN based deep learning
approach, in which a model learns from the raw materials data. These models
can possibily lead to better transferability across different materials
composition, crystal symmetry, atomic arrangements and environments. In this
work, we provide an initial demonstration of such transferability and propose
Landscape Crystal Graph Convolution Network (LCGCN), an extension of CGCNN
architecture, to predict the PES of alumina $\left(Al_{2}O_{3}\right)$, a
material that serves as the base material for various adsorbent and catalytic
applications [27, 28, 29, 30]. Specifically, we demonstrate the
generalizability of this approach by developing LCGCN models that not only
accurately predict PES, but are also transferable across scales, environmental
conditions (temperature) and different crystallographic polymorphs, trigonal
phase (space group R$\overline{3}c$), cubic phase (space group Ia3),
Orthorhombic (space groups Pbca and Pna21) of alumina.

##  III Methodology

###  III.1 Molecular Dynamics

The data required to train, validate and test the CGCNN models were obtained
by performing Molecular Dynamics (MD) simulations. We have performed MD
simulations under $(N,V,T)$ conditions with the LAMMPS code [31]. In these
simulations the temperature is kept fluctuating around a setpoint value by
using Nosé - Hoover thermostats. We employ simulation boxes, typically varying
from $30$ to $480$ atoms, and apply periodic boundary conditions along the
three Cartesian directions. Newton’s equations of motion are integrated using
the customary Verlet’s algorithm with a time-step length of $10^{-3}$ ps. The
typical duration of a MD run is of $10$ ps. We used interatomic potentials of
the Second-Moment tight-Binding-QEq (SMTB-Q) form [32].

In this work, employing the above MD protocol, we simulate $Al_{2}O_{3}$
systems comprising of different sizes, environmental temperatures and
polymorphs. We pick 10,000 configurations of these systems after thermal
equilibriation. Unless specified otherwise, the 10,000 configurations were
divided into a ratio of 6:2:2 for training, validation and testing of the
LCGCN model.

Figure 1:  Schematic representation of using LCGCN to predict potential energy
of distorted alumina crystals. Initially, the crystal structure is represented
as a graph with nodes and edges representing the atoms and bonds in the
crystal respectively. Here, $\bm{v_{i}},\bm{v_{j}}$ and $\bm{u_{ij}}$ denote
the node, neighbor and edge embeddings respectively. Then, CGCNN is used to
learn a crystal level embedding $\bm{v}_{\mathcal{G}}$ which is used to
predict the potential energy. Please refer to Sec. III.2 for more details.

###  III.2 Landscape Crystal Graph Convolution Network(LCGCN)

CGCNN [13] is a graph neural network which focuses on building a crystal graph
from a given crystal structure and utilizes graph convolution networks (GCNs)
to model the crystal graph and predict their properties with accuracy
comparable to ab initio physics models. We propose Landscape Crystal Graph
Convolution Network (LCGCN), an extension of CGCNN, that can be used to
predict the potential energy surface (PES) of crystal structures. Next, we
discuss the details of proposed method LCGCN in detail.

A crystal graph $\mathcal{G}$ is an undirected multigraph defined by nodes
representing atoms and edges representing bonds in a crystal. Formally, let
$\mathcal{G}\mathord{=}(\mathcal{A},\mathcal{E},\mathcal{V},\mathcal{U})$
denote a crystal graph. Here $\mathcal{A}$ represents the set of atoms in the
crystal structure,
$\mathcal{E}\mathord{=}\\{(i,j)_{k}\mathord{:}~{}k^{th}~{}\text{bond between
atom}~{}i~{}\text{and}~{}j\\}$, is the set of undirected edges denoting the
bonds and $|\mathcal{A}|\mathord{=}N$ is the number of atoms in the crystal
graph. $\bm{v_{i}}\in\mathcal{V}$ contains the features of the $i^{th}$ atom
which encodes various properties of the atom.
$\bm{u}_{(i,j)_{k}}\in\mathcal{U}$ is the feature vector for the $k^{th}$ bond
between atoms $i$ and $j$ that captures properties of the edges. To capture
the influence of local neighborhood in the crystal structure, a graph
convolution formulation similar to CGCNN [13] is used defined as follows,

$\bm{v}_{i}^{(t+1)}=\bm{v}_{i}^{(t)}+\sum_{j,k}\sigma(\bm{z}_{(i,j)_{k}}^{(t)}\bm{W}_{c}^{(t)}+\bm{b}_{c}^{(t)})\odot
g(\bm{z}_{(i,j)_{k}}^{(t)}\bm{W}_{s}^{(t)}+\bm{b}_{s}^{(t)}),$

where
$\bm{z}_{(i,j)_{k}}^{(t)}=\bm{v}_{i}^{(t)}\oplus\bm{v}_{j}^{(t)}\oplus\bm{u}_{(i,j)_{k}}$
denotes the concatenation of atom and bond feature vectors of the neighbors of
$i^{th}$ atom. $\odot$ denotes element-wise multiplication and $\sigma$
denotes a sigmoid function. The $\sigma(\bm{\cdot})$ factor acts as a learned
weight matrix to incorporate different interaction strengths between
neighbors. $\bm{W}_{c}^{(t)}$, $\bm{W}_{s}^{(t)}$, $\bm{b}_{c}^{(t)}$ and
$\bm{b}_{s}^{(t)}$ are the convolution weight matrix, self weight matrix,
convolution bias and self bias of the $t$-th layer of GCN respectively, and
$g(\cdot)$ is the activation function for introducing non-linear coupling
between layers.

Graph pooling [33] is a technique used to learn graph representations given
the learnt node and edge features. The learnt atom features are pooled to get
a vector representation of the crystal $\left(\bm{v}_{\mathcal{G}}\right)$.
The crystal representation $\bm{v}_{\mathcal{G}}$ is then fed to a network of
fully-connected layers with non-linearities that learns to predict a property
of the crystal in a supervised manner. The workflow is depicted in Fig. 1. To
train the network, we calculate the loss using mean squared error between the
true values obtained using MD simulations and the model predictions followed
by back-propagation [34]. The accuracy of the model is quantified through mean
absolute error per atom (MAE/atom) by comparing the model predictions on test
configurations against the true values obtained using MD simulations.

##  IV Results and Discussion

Pooling Operator | Formulation
---|---
Average | $\frac{1}{N}\sum_{i=1}^{N}\bm{v_{i}}$
Max | $max(\bm{v_{i}}),\ \ \forall i\in[1,N]$
FC-Max | $max(\bm{v_{i}}\bm{W_{pool}}+\bm{b_{pool}}),\ \ \forall i\in[1,N]$
Table 1:  Formulation of different graph pooling operators. Here, $\bm{v_{i}}$
denotes the atom representation and $max$ is an elementwise max operator.
$\bm{W_{pool}}$ and $\bm{b_{pool}}$ are learnable pooling parameters. Please
refer to Sec. IV.3 for more details.

###  IV.1 Effect of different pooling operators

In this section, we explore some of the existing graph pooling operators -
Average, Max and FC-max. In Average pooling, the graph embedding is defined as
the mean of all the node embeddings. Similarly, for Max pooling, it is defined
as the elementwise max of all the node embeddings. Hamilton $et\,\,al.$ [33]
proposed FC-max, a learnable pooling operator, which transforms node
embeddings using a fully-connected neural network and then applies Max pooling
operator. The formulation of these pooling operators are defined in Table 1.
To analyze the effect of pooling, we perform the experiment where we use
alumina trigonal polymorph data for training and predict performance on other
alumina polymorphs. The results are shown in Table 2. We observe that FC-max
performs the best among all of the pooling operators. And hence, we use FC-max
as the pooling operator in all further analysis of our LCGCN architecture.

Pooling Operator | Trigonal | Orthorhombic (Pbca) | Orthorhombic (Pna21)
---|---|---|---
Average | 0.0009 | 0.0349 | 0.1796
Max | 0.0007 | 0.4256 | 0.0789
FC-max | 0.0007 | 0.0112 | 0.047
Table 2:  MAE/atom of LCGCN model for different pooling operators on test data of 240 atom alumina orthorhombic (Pbca) and Orthorhombic (Pna21) polymorph trained on trigonal data. Please refer to Sec. IV.1 for details. | CGCNN | LCGCN
---|---|---
Phase |  Error/atom |  Relative Error |  Error/atom |  Relative Error
Trigonal |  -0.0466 |  0 |  -0.0018 |  0
Orthorhombic (Pbca) |  -0.0208 |  0.0258 |  0.0001 |  0.0019
Orthorhombic (Pna21) |  -0.0302 |  0.0164 |  0.0031 |  0.0049
Table 3:  Error/atom (eV/atom) and Relative Error (w.r.t. trigonal alumina
phase) of LCGCN and CGCNN model on test dataset of different polymorphs using
models trained on setups as defined in Sec. IV.2. We observe that LCGCN
consistently outperforms CGCNN for all alumina phase equilibrium predictions.

###  IV.2 Prediction of energy at equilibrium

Xie $et\,\,al.$ [13] train their CGCNN model to predict formation energy per
atom $\left(\Delta E_{f}\right)$ using a large number of equilibrium crystal
structures obtained from the Materials Project [5]. Since $\left(\Delta
E_{f}\right)$ for a large number of crystals can span across a wide range
(approximately 8.5 eV/atom), such a model will not be able to accurately
predict the small energy differences across different polymorphs of the same
composition. However, in this work we train LCGCN model to predict the PES of
various polymorphs of alumina through training on distorted input structures
obtained from MD simulations. The results are shown in table 3. Due to this
training data, we expect our model to perform better in predicting the energy
differences of the alumina polymorphs. The errors of predictions from the two
approaches are shown in table 3. We report both ’signed’ error and relative
error. For relative error, we choose the energy of the trigonal phase as the
reference value. As expected, LCGCN significantly outperforms the previous
CGCNN model  suggesting that LCGCN is well-equipped than CGCNN to predict
small energy changes typically seen in distortions. In the next upcoming
subsections, we will discuss the PES predictions of alumina across supercell
sizes, temperature and polymorphs.

Figure 2:  Plotting test MAE/atom of LCGCN model in predicting potential
energy of 120 and 480 atoms systems of alumina trigonal polymorph with
increasing size of training dataset. Please refer to Sec. IV.3 for more
details.

###  IV.3 Prediction of potential energy

The LCGCN model was used to predict the potential energy of the trigonal
alumina polymorph of different supercell sizes – namely 120 and 480 atoms. The
data to train, validate and test was obtained from MD simualtions as described
in Sec. III.1 under a temperature ramp of 300K-500K. We quantify the error by
reporting mean absolute error per atom (MAE/atom) (eV). The results for the
two different supercells are shown in Fig. 2. As expected, increasing the
amount of training data reduces the error. We observe that with 0.6 times the
total alumina data (i.e. 6000 samples) for training, our model is able to
achieve an MAE/atom which is close to the accuracy of DFT (1 meV/atom).

Training | Validation |  Testing |
---|---|---|---
120 atoms |  480 atoms |  480 atoms |  480 atoms | MAE/atom
1000 |  2400 |  600 |  2000 | 0.0017
3000 |  2400 |  600 |  2000 | 0.0017
6000 |  2400 |  600 |  2000 | 0.0011
10000 |  2400 |  600 |  2000 | 0.0011
Table 4:  MAE/atom of LCGCN model on test data of 480 atoms with increasing training dataset of 120 atoms. Please refer to Sec. IV.4 for details. Training | Validation |  Testing |
---|---|---|---
120 atoms |  480 atoms |  120 atoms |  480 atoms |  480 atoms | MAE/atom
2400 |  0 |  600 |  0 |  2000 | 0.0023
2400 |  800 |  600 |  200 |  2000 | 0.0014
2400 |  1600 |  600 |  400 |  2000 | 0.0013
2400 |  2400 |  600 |  600 |  2000 | 0.0011
Table 5:  MAE/atom of LCGCN model on test data of 480 atoms with increasing
training data of 480 atom configurations. Please refer to Sec. IV.4 for
further details.

###  IV.4 Transferable LCGCN models across scale

In many scenarios, the property of interest can only be obtained from
performing large scale atomistic and quantum mechanical simulations. However,
the scaling of such simulations with system size can be non-linear. For
example in case of DFT, the cpu time of a simulation increases as $O(N^{3})$
where $N$ is the number of atoms that limits the system size. Hence, an ML
model that is transferable from smaller scale systems to larger scale systems
are of importance in the field. Firstly, we note that the cpu time for
inference using our model increases as $O(N)$ with the number of nodes $N$. We
perform a time analysis of the LCGCN model. For this, we use a pre-trained
LCGCN model to predict potential energy for a set of 100 configurations of 120
atom alumina trigonal polymorph. Next, we repeat this for 240, 360, 480 and
600 atoms supercells. The average cpu time for each setup is depicted in Fig.
3. We observe that LCGCN is significantly faster than DFT and MD simulations,
and in contrast to MD and DFT the LCGCN scales linearly with the system size
and easily parallelized [35, 36].

Figure 3:  Time analysis of the LCGCN code with respect to the CPU run, as a
function of number of atoms.

Next, we demonstrate how increasing the training data of the 120 atoms
supercell of alumina trigonal polymorph can reduce the test error of
predicting the potential energy of the 480 atom supercell configurations. For
this, we consider a setup where there are a fixed 2400 configurations of 480
atoms for model training. Also, we have 10000 configurations of 120 atoms that
we progressively add to the training set. As observed from Table 4, an
increase in the number of 120 atom configurations in training data leads to a
reduction in the test error for the 480 atom configuration. After adding all
the 10000 configurations of 120 atoms in training data, the error in MAE/atom
for 480 atom configuration is 0.0011 eV/atom. This is shown in Fig. 4a. In a
second set of experiments, the data of 120 atoms is fixed in training set,
while the number of 480 atoms is increased in training. The results are
reported in Table 5 and depicted in Fig. 4b. We find that an increase in the
480 atom configuration in training set leads to significant improvements in
model performance on predicting potential energy for 480 atoms. From these set
of experiments we can conclude that although it is not trivial for LCGCN to
generalize, although with comparitively higher error to predict PES of the
larger supercell (480 atoms) even when these larger supercell configurations
are not included in the training data (1st row of Table 5). However, this
error goes down significantly even if a small amount of the 480 atoms
configurations are included into the training data (4th row of Table 5).

Figure 4:  (a) Test MAE/atom of LCGCN model in predicting potential energy of
480 atoms of alumina trigonal polymorph as a function of increasing training
dataset of 120 atoms. Refer Sec. IV.4 for details. (b) Test MAE/atom of LCGCN
model in predicting potential energy of 480 atoms of alumina trigonal
polymorph as a function of increasing number of 480 atom configurations in
training dataset. Refer to Sec. IV.4 for details.

###  IV.5 Transferable LCGCN models across environmental conditions

Envrionmental conditions, in particular temperature, play an important role in
determining the dynamics of the materials and the resultant phase stability
and material properties  [1, 37]. Higher temperature molecular dynamics has
been employed as a tool in different problems to sample the configurations
rapidly and to predict the resultant low energy configuration and their
corresponding properties [38]. In this section, we analyze LCGCN’s
transferability in predicting potential energies of configurations of trigonal
alumina across different environmental conditions. We specifically focus on
variation of temperature and test the effectiveness of the model in three
temperature values for alumina - 300K, 500K and 800K. We pretrain a LCGCN
model using all the data (10,000 configurations) from one temperature and
predict potential energy at a different temperature. The results are shown in
Table 6.

In $\alpha$-Al2O3, two types of Al-O bonds (1.87 angstrom and 1.98 angstrom)
are present at 0 K. But under classical MD simulation, these two bonds vary
drastically at a given temperature T. At low temperature (say 300 K), bond-
bond distribution function shows two peaks (at 1.85 and 2.03 angstrom) while
at high temperatures (500 K $\&$ 800 K), the peak corresponding to longer
(weaker) Al-O bond (2.03 angstrom peak in 300 K configuration) gets
disappeared, see Fig. 5a. It means that weakly bonded Al and O got separated
much compared to strongly bonded Al and O due to thermal stress.

Next, we observe the effect of mixing the data across two different
temperatures (resulting in a total of 20,000 configurations split into
training, validation and test set) in predicting energies for data at an
unseen temperatures. The results are shown in Table 7. We see when the model
is fed with a mixture of two temperatures data for training the prediction
error reduces significantly for test data at unseen temperatures.
Specifically, we observe a significant decrease in the MAE for prediction at
800K, when data for both 300K and 500K are mixed. In order to delineate if the
improvement in model prediction is due to training data volume, or due to data
from two different temperatures, we train a model with 5000 configurations of
300K and 500K each and test the model on data from 800K. The resulting
MAE/atom ($4^{th}$ row of Table 7) is still significantly lower than
predictions at 800K using only 300K or 500K data for training. This suggests
that the model is able to learn better when it is trained with heterogeneous
data from two different temperatures.

Temp. |  300K |  500K |  800K
---|---|---|---
300K |  0.0003 |  0.0212 |  0.0766
500K |  0.0151 |  0.0004 |  0.0254
800K |  0.0372 |  0.0192 |  0.0007
Table 6:  MAE/atom of LCGCN model on test data of 240 atom alumina trigonal polymorph at three different temperatures. Row temperature indicates that model was trained using the data of that particular temperature and the column temperature denotes that the trained model was used to predict potential energy at that temperature. Please refer to Sec. IV.5 for more details. Mixture |  Target |  MAE/atom
---|---|---
300K & 500K |  800K |  0.02
300K & 800K |  500K |  0.0048
500K & 800K |  300K |  0.0035
300K & 500K(5000 each) |  800K |  0.0169
Table 7:  MAE/atom of LCGCN model on test data of 240 atom alumina trigonal polymorph at a target temperature using a mixture of data for two different temperatures in training data. Please refer to Sec. IV.5 for details. Figure 5:  (a) Al-O bond distribution of $\alpha$-Al2O3 at 300 K , 500 K and 800 K temperatures.(b)Al-O bond distribution of Alumina polymorphs. (c) Density of energy distribution of Alumina polymorphs. Phase | Orthorhombic (Pbca) | Trigonal | Orthorhombic (Pna21)
---|---|---|---
Orthorhombic (Pbca) | 0.0007 | 0.0265 | 0.037
Trigonal | 0.0112 | 0.0007 | 0.047
Orthorhombic (Pna21) | 0.9689 | 0.1542 | 0.0006
Table 8:  MAE/atom of LCGCN model on test data of 240 atom alumina for three different polymorphs. The row phase denotes that the model was trained using the data for that particular polymorph and column phase indicates that the trained model was used to predict the potential energy for that polymorph. Refer to Sec. IV.6 for more details. Mixture | MAE/atom
---|---
Trigonal + Orthorhombic (Pna21)(5000 each) | 0.1804
Trigonal + Orthorhombic (Pna21) | 0.1312
Table 9:  MAE/atom of LCGCN model on test data of 240 atom alumina othorhombic
(Pbca) polymorph using a mixture of data for different alumina polymorphs in
training data. Please refer to Sec. IV.6 for details.

###  IV.6 Transferable LCGCN models across different Polymorphs

One of the major challenges in materials and chemistry arises from the
difficulty in transfering information, model and knowledge across different
crystal structures. A model that can learn properties such as potential
energy, from equilibrium and distorted configurations of some crystal
structures and predict for new crystal structures with minimal data would be a
very useful tool in accelerating materials discovery. As a preliminary step
toward this goal, in this work, we study the transferability of LCGCN models
across different polymorphs of alumina. We consider three phases of alumina -
orthorhombic (Pbca and Pna21) and trigonal. For the first set of experiments,
we predict the potential energy of the different polymorphs based on the
training data from another polymorph (aka zero shot learning). Similar to
temperature study, here we use all the data pertaining to a polymorph (10,000
configurations) to train the model. The results are shown in Table 8.

To understand the trends in the transferability of LCGCN across different
polymorphs, we plotted the $Al-O$ bond distribution and configurational
energies for the three alumina polymorphs as depicted in Fig. 5b respectively.
We observe that, the more the overlap between $Al-O$ bond distribution of
different polymorphs, the better transferability, and hence lower the error in
LCGCN’s prediction. For example, the overlap between the Orthorhombic (Pbca)
and Trigonal polymorphs is greater than that between the Orthorhombic (Pbca)
and Orthorhombic (Pna21) polymorphs. As seen in Table 8, our LCGCN model, when
trained on Orthorhombic (Pbca) data, predicts the Trigonal polymorph with a
lower error compared to when it is trained on Orthorhombic (Pna21). This is
likely due to the fact that the LCGCN model’s accuracy is highly dependent on
bond length, and transferability is better when predicting systems with
similar bond lengths.

In order to understand how the model accuracy improves when trained with
heterogeneous data, we train the model with data from multiple polymorphs to
predict potential energy of orthorhombic (Pbca) polymorph. We choose
orthorhombic (Pbca) for our analysis because the error in PES prediction was
the highest as shown in Table 8. The results for this analysis is shown in
Table 9. We observe that if we restrict the training data to a total of 10,000
configurations (each polymorph has 5,000 data points), the error in prediction
is about 0.1804 eV/atom suggesting that training on diverse dataset can
improve the transferability of the model. However, if we incorporate all the
data from the two polymorphs (20,000 configurations in total), the model error
further reduces to 0.1312 eV/atom suggesting that the transferability of the
model can be improved by increasing both the volume and variety of training
data.

##  V Conclusion

In this work, we proposed LCGCN, a graph neural network framework to directly
predict high-dimensional potential energy surface for alumina
$\left(Al_{2}O_{3}\right)$, a base material for various adsorbent and
catalytic applications, from atomic configurations. This approach enables to
bypass calculating the derivatives of the PES which are often the
computational bottlenecks of existing ab-initio methods.

We further demonstrated the transferability of LCGCN across scales,
environmental conditions and different alumina polymorphs. Finally, We
validated that, with appropriate conditioning, LCGCN can obtain accuracies
that are close to DFT (1 meV/atom). Therefore, we anticipate that LCGCN will
significantly accelerate the atomistic studies of other materials with its
scalability, versatility, and excellent potential energy surface prediction
accuracy.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author (J.B. and A.K.S) upon reasonable request.

## ACKNOWLEDGMENTS

The authors would like to acknowledge Shell, JNCASR and IISc for their
support.

## AUTHOR CONTRIBUTIONS

J.B, U.V.W and P.T conceived the study and planned the research. A.K.S. and
N.K. performed the MD simulations to generate the dataset. S.S , M. K, S. R
and A.K.S implemented and optimized the LCGCN model. All authors discussed the
results and their implications and contributed to the writing of the article.

## COMPETING INTERESTS

The authors declare no competing interests.

## References

  * Sagotra et al. [2019] Arun K Sagotra, Dewei Chu, and Claudio Cazorla.  Influence of lattice dynamics on lithium-ion conductivity: A first-principles study.  _Physical Review Materials_ , 3(3):035405, 2019\.
  * Eisa et al. [2022] Mohamed Eisa, Dovile Ragauskaite, Sushil Adhikari, Federico Bella, and Jonas Baltrusaitis.  Role and responsibility of sustainable chemistry and engineering in providing safe and sufficient nitrogen fertilizer supply at turbulent times.  _ACS Sustainable Chemistry & Engineering_, 10(28):8997–9001, 2022.
  * Science and [US] National Science and Technology Council (US).  _Materials genome initiative for global competitiveness_.  Executive Office of the President, National Science and Technology Council, 2011.
  * Kumar et al. [2018] Narendra Kumar, Padmini Rajagopalan, Praveen Pankajakshan, Arnab Bhattacharyya, Suchismita Sanyal, Janakiraman Balachandran, and Umesh V Waghmare.  Machine learning constrained with dimensional analysis and scaling laws: simple, transferable, and interpretable models of materials from small datasets.  _Chemistry of Materials_ , 31(2):314–321, 2018\.
  * Jain et al. [2013] Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, and Kristin a. Persson.  The Materials Project: A materials genome approach to accelerating materials innovation.  _APL Materials_ , 1(1):011002, 2013.  ISSN 2166532X.  doi: 10.1063/1.4812323.  URL http://link.aip.org/link/AMPADS/v1/i1/p011002/s1&Agg=doi.
  * Curtarolo et al. [2012] Stefano Curtarolo, Wahyu Setyawan, Gus LW Hart, Michal Jahnatek, Roman V Chepulskii, Richard H Taylor, Shidong Wang, Junkai Xue, Kesong Yang, Ohad Levy, et al.  Aflow: an automatic framework for high-throughput materials discovery.  _Computational Materials Science_ , 58:218–226, 2012.
  * [7] Computational materials database and repositories-nomad.  URL https://nomad-coe.eu/.
  * Mannodi-Kanakkithodi et al. [2016] Arun Mannodi-Kanakkithodi, Gregory M. Treich, Tran Doan Huan, Rui Ma, Mattewos Tefferi, Yang Cao, Gregory A. Sotzing, and Rampi Ramprasad.  Rational Co-Design of Polymer Dielectrics for Energy Storage.  _Adv. Mater. (Weinheim, Ger.)_ , 28(30):6277–6291, August 2016.  ISSN 0935-9648.  doi: 10.1002/adma.201600377.
  * Saal et al. [2013] James E. Saal, Scott Kirklin, Muratahan Aykol, Bryce Meredig, and C. Wolverton.  Materials Design and Discovery with High-Throughput Density Functional Theory: The Open Quantum Materials Database (OQMD).  _JOM_ , 65(11):1501–1509, November 2013.  ISSN 1047-4838, 1543-1851.  doi: 10.1007/s11837-013-0755-4.
  * Pilania et al. [2013] Ghanshyam Pilania, Chenchen Wang, Xun Jiang, Sanguthevar Rajasekaran, and Ramamurthy Ramprasad.  Accelerating materials property predictions using machine learning.  _Scientific reports_ , 3:2810, 2013.
  * Ward et al. [2016] Logan Ward, Ankit Agrawal, Alok Choudhary, and Christopher Wolverton.  A general-purpose machine learning framework for predicting properties of inorganic materials.  _npj Computational Materials_ , 2:16028, 2016.
  * Kim et al. [2016] Chiho Kim, Ghanshyam Pilania, and Ramamurthy Ramprasad.  From organized high-throughput data to phenomenological theory using machine learning: the example of dielectric breakdown.  _Chemistry of Materials_ , 28(5):1304–1311, 2016\.
  * Xie and Grossman [2018] Tian Xie and Jeffrey C. Grossman.  Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties.  _Phys. Rev. Lett._ , 120:145301, Apr 2018.  doi: 10.1103/PhysRevLett.120.145301.  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.145301.
  * Bruna et al. [2014] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  In _International Conference on Learning Representations (ICLR)_ , 2014.
  * Bronstein et al. [2017] Michael M. Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst.  Geometric deep learning: Going beyond euclidean data.  _IEEE Signal Process. Mag._ , 2017.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations (ICLR)_ , 2017.
  * Duvenaud et al. [2015] David K Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alan Aspuru-Guzik, and Ryan P Adams.  Convolutional networks on graphs for learning molecular fingerprints.  In _Advances in Neural Information Processing Systems (NIPS) 28_ , pages 2224–2232. Curran Associates, Inc., 2015.
  * Kearnes et al. [2016] Steven Kearnes, Kevin McCloskey, Marc Berndl, Vijay Pande, and Patrick Riley.  Molecular graph convolutions: moving beyond fingerprints.  _Journal of Computer-Aided Molecular Design (CAMD)_ , 30(8):595–608, 2016.
  * Gilmer et al. [2017] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural message passing for quantum chemistry.  In _Proceedings of the 34th International Conference on Machine Learning (ICML)_ , pages 1263–1272, 2017.
  * Altae-Tran et al. [2017] Han Altae-Tran, Bharath Ramsundar, Aneesh S. Pappu, and Vijay Pande.  Low data drug discovery with one-shot learning.  _ACS Central Science_ , 3(4):283–293, 2017.
  * Sanyal et al. [2018] Soumya Sanyal, Janakiraman Balachandran, Naganand Yadati, Abhishek Kumar, Padmini Rajagopalan, Suchismita Sanyal, and Partha Talukdar.  MT-CGCNN: Integrating crystal graph convolutional neural network with multitask learning for material property prediction.  _arXiv preprint arXiv:1811.05660_ , 2018.
  * Deringer and Csányi [2017] Volker L. Deringer and Gábor Csányi.  Machine learning based interatomic potential for amorphous carbon.  _Phys. Rev. B_ , 95:094203, Mar 2017.  doi: 10.1103/PhysRevB.95.094203.  URL https://link.aps.org/doi/10.1103/PhysRevB.95.094203.
  * Rowe et al. [2018] Patrick Rowe, Gábor Csányi, Dario Alfè, and Angelos Michaelides.  Development of a machine learning potential for graphene.  _Phys. Rev. B_ , 97:054303, Feb 2018.  doi: 10.1103/PhysRevB.97.054303.  URL https://link.aps.org/doi/10.1103/PhysRevB.97.054303.
  * Deringer et al. [2018] Volker L. Deringer, Chris J. Pickard, and Gábor Csányi.  Data-driven learning of total and local energies in elemental boron.  _Phys. Rev. Lett._ , 120:156001, Apr 2018.  doi: 10.1103/PhysRevLett.120.156001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.156001.
  * Veit et al. [2019] Max Veit, Sandeep Kumar Jain, Satyanarayana Bonakala, Indranil Rudra, Detlef Hohl, and Gábor Csányi.  Equation of state of fluid methane from first principles with machine learning potentials.  _Journal of Chemical Theory and Computation_ , 15(4):2574–2586, 2019.  doi: 10.1021/acs.jctc.8b01242.  URL https://doi.org/10.1021/acs.jctc.8b01242.  PMID: 30794393.
  * Bartók et al. [2010] Albert P. Bartók, Mike C. Payne, Risi Kondor, and Gábor Csányi.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  _Phys. Rev. Lett._ , 104:136403, Apr 2010.  doi: 10.1103/PhysRevLett.104.136403.  URL https://link.aps.org/doi/10.1103/PhysRevLett.104.136403.
  * Hlavay and Polyák [2005] József Hlavay and Klára Polyák.  Determination of surface properties of iron hydroxide-coated alumina adsorbent prepared for removal of arsenic from drinking water.  _Journal of Colloid and Interface Science_ , 284(1):71–77, 2005.
  * Gupta et al. [2011] Vinod K Gupta, Shilpi Agarwal, and Tawfik A Saleh.  Synthesis and characterization of alumina-coated carbon nanotubes and their application for lead removal.  _Journal of hazardous materials_ , 185(1):17–23, 2011.
  * Pines and Haag [1960] Herman Pines and Werner O Haag.  Alumina: catalyst and support. i. alumina, its intrinsic acidity and catalytic activity1.  _Journal of the American Chemical Society_ , 82(10):2471–2483, 1960.
  * Trueba and Trasatti [2005] Monica Trueba and Stefano P Trasatti.  $\gamma$-alumina as a support for catalysts: a review of fundamental aspects.  _European journal of inorganic chemistry_ , 2005(17):3393–3403, 2005.
  * Plimpton [1995] Steve Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  _Journal of Computational Physics_ , 117(1):1 – 19, 1995.  ISSN 0021-9991.  doi: https://doi.org/10.1006/jcph.1995.1039.  URL http://www.sciencedirect.com/science/article/pii/S002199918571039X.
  * Salles et al. [2016] N. Salles, O. Politano, E. Amzallag, and R. Tétot.  Molecular dynamics study of high-pressure alumina polymorphs with a tight-binding variable-charge model.  _Computational Materials Science_ , 111:181 – 189, 2016\.  ISSN 0927-0256.  doi: https://doi.org/10.1016/j.commatsci.2015.09.017.  URL http://www.sciencedirect.com/science/article/pii/S0927025615005868.
  * Hamilton et al. [2017] William L. Hamilton, Rex Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In _NIPS_ , 2017.
  * Rumelhart et al. [1986] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning Representations by Back-propagating Errors.  _Nature_ , 323(6088):533–536, 1986.  doi: 10.1038/323533a0.  URL http://www.nature.com/articles/323533a0.
  * Smith et al. [2017] Justin S Smith, Olexandr Isayev, and Adrian E Roitberg.  Ani-1: an extensible neural network potential with dft accuracy at force field computational cost.  _Chemical science_ , 8(4):3192–3203, 2017.
  * Behler et al. [2008] Jörg Behler, Roman Martoňák, Davide Donadio, and Michele Parrinello.  Pressure-induced phase transitions in silicon studied by neural network-based metadynamics simulations.  _physica status solidi (b)_ , 245(12):2618–2629, 2008.
  * Klarbring and Simak [2018] Johan Klarbring and Sergei I Simak.  Phase stability of dynamically disordered solids from first principles.  _Physical review letters_ , 121(22):225702, 2018\.
  * Behler and Parrinello [2007] Jörg Behler and Michele Parrinello.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  _Physical review letters_ , 98(14):146401, 2007\.