# Constraining the Surface Curvature of an Anisotropic Neutron Star

H. C. Das1,2<EMAIL_ADDRESS>Jeet Amrit Pattnaik3 S. K. Patra1,2
1Institute of Physics, Sachivalaya Marg, Bhubaneswar 751005, India 2Homi
Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai
400094, India 3Department of Physics, Siksha ${}^{\prime}O^{\prime}$
Anusandhan, Deemed to be University, Bhubaneswar -751030, India

###### Abstract

The anisotropy of pressure arises due to the various complex phenomena that
happen inside the neutron star (NS). In this study, we calculate the degree of
anisotropy inside the NS using the scalar pressure anisotropy model.
Macroscopic properties such as mass, radius, compactness, redshift, tidal
deformability, the moment of inertia, and surface curvature (SC) are computed
for the anisotropic NS with the equation of states spanning from relativistic
to nonrelativistic cases. The variation of SC as the functions of the above-
mentioned quantities are computed by changing the degree of anisotropy.
Pressure anisotropy has significant effects on the magnitude of SC. The
relations between the canonical SC$-\Lambda$ and SC$-\bar{I}$ are studied.
From the GW170817 tidal deformability data, we constraints the magnitude of SC
are found to be
SC${}_{1.4}(10^{14})=3.44_{-1.0}^{+0.4},2.85_{-1.20}^{+0.62},2.52_{-1.02}^{+0.61}$
for $\lambda_{\rm BL}=0.0,1.0$, and $2.0$ respectively.

##  I Introduction

One of the most difficult tasks is to examine the internal structure of
neutron stars (NS), which requires knowledge from many different branches of
physics. Due to its intricate internal structure and powerful gravitational
pull, we haven’t yet developed a comprehensive theoretical understanding of
this object. Several exotic processes are happening inside the NS, which makes
the system more complex. In addition, we consider the pressure anisotropy
inside the star, which is a realistic phenomenon. There are different sources
of pressure anisotropy inside the NS, such as a strong magnetic field [1, 2,
3, 4, 5, 6, 7, 8], pion condensation [9], phase transitions [9], core
crystallization [10], and a superfluid core [11, 12, 13]. For details
regarding the sources of anisotropy and its significance, see Ref. [14].

In literature, few anisotropic models have been developed, such as Bowers-
Liang (BL) [15], Horvat et al. [16], and Cosenza et al. [17] models. The BL
model is based on the assumption that (i) the anisotropy quadratically
disappears at the origin, (ii) it varies non-linearly on radial pressure, and
(iii) the anisotropy is gravitationally driven. According to the hypothesis by
Horvat et al. [16], the anisotropy results from the quasi-local equation as
described in Ref. [16]. According to several research, the limits of the
anisotropic parameter are $-2\leq\lambda_{\rm BL}\leq+2$ for BL model [18],
and $-2\leq\lambda_{\rm H}\leq+2$ in the Horvat model [19]. In this scenario,
the BL model is taken into account to calculate the various NS properties.

Numerous research clarified the impacts of anisotropic pressure on the
macroscopic characteristics of compact objects, such as their mass, radius,
the moment of inertia, tidal deformability, and non-radial oscillation [15,
16, 17, 20, 21, 22, 23, 24, 18, 25]. The magnitudes of macroscopic
characteristics often rise when the magnitude of the anisotropy parameter
increases, and vice versa. In this study, we intend to determine the surface
curvature (SC) as a function of different NS observables for an anisotropic
NS. It is an important quantity to measure the curvature at the surface of the
star. Also, the SC of the NS is approximately $10^{14}$ orders of magnitude
more than our solar counterpart [26]. In the near future, we may be able to
measure such curvature using the Shapiro delay measurement technique [27].

Here, we explore the relations between SC-tidal deformability and SC-moment of
inertia. Universal/approximate relations play a great role to determine the
quantity that is not accessible to detect if there exists a relation between
them. For example, the moment of inertia-tidal Love number and quadrupole
($I-$Love$-Q$) relation [28]. If we detect the Love number, we can calculate
the other two using that relation. Different types of relations exist in the
literature [28, 29, 30, 31, 32, 33]. In our previous calculations, we
calculated the $I-$Love$-C$ (where $C$ is the compactness of a star) for
different degrees of anisotropy [30]. We have put a limit on the moment of
inertia of the anisotropic NS using various observational data. In this work,
we want to find the relations between SC$-\Lambda$ and SC$-\bar{I}$ of the
anisotropic NS. However, the approximate relations are different than the
universal relations. Sometimes approximate relations are helpful in
determining the behavior of the observables. An approximate relation is used
when it is not feasible or practical to obtain the complete and accurate
nature of the data. For example, in data mining, it may be acceptable to have
a relation with missing values or errors if the goal is to identify patterns
or trends in the data. In the present scenario, the relations between SC,
tidal deformability, and the moment of inertia become weaker in case of the
maximum-mass star. However, the relations are strong between them in the case
of canonical stars, which will be very helpful in future aspects.

Hence, to explore those properties of the NS, one needs EOSs which can be
taken either from the relativistic mean-field model (RMF) or the Skyrme-
Hartree-Fock (SHF). In this calculation, we have taken RMF, SHF, and density-
dependent (DD-RMF) EOSs as described in our previous study [30]. Here, we
choose the unified EOSs to explore the anisotropic NS properties. The EOSs
which we have taken in this calculation produce the mass of the NS $\sim 2\
M_{\odot}$. The manuscript is categorized as follows. After a brief
introduction, we provide the formalism to calculate the pressure anisotropy
inside the star. The expressions and results for SC, the moment of inertia,
and tidal deformability are given in Sec-II. The numerical results are given
in Sec-VI. Finally, we enumerate the summary and conclusion in Sec-VII.

##  II Anisotropic configurations

For an isotropic star with the perfect fluid condition, the Maxwell stress-
energy tensor is defined as [34]

$\displaystyle T_{\mu\nu}=({\cal{E}}+P)u_{\mu}u_{\nu}+Pg_{\mu\nu},$ (1)

where ${\cal{E}}$, $P$, and $u_{\mu}$ are the energy density, pressure, and
4-velocity of the perfect fluid, respectively. However, various phenomena
inside the star make the pressure anisotropic in the tangential direction.
That means the magnitudes of radial pressure ($P_{r}$) differ from the
tangential pressure ($P_{t}$). For that case, the stress-energy tensor is
defined as [19, 18, 35]

$\displaystyle
T_{\mu\nu}=({\cal{E}}+P_{t})u_{\mu}u_{\nu}+(P_{r}-P_{t})k_{\mu}k_{\nu}+P_{t}g_{\mu\nu},$
(2)

where $k_{\mu}$ is the unit radial vector ($k^{\mu}k_{\mu}=1$) with
$u^{\mu}k_{\mu}=0$. The Tolman-Oppenheimer-Volkoff (TOV) equations for an
anisotropic star is defined as [19]

$\displaystyle\frac{dP_{r}}{dr}=-\frac{\left({\cal E}+P_{r}\right)\left(m+4\pi
r^{3}P_{r}\right)}{r\left(r-2m\right)}+\frac{2}{r}(P_{t}-P_{r})\,,$ (3)
$\displaystyle\frac{dm}{dr}=4\pi r^{2}{{\cal E}}\,,$ (4)

where $\sigma=P_{t}-P_{r}$ is the anisotropy parameter. The ‘$m$’ is the
enclosed mass correspond to radius $r$. To solve these TOV equations, we need
two separate EOSs for radial and tangential pressure. We use the various RMF,
DD-RMF, and SHF equation of states for radial pressure $P_{r}({\cal E})$, and
the BL model for transverse pressure $P_{t}$ [15]. The TOV equations can be
solved using the boundary conditions $r=0,m=0,P_{r}=P_{c}$, and $r=R,m=M,\
{\rm and}\ P_{r}=0$ for a particular choice of anisotropy.

###  II.1 BL Model

According to the BL model, the local anisotropy inside the star is based on
two main assumptions (i) The anisotropy quadratically disappears at the
origin, and (ii) At least a part of the anisotropy is gravitationally driven,
and it is nonlinear with the pressure. The anisotropy pressure along the
tangential direction $P_{t}$ is defined as [15]

$\displaystyle P_{t}=P_{r}+\frac{\lambda_{\rm BL}}{3}\frac{({\cal
E}+3P_{r})({\cal E}+P_{r})r^{2}}{1-2m/r}\,,$ (5)

where the factor $\lambda_{\rm BL}$ measures the degree of anisotropy in the
fluid. We assume that the value of $\lambda_{\rm BL}$ is in the range -2 to +2
following the Silva et al. [18]. This range is purely model dependent. One can
take different ranges and put constraints using different observational data.
Like, Silva et al. [18] has taken $-2<\lambda_{\rm BL}<+2$. Biswas et al. [20]
has constrained the $\lambda_{\rm BL}$ limit using GW170817 tidal
deformability data. In our previous study [30], we found that $\lambda_{\rm
BL}=+2$ satisfied the secondary mass limit of the GW190814 event. Also, the
canonical radius for an anisotropic star satisfied the NICER constraints for
$-2<\lambda_{\rm BL}<+2$. Therefore, in this study, we vary the degree of
anisotropicity in the range of $-2$ to $+2$.

##  III Surface Curvature

In this Sub-Sec., we adopt the curvature quantity from the Refs. [26, 36]. The
Kretschmann scalar (full contraction of the Riemann tensor) is defined as

$\displaystyle{\cal{K}}(r)$ $\displaystyle\equiv$
$\displaystyle\sqrt{{\cal{R}}^{\mu\nu\rho\sigma}{\cal{R}}_{\mu\nu\rho\sigma}}$
(6) $\displaystyle=$
$\displaystyle\bigg{[}(8\pi)^{2}\left\\{3{\cal{E}}^{2}(r)+3P^{2}(r)+2P(r){\cal{E}}(r)\right\\}$
$\displaystyle-\frac{128{\cal{E}}(r)m(r)}{r^{3}}+\frac{48m^{2}(r)}{r^{6}}\bigg{]}^{1/2},$

At the surface, $m\rightarrow M$ as $r\rightarrow R$. Except for the last
term, all other terms vanish outside the star because they depend on
${\cal{E}}_{\rm NS}(r)$, and $P_{\rm NS}(r)$, which are zero outside the star.
But, there is a non-vanishing component of the Riemann tensor that does not
vanish; $\tensor{{\cal R}}{{}^{1}_{010}}=-\frac{2M}{R^{3}}=-\xi$, even in the
outside of the star [26, 37]. Therefore, the Riemann tensor is a more relevant
quantity to measure the curvature of the stars. Kretschmann scalar is the
square root of the full contraction of the Riemann tensor. The vacuum value
for both $\cal{K}$ is $\frac{4\sqrt{3}M}{R^{3}}$, one can easily see from Eq.
(6). Therefore, one can take $\cal K$ reasonable measures for the curvature
within the star. The SC is defined as the ratio of curvature at the surface of
the NS ${\cal{K}}(R)$ to the curvature of the Sun ${\cal{K}}_{\odot}$, SC
$={\cal{K}}(R)/{\cal{K}}_{\odot}$. This ratio
${\cal{K}}(R)/{\cal{K}}_{\odot}\approx 10^{14}$ i.e, the NS curvature is
$10^{14}$ times more than the Sun.

##  IV Moment of Inertia

For a slowly rotating NS, the system’s equilibrium position can be obtained by
solving Einstein’s equation in the Hartle-Throne metric as [38, 39, 40]

$\displaystyle ds^{2}=$ $\displaystyle-e^{2\nu}\ dt^{2}+e^{2\lambda}\
dr+r^{2}\ (d\theta^{2}+\sin^{2}\theta d\phi^{2})$ (7)
$\displaystyle-2\omega(r)r^{2}\sin^{2}\theta\ dt\ d\phi.$

The MI of the slowly rotating anisotropic NS is calculated in Ref. [22]

$\displaystyle I$
$\displaystyle=\frac{8\pi}{3}\int_{0}^{R}\frac{r^{5}J\tilde{\omega}}{r-2M}({\cal{E}}+P)\left[1+\frac{\sigma}{{\cal{E}}+P}\right]\,dr,$
(8)

where $\tilde{\omega}=\bar{\omega}/\Omega$, where $\bar{\omega}$ is the frame
dragging angular frequency, $\bar{\omega}=\Omega-\omega(r)$. $J$ is defined as
$e^{-\nu}(1-2m/r)^{1/2}$.

Figure 1: Left: Mass-radius profiles for anisotropic NS with
$-2.0<\lambda_{\rm BL}<+2.0$ for IOPB-I EOS. Different color bands signify the
masses of the NS observed from the various pulsars, such as PSR J0348+0432
[41], PSR J0740+6620 [42], heaviest pulsars J0952-0607 [43], and GW190814
[44]. The NICER results are shown with two green boxes from two different
analyses [45, 46]. The revised NICER results are also shown for the canonical
star and $2.08\ M_{\odot}$ (red horizontal error bars) given by Miller et al.
[47]. Right: Variation of mass as a function of central density. The black
line represents the isotropic star.

##  V Tidal Deformability

The shape of the NS is deformed when it is present in the external field
($\epsilon_{ij}$) of its companion. Hence the stars develop the quadrupole
moment ($Q_{ij}$), which is linearly dependent on the tidal field and is
defined as [48, 49]

$\displaystyle Q_{ij}=-\lambda\epsilon_{ij}\,,$ (9)

where $\lambda$ is defined as the tidal deformability of a star. It has
relation to the dimensionless tidal Love number $k_{2}$ as
$\lambda=\frac{2}{3}k_{2}R^{5}$, where $R$ is the radius of the star. The
dimensionless tidal deformability is defined as
$\Lambda=\lambda/M^{5}=2k_{2}/3C^{5}$, where $C$ is the compactness of the
star defined as $M/R$.

To determine $k_{2}$, we use the linear perturbation in the Throne and
Campolattaro metric [50]. We have solved the Einstein equation and obtained
the following second-order differential equation for the anisotropic star [20]

$\displaystyle H^{{}^{\prime\prime}}$
$\displaystyle+H^{{}^{\prime}}\bigg{[}\frac{2}{r}+e^{\lambda}\left(\frac{2m(r)}{r^{2}}+4\pi
r(P-{\cal E})\right)\bigg{]}$ $\displaystyle+H\left[4\pi
e^{\lambda}\left(4{\cal E}+8P+\frac{{\cal E}+P}{dP_{t}/d{\cal
E}}(1+c_{s}^{2})\right)-\frac{6e^{\lambda}}{r^{2}}-{\nu^{\prime}}^{2}\right]$
$\displaystyle=0\,.$ (10)

The term $dP_{t}/d{\cal E}$ represents the change of $P_{t}$ (see Eq. (5) for
the $P_{t}$) with respect to energy density for a fixed value of $\lambda_{\rm
BL}$.

The internal and external solutions to the perturbed variable $H$ at the
star’s surface can be matched to get the tidal Love number [51, 48]. The value
of the tidal Love number can be calculated using the $y_{2}$, and compactness
parameter $C$ is defined as [48, 49, 52]

$\displaystyle k_{2}$
$\displaystyle=\,\frac{8}{5}C^{5}(1-2C)^{2}\big{[}2(y_{2}-1)C-y_{2}+2\big{]}$
$\displaystyle\times\Big{\\{}2C\big{[}4(y_{2}+1)C^{4}+2(3y_{2}-2)C^{3}-2(11y_{2}-13)C^{2}$
$\displaystyle+3(5y_{2}-8)C-3(y_{2}-2)\big{]}+3(1-2C)^{2}$
$\displaystyle\times\big{[}2(y_{2}-1)C-y_{2}+2\big{]}\log(1-2C)\Big{\\}}^{-1}\,,$
(11)

where $y_{2}$ depends on the surface value of $H$ and its derivative

$y_{2}=\frac{rH^{{}^{\prime}}}{H}\Big{|}_{R}.$ (12)

Figure 2: Left: SC as a function of mass for anisotropic NS with
$-2.0<\lambda_{\rm BL}<+2.0$ for IOPB-I equation of state. Right: SC as a
function of radius. Figure 3: Left: SC as a function of red-shift for
anisotropic NS with $-2.0<\lambda_{\rm BL}<+2.0$ for IOPB-I equation of state.
Right: SC as a function of compactness.

##  VI Results and Discussions

###  VI.1 Mass-Radius relations

For IOPB-I equation of state, the mass-radius profiles of the anisotropic NS
are solved for various values of BL, as illustrated in Fig. 1. The magnitude
of the maximum masses and their corresponding radii increases with increasing
positive values of $\lambda_{\rm BL}$, and vice versa. Different observational
data, including x-ray, NICER, and GW, may constrain the degree of anisotropy
within the NS. Recently, it is discovered that PSR J0952-0607, the fastest and
heaviest Galactic NS in the Milky Way disc, has a mass $M=2.35\pm
0.17\,M_{\odot}$ [43]. Additionally, we can put this restriction to control
the degree of anisotropy.

The GW190814 event sparked a debate about whether the secondary component is
the lightest black hole or the heaviest neutron star [44]. To explain this
behavior, several models have already been proposed in the literature [53, 54,
55, 56, 57]. However, Roupas et al. Ref. [53] asserted that the secondary
component might be an anisotropic NS. Hence, we investigated the mass-radius
diagram with the secondary component mass limit $M=2.50-2.67\ M_{\odot}$ to
see if it replicated the limit for anisotropic stars in the BL model. The mass
$2.50–2.67\ M_{\odot}$ can be reproduced for $\lambda_{\rm BL}=1.8–2$,
although those values of $\lambda_{\rm BL}$ do not adhere to the revised NICER
limits [47].

Figure 4: SC, $y_{2}$, and $k_{2}$ with the variation of $0.0<\lambda_{\rm
BL}<+2.0$ for canonical anisotropic NS for G3 (soft), IOPB-I (moderate stiff),
and NL3 (stiff) EOSs.

Figure 5: Left: SC as a function of $\Lambda$ for anisotropic NS with
$0.0<\lambda_{\rm BL}<+2.0$ for IOPB-I equation of state. The error bars are
the observational constraints given by LIGO/Virgo events GW170817 (NS-NS
merger) [58] and GW190814 (assuming BH-NS merger) [44]

. Right: SC as a function of $\bar{I}$ for $-2.0<\lambda_{\rm BL}<+2.0$.

###  VI.2 Surface Curvatures

We calculate the SC with the IOPB-I equation of state by changing the value of
$\lambda_{\rm BL}$, which is shown in Fig. 2. The magnitude of SC increases
both with the mass of the star as well as with the positive value of
$\lambda_{\rm BL}$ and vice-versa. This is due to the fact that the variation
of positive values $\lambda_{\rm BL}$ enhances the magnitude of the tangential
pressure of the star, which supports large gravitational mass. In the previous
section, we explained that the SC is directly proportional to the mass of the
star. Therefore, with increasing the mass of the star, the SC increases.
However, the negative values of $\lambda_{\rm BL}$ follow the opposite trend,
which is clearly seen in the figure. Also, the radial variation of SC is
depicted in the same figure. The variation in SC mainly follows the $1/R^{3}$
trend. The magnitude of the radius increases towards the surface, which
decreases SC. However, it enhances its magnitude near the center. Various
observational bands are the data from which one can constrain the magnitude of
curvature, as well as the amount of anisotropy inside the star [30]. Recently,
the gravity-field curvature has been measured using the three conjugated atom
interferometers [27]. In the future, we may be able to measure the SC of the
NS by either using interferometers or Shapiro delay measurements.

The parametric dependence of the variation of SC has been discussed in our
previous work [36]. The softer EOS, such as G3, provides more SC in comparison
to IOPB-I (see Fig. 5 of Ref. [36]). This is due to the fact that the
magnitude of both $M$ and $R$ decreases for softer EOS in comparison to
stiffer ones. Hence, the quantity $M/R^{3}$ for soft EOS is more. Also, it
increases more for dark matter (DM) admixed star due to the fact that with the
addition of DM, the EOS becomes softer [59, 36, 60, 61, 57].

In a similar way, we show the variations of SC with red-shift
$\left(Z_{s}=\frac{1}{\sqrt{1-2C}}-1\right)$ and compactness in Fig. 3 for
IOPB-I equation of state as a representative case. The magnitude of SC
increases with increasing both for $Z_{s}$ and $C$. Since the redshift is a
function of compactness only, therefore, the trend of changing SC is the same
as compactness. However, the magnitude of both $Z_{s}$ and $C$ increases with
increasing the value of $\lambda_{\rm BL}$. From different observational data,
we can put constrain on the degree of anisotropy inside the star. Until now,
just one value of $Z_{s}=0.35$ from the analysis of stacked bursts in Ref.
[62] has been reported in the low-mass x-ray binary EXO 0748-676, which is
also discarded by the subsequent observation [63].

The variation of SC, $y_{2}$, and $k_{2}$ by varying only positive values of
$\lambda_{\rm BL}$ are shown in Fig. 4 for the canonical star. In this case,
we take three different types of EOSs, such as G3 (soft) [64], IOPB-I
(moderately stiff) [64], and NL3 (stiff) [65], for representative cases. The
magnitude of SC is larger for G3 in comparison to IOPB-I and NL3. However, the
magnitude of $y_{2}$ and $k_{2}$ are almost the same for the three EOSs with
the variation of $\lambda_{\rm BL}$. This simply indicates the model-
independent nature of $y_{2}$ and $k_{2}$ for different values of
$\lambda_{\rm BL}$ for the canonical anisotropic NS.

The variation of SC both for dimensionless tidal deformability and the moment
of inertia for the IOPB-I equation of state with different values of
$\lambda_{\rm BL}$ is shown in Fig. 5. We observe that the magnitude of SC
decreases with increasing either $\Lambda$ or $\bar{I}$. Also, the relations
seem to be correlated for a fixed $\lambda_{\rm BL}$ as well as the mass. The
relation becomes tighter in the case of SC$-\bar{I}$ as compared to
SC$-\Lambda$. One can put a direct constraint on the SC from the various
observational data if any relation exists between them. Therefore, in the next
sub-section, we try to get the functional relations between them by changing
the degree of anisotropy.

Figure 6: Relations between SC${}_{1.4}-\Lambda_{1.4}$ (left), and
SC${}_{1.4}-\bar{I}_{1.4}$ (right) by varying anisotropy parameter
$\lambda_{\rm BL}$. The lower panels in each figure are for the residual of
the fitting using Eq. (14). Table 1: The fitting coefficients are listed for
SC${}_{1.4}-\Lambda_{1.4}$, and SC${}_{1.4}-\bar{I}_{1.4}$ relations with
$\lambda_{\rm BL}=0.0,1.0,$ and $2.0$ respectively.

SC${}_{1.4}-\Lambda_{1.4}$ |  SC${}_{1.4}-\bar{I}_{1.4}$
---|---
$\lambda_{\rm BL}=$ | 0.0 | 1.0 | 2.0 | $\lambda_{\rm BL}=$ | 0.0 | 1.0 | 2.0
$a_{0}=$ | $4.0908$ | $3.9029$ | $3.5629$ | $b_{0}=$ | $12.3498$ | $10.8343$ | $9.7856$
$a_{1}(10^{-3})=$ | $-3.6762$ | $-6.3974$ | $-6.3987$ | $b_{1}=$ | $-1.2453$ | $-0.9840$ | $-0.8173$
$a_{2}(10^{-6})=$ | $1.4463$ | $4.57661$ | $4.9067$ | $b_{2}=$ | $0.0358$ | $0.0254$ | $0.0194$

###  VI.3 Relations between SC, $\Lambda$, and $\bar{I}$ for a canonical star

Various universal relations have already been proposed, such as $I-$Love$-Q$,
$I-$Love-$C$, etc., both for static as well as rotating NS. In this case, we
want to calculate the approximate/universal functional relations between
SC$-\Lambda$ and SC$-I$ for an anisotropic star.

The canonical value for SC$-\Lambda$ with spans of almost 60 EOSs are
calculated and shown in left panels of Fig. 6 for the anisotropic star cases
for $\lambda_{\rm BL}=0.0,1.0$, and $2.0$. In this calculation, we take only
positive values of $\lambda_{\rm BL}$. This is because the higher negative
values of $\lambda_{\rm BL}$ give unphysical solutions of $\Lambda$, and lower
negative values don’t change the magnitude of $\Lambda$ significantly, as
discussed in Refs. [20, 30]. The calculated data are fitted with the
approximate formula using the least-square method [30]

${\rm SC}_{1.4}=\sum_{n=0}^{2}a_{n}(\Lambda_{1.4})^{n},$ (13)

where $a_{n}$ is the fitting coefficient given in Table 1. The lower panel of
the figure represents the residuals for the fit calculated using the formula
[30]

$\displaystyle\Delta{\rm SC}=\frac{|{\rm SC}-{\rm SC}_{\rm fit}|}{{\rm
SC}_{\rm fit}}.$ (14)

Another relation between SC${}_{1.4}-\bar{I}_{1.4}$ is also shown in the right
panels of Fig. 6 for $\lambda_{\rm BL}=0.0,1.0$, and $2.0$. We use the same
fitting formula as given in Eq. (13), but the coefficients are different as
given in Table 1 for $b_{n}$. We also checked for this calculation for the
maximum mass star; however, the relationships were found to be weak in
comparison to the canonical mass. Hence, in this study, we mainly focus on the
relations between SC$-\Lambda$ and SC$-\bar{I}$ for the canonical star.

From these approximate relations, one can put constrain on the SC of the NS
with the help of observational data of $\Lambda$ and $\bar{I}$. The tidal
deformability limit given by GW170817 from the merger of two binary NS is
$\Lambda_{1.4}=190_{-120}^{+390}$. We obtain the value of SC for the canonical
star
SC${}_{1.4}(10^{14})=3.44_{-1.0}^{+0.4},2.85_{-1.2}^{+0.62},2.52_{-1.02}^{+0.61}$
for $\lambda_{\rm BL}=0.0,1.0$, and $2.0$ respectively. To our knowledge, till
now, we don’t have any observational data for the moment of inertia of the NS.
However, various theoretical studies have been dedicated to constraining the
value of the moment of inertia using various observational data. From these
limits, one can also constrain the magnitude of surface curvature.

##  VII Summary and Conclusions

In this study, we have explored various properties of anisotropic NS within
the scalar anisotropic BL model. The magnitude of various macroscopic
properties changes with the degree of anisotropy. This is due to the fact that
the pressure difference between the radial and transverse components has an
additional contribution to balance the hydrostatic equilibrium, and it is
purely model-dependent. Any anisotropic models should satisfy the different
perfect fluid conditions as mentioned in Refs. [30].

We have calculated various macroscopic properties such as mass, radius, tidal
deformability, the moment of inertia, and surface curvature for a variety of
EOSs spanning from relativistic to non-relativistic cases. We observed that
the magnitude of maximum mass and its corresponding radius increases
(decreases) with the positive (negative) values of $\lambda_{\rm BL}$.
However, the magnitude of canonical tidal deformability value decreases with
positive values of $\lambda_{\rm BL}$ and vice-versa. Again, the effect of
$\lambda_{\rm BL}$ on $y_{2}$ and $k_{2}$ are almost similar for all those
considered sets, showing the model-independent nature for the anisotropic NS.
The surface curvature of the NS is a crucial quantity which is $10^{14}$
orders of magnitude more than the sun. Here, we mainly focus on the variation
of SC with some NS properties by varying degrees of anisotropy inside the NS.
The SC plays a great role to describe the curvature made by the star in the
space-time metric. Using different observational data, we can put stringent
constrain on its magnitude if we can find any relationship between known
observables such as mass, radius, and tidal deformability. In another way, we
can put constraints on the direct observational data such as Shapiro delay
measurement. Therefore, in this study, we have explored the relationships
between SC and other quantities for the anisotropic NSs.

The magnitude of SC as a function of mass increases (decreases) with the
positive (negative) value of $\lambda_{\rm BL}$. This is because the SC
directly depends on mass. However, the radial variation depends on
($1/R^{3}$), which means, at the core, the SC is more in comparison to the
surface of the star. Since almost all mass of the NS is mainly concentrated in
the core. In a similar fashion, the SC follows the same trends for the
variation of redshift and compactness with the degree of anisotropy. The SC as
a function of either tidal deformability or the moment of inertia decreases by
varying anisotropicity. From these relations, we observed that there might be
some relationship between the SC$-\Lambda$ and SC$-\bar{I}$.

We have found the approximate relations between SC$-\Lambda$ and SC$-\bar{I}$
for the canonical star. However, the universality becomes weaker for the
maximum mass of the star. From the relation between
SC${}_{1.4}-\Lambda_{1.4}$, we put a constraint on the canonical SC using the
GW170817 tidal deformability data. The obtained values are
SC${}_{1.4}(10^{14})=3.44_{-1.0}^{+0.4},2.85_{-1.2}^{+0.62},\ {\rm and}\
2.52_{-1.02}^{+0.61}$ for $\lambda_{\rm BL}=0.0,1.0$, and $2.0$ respectively.
However, till now, we don’t have any observational data on the dimensionless
moment of inertia. Maybe in the future, by the observation of double pulsars,
we can find its value and put stringent constraints on the surface curvature
as well as the degree of anisotropy.

## Acknowledgement

JAP acknowledges to Institute Of Physics (IOP), Bhubaneswar, for providing the
computer facilities during the work. SERB (Project Nos. CRG/2019/002691)
partly reinforces this work.

## References

  * Yazadjiev [2012] S. S. Yazadjiev, Phys. Rev. D 85, 044030 (2012).
  * Cardall  _et al._ [2001] C. Y. Cardall, M. Prakash, and J. M. Lattimer, The Astrophysical Journal 554, 322 (2001).
  * Ioka and Sasaki [2004] K. Ioka and M. Sasaki, The Astrophysical Journal 600, 296 (2004).
  * Ciolfi  _et al._ [2010] R. Ciolfi, V. Ferrari, and L. Gualtieri, Monthly Notices of the Royal Astronomical Society 406, 2540 (2010).
  * Ciolfi and Rezzolla [2013] R. Ciolfi and L. Rezzolla, Monthly Notices of the Royal Astronomical Society: Letters 435, L43 (2013).
  * Frieben and Rezzolla [2012] J. Frieben and L. Rezzolla, Monthly Notices of the Royal Astronomical Society 427, 3406 (2012).
  * Pili  _et al._ [2014] A. G. Pili, N. Bucciantini, and L. Del Zanna, Monthly Notices of the Royal Astronomical Society 439, 3541 (2014).
  * Bucciantini  _et al._ [2015] N. Bucciantini, A. G. Pili, and L. Del Zanna, Monthly Notices of the Royal Astronomical Society 447, 3278 (2015).
  * Carter and Langlois [1998] B. Carter and D. Langlois, Nuclear Physics B 531, 478 (1998).
  * Nelmes and Piette [2012] S. Nelmes and B. M. A. G. Piette, Phys. Rev. D 85, 123004 (2012).
  * Kippenhahn Rudolf [1990] W. A. Kippenhahn Rudolf,  _Stellar Structure and Evolution_, Vol. XVI (Springer Berlin, Heidelberg, 1990).
  * Glendenning [1997] N. K. Glendenning,  _Compact stars: Nuclear physics, particle physics, and general relativity_ (Springer-Verlag New York, 1997).
  * Heiselberg and Hjorth-Jensen [2000] H. Heiselberg and M. Hjorth-Jensen, Physics Reports 328, 237 (2000).
  * Herrera and Santos [1997] L. Herrera and N. Santos, Physics Reports 286, 53 (1997).
  * Bowers and Liang [1974] R. L. Bowers and E. P. T. Liang, Astrophys. J. 188, 657 (1974).
  * Horvat  _et al._ [2010] D. Horvat, S. Ilijić, and A. Marunović, Classical and Quantum Gravity 28, 025009 (2010).
  * Cosenza  _et al._ [1981] M. Cosenza, L. Herrera, M. Esculpi, and L. Witten, J Math Phys (NY) 22, 118 (1981).
  * Silva  _et al._ [2015] H. O. Silva, C. F. B. Macedo, E. Berti, and L. C. B. Crispino, Classical and Quantum Gravity 32, 145008 (2015).
  * Doneva and Yazadjiev [2012] D. D. Doneva and S. S. Yazadjiev, Phys. Rev. D 85, 124023 (2012).
  * Biswas and Bose [2019] B. Biswas and S. Bose, Phys. Rev. D 99, 104002 (2019).
  * Sulaksono [2015] A. Sulaksono, International Journal of Modern Physics E 24, 1550007 (2015).
  * Rahmansyah, A.  _et al._ [2020] Rahmansyah, A., Sulaksono, A., Wahidin, A. B., and Setiawan, A. M., Eur. Phys. J. C 80, 769 (2020).
  * Rahmansyah  _et al._ [2020] A. Rahmansyah, A. Sulaksono, A. B. Wahidin, and A. M. Setiawan, The European Physical Journal C 80, 769 (2020).
  * Rahmansyah and Sulaksono [2021] A. Rahmansyah and A. Sulaksono, Phys. Rev. C 104, 065805 (2021).
  * Deb  _et al._ [2021] D. Deb, B. Mukhopadhyay, and F. Weber, The Astrophysical Journal 922, 149 (2021).
  * Ekşi  _et al._ [2014] K. Y. Ekşi, C. Güngör, and M. M. Türkoğlu, Phys. Rev. D 89, 063003 (2014).
  * Rosi  _et al._ [2015] G. Rosi, L. Cacciapuoti, F. Sorrentino, _et al._ , Phys. Rev. Lett. 114, 013001 (2015).
  * Yagi and Yunes [2015] K. Yagi and N. Yunes, Phys. Rev. D 91, 123008 (2015).
  * Haskell  _et al._ [2014] B. Haskell, R. Ciolfi, F. Pannarale, and L. Rezzolla, Mon. Not. R. Astron. Soc. Lett. 438, L71 (2014).
  * Das [2022] H. C. Das, Phys. Rev. D 106, 103518 (2022).
  * Landry and Kumar [2018] P. Landry and B. Kumar, APJ 868, L22 (2018).
  * Kumar and Landry [2019] B. Kumar and P. Landry, Phys. Rev. D 99, 123026 (2019).
  * Pradhan  _et al._ [2023] B. K. Pradhan, A. Vijaykumar, and D. Chatterjee, Phys. Rev. D 107, 023010 (2023).
  * Walecka [1974] J. Walecka, Ann. Phys. 83, 491 (1974).
  * Estevez-Delgado and Estevez-Delgado [2018] G. Estevez-Delgado and J. Estevez-Delgado, The European Physical Journal C 78, 673 (2018).
  * Das  _et al._ [2021a] H. C. Das, A. Kumar, B. Kumar, S. K. Biswal, and S. K. Patra, JCAP 2021 (01), 007.
  * He  _et al._ [2015] X.-T. He, F. J. Fattoyev, B.-A. Li, and W. G. Newton, Phys. Rev. C 91, 015810 (2015).
  * Hartle [1967] J. B. Hartle, The Astrophysical Journal 150, 1005 (1967).
  * Hartle and Thorne [1968] J. B. Hartle and K. S. Thorne, The Astrophysical Journal 153, 807 (1968).
  * Hartle [1973] J. B. Hartle, Astrophysics and Space Science 24, 385 (1973).
  * Antoniadis  _et al._ [2013] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, _et al._ , Science 340, 1233232 (2013).
  * Fonseca  _et al._ [2021] E. Fonseca, H. T. Cromartie, T. T. Pennucci, _et al._ , The Astrophysical Journal Letters 915, L12 (2021).
  * Romani  _et al._ [2022] R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, The Astrophysical Journal Letters 934, L17 (2022).
  * Abbott  _et al._ [2020] R. Abbott, T. D. Abbott, S. Abraham, _et al._ , The Astrophysical Journal 896, L44 (2020).
  * Miller  _et al._ [2019] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, _et al._ , APJ 887, L24 (2019).
  * Riley  _et al._ [2019] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, _et al._ , APJ 887, L21 (2019).
  * Miller  _et al._ [2021] M. C. Miller, F. K. Lamb, A. J. Dittmann, _et al._ , The Astrophysical Journal Letters 918, L28 (2021).
  * Hinderer [2008] T. Hinderer, The Astrophysical Journal 677, 1216 (2008).
  * Hinderer [2009] T. Hinderer, The Astrophysical Journal 697, 964 (2009).
  * Thorne and Campolattaro [1967] K. S. Thorne and A. Campolattaro, APJ 149, 591 (1967).
  * Damour and Nagar [2009] T. Damour and A. Nagar, Phys. Rev. D 80, 084035 (2009).
  * Das  _et al._ [2021b] H. C. Das, A. Kumar, B. Kumar, S. K. Biswal, and S. K. Patra, International Journal of Modern Physics E 30, 2150088 (2021b).
  * Roupas [2021] Z. Roupas, Astrophysics and Space Science 366, 9 (2021).
  * Fattoyev  _et al._ [2020] F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and B. Reed, Phys. Rev. C 102, 065805 (2020).
  * Huang  _et al._ [2020] K. Huang, J. Hu, Y. Zhang, and H. Shen, The Astrophysical Journal 904, 39 (2020).
  * Lim  _et al._ [2021] Y. Lim, A. Bhattacharya, J. W. Holt, and D. Pati, Phys. Rev. C 104, L032802 (2021).
  * Das  _et al._ [2021c] H. C. Das, A. Kumar, and S. K. Patra, Phys. Rev. D 104, 063028 (2021c).
  * Abbott  _et al._ [2017] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * Das  _et al._ [2020] H. C. Das, A. Kumar, B. Kumar, _et al._ , MNRAS 495, 4893 (2020).
  * Das  _et al._ [2021d] H. C. Das, A. Kumar, and S. K. Patra, MNRAS 507, 4053 (2021d).
  * Das  _et al._ [2021e] H. C. Das, A. Kumar, S. K. Biswal, and S. K. Patra, Phys. Rev. D 104, 123006 (2021e).
  * Cottam  _et al._ [2002] J. Cottam, F. Paerels, and M. Mendez, Nature 420, 51 (2002).
  * Cottam  _et al._ [2008] J. Cottam, F. Paerels, M. Mendez, _et al._ , The Astrophysical Journal 672, 504 (2008).
  * Parmar  _et al._ [2022] V. Parmar, H. C. Das, A. Kumar, A. Kumar, M. K. Sharma, P. Arumugam, and S. K. Patra, Phys. Rev. D 106, 023031 (2022).
  * Lalazissis  _et al._ [1997] G. A. Lalazissis, J. König, and P. Ring, Phys. Rev. C 55, 540 (1997).


# Calibrated Uncertainty for Molecular Property Prediction using Ensembles of
Message Passing Neural Networks

Jonas Busk1,∗ Peter Bjørn Jørgensen1 Arghya Bhowmik1 Mikkel N. Schmidt2 Ole
Winther2,3,4 Tejs Vegge1

( 1Department for Energy Conversion and Storage, Technical University of
Denmark, Lyngby, Denmark
2Department of Applied Mathematics and Computer Science, Technical University
of Denmark, Lyngby, Denmark
3Center for Genomic Medicine, Rigshospitalet, Copenhagen University Hospital,
Copenhagen, Denmark
4Bioinformatics Centre, Department of Biology, University of Copenhagen,
Copenhagen, Denmark
∗Corresponding author: Jonas Busk<EMAIL_ADDRESS>)

###### Abstract

Data-driven methods based on machine learning have the potential to accelerate
analysis of atomic structures. However, machine learning models can produce
overconfident predictions and it is therefore crucial to detect and handle
uncertainty carefully. Here, we extend a message passing neural network
designed specifically for predicting properties of molecules and materials
with a calibrated probabilistic predictive distribution. The method presented
in this paper differs from the previous work by considering both aleatoric and
epistemic uncertainty in a unified framework, and by re-calibrating the
predictive distribution on unseen data. Through computer experiments, we show
that our approach results in accurate models for predicting molecular
formation energies with calibrated uncertainty in and out of the training data
distribution on two public molecular benchmark datasets, QM9 and PC9. The
proposed method provides a general framework for training and evaluating
neural network ensemble models that are able to produce accurate predictions
of properties of molecules with calibrated uncertainty.

##  1 Introduction

Autonomous high-throughput computational analysis of atomic structures has the
potential to speed up the discovery of novel materials and chemical reactions
dramatically with applications in a wide range of research areas like
biotechnology and conversion and storage of renewable energy. This process can
be enabled and accelerated by data-driven methods based on machine learning
that are less computationally demanding than traditional quantum mechanical
methods such as density functional theory (DFT) [1]. In recent years, graph-
based models such as message passing neural networks (MPNNs), that operate on
atomic structures represented as graphs, have shown impressive capabilities at
predicting properties of molecules and materials with high accuracy [2].
However, deep neural networks are known to produce overconfident predictions
[3], especially outside the training data distribution, which can lead to sub-
optimal or incorrect results. Because chemical space is too vast to represent
in any training dataset [4, 5], in an autonomous high-throughput setting it is
crucial to detect and handle predictive uncertainty, for example by falling
back to more accurate but computationally demanding methods like DFT [6].
Consequently, predictive algorithms that express probabilistic uncertainty
about predictions can help identify problematic instances and thereby enable
the design of new robust workflows [7] and applications in computational
materials research.

When quantifying uncertainty it is often useful to distinguish between
_epistemic_ and _aleatoric_ uncertainty [8, 9]. Epistemic uncertainty arises
from the model’s inability to fit the data distribution and can in principle
be reduced by observing more data or improving the model. Aleatoric
uncertainty on the other hand comes from inherent noise in the data and can
therefore not be reduced by observing more data. When the aleatoric
uncertainty is constant across all observations it is called _homoscedastic_
aleatoric uncertainty and is often not modelled explicitly. If the aleatoric
uncertainty depends on the input, and thus varies across the data
distribution, it is called _heteroscedastic_ aleatoric uncertainty and can be
estimated from the data by explicitly including it in the model. Thus
epistemic uncertainty is important for understanding when predictions can be
trusted and aleatoric uncertainty captures noise in the data. Consequently,
both types of uncertainty are necessary to obtain a complete picture of the
predictive uncertainty and provide important information for interpretability
and decision making.

In this work, we extend a message passing neural network regression model
designed specifically for predicting properties of molecules and materials
[10] with a probabilistic predictive distribution and consider a deep ensemble
of models [11] to express the aleatoric and epistemic uncertainty about
predictions of molecular formation energies. The predicted uncertainty is re-
calibrated to fit the error distribution on unseen data to address model
overconfidence from training and the expected improved accuracy from using an
ensemble approximation. We show through computer experiments that our approach
results in well calibrated models on two public benchmark datasets for
molecular property prediction, QM9 [12] and PC9 [13]. The main contribution of
this paper is a complete framework for training and evaluating neural network
models that are able to produce accurate predictions of properties of
molecules with calibrated aleatoric and epistemic uncertainty estimates.

Uncertainty quantification for property prediction of atomic structures with
graph neural networks has received increasing interest recently. Scalia et al.
[14] evaluated and compared scalable uncertainty estimation methods based on
graph neural networks for molecular property prediction and found that deep
ensembles [11] and bootstrapping consistently outperformed Monte Carlo Dropout
[15] on multiple public benchmark datasets in terms of error and uncertainty
calibration. Tran et al. [16] highlight the importance of predictive
uncertainty in materials screening applications and review methods for
uncertainty quantification and procedures for evaluating the quality of
uncertainty estimates including accuracy, calibration and sharpness. Nigam et
al. [17] provide an extensive overview of different sources of uncertainty in
molecular property prediction in the context of drug discovery, many of which
are also relevant in materials research, and describe the importance and
perspectives of having good uncertainty estimates in data driven decision
making. Related work has studied the use of Gaussian Process Regression models
for molecular property prediction [18] and molecular dynamics [19]. The method
presented in this paper differs from the previous work by considering both
aleatoric and epistemic uncertainty in a unified framework, and by re-
calibrating the uncertainty estimates to obtain more accurate uncertainties on
unseen data.

The rest of the paper is organised as follows. The proposed method is
introduced in Section 2 and experiments and results are presented in Section
3. The main findings and perspectives are discussed in Section 4 and we
conclude in Section 5.

##  2 Method

###  2.1 Message passing neural network model

In general, a message passing neural network (MPNN), as described in [2],
operates on a graph structure $g$ with node features $x_{v}$ and edge features
$e_{vw}$, where $v$ and $w$ denote vertices in the graph. A forward pass
through the neural network consists of two phases: i) a message passing phase
with $T$ interaction steps where messages are passed along the edges of the
graph to update the internal graph embedding, and ii) a readout phase where an
output value $\hat{y}$ is computed from the final graph embedding.

We base our work on the SchNet with edge updates MPNN model, which was
previously introduced by the authors [10]. This model is in turn is based on
the SchNet model, that was designed specifically for predicting properties of
molecules and materials [20]. We refer the reader to the cited literature for
specific details about the neural network architecture. It is worth noting
that the uncertainty quantification method proposed in the following sections
does not depend on the particular choice of neural network base model and can
thus be adapted to use other models based on the specific application.

###  2.2 Extended model with predictive uncertainty

To capture both epistemic and heteroscedastic aleatoric uncertainty, we extend
the MPNN described in the previous section by constructing a deep ensemble of
networks [11] (without adversarial training) in the following way. Given a
regression task with a training dataset
$\mathcal{D}=\\{g_{n},y_{n}\\}_{n=1}^{N}$ consisting of $N$ datapoints with
real-valued targets $y\in\mathbb{R}$, we consider an ensemble of $M$ neural
network models with parameters $\\{\theta_{m}\\}_{m=1}^{M}$, each with
probabilistic predictive distribution:

$p_{\theta}(y|g)=\mathcal{N}\big{(}\mu_{\theta}(g),\sigma_{\theta}^{2}(g)\big{)}\,,$
(1)

assuming a normal distribution of errors. Each network is constructed with two
outputs corresponding to the predicted mean $\mu_{\theta}(g)$ and variance
$\sigma_{\theta}^{2}(g)$, where the latter represents the predicted
heteroscedastic aleatoric uncertainty [21]. The predicted variance is
constrained to be positive by passing the second network output through the
softplus function, $\log(1+\exp(\cdot))$, and adding a small minimum variance
for numerical stability (e.g. $10^{-6}$).

###  2.3 Model training procedure

Each network of the ensemble is initialized with random parameters and trained
individually (on the same training set) using stochastic gradient descent to
minimise the negative log likelihood (NLL) loss:

$\displaystyle\text{NLL}(\theta)$
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}-\log p_{\theta}(y_{n}|g_{n})$ (2)
$\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\frac{1}{2}\bigg{(}\frac{1}{\sigma_{\theta}^{2}(g_{n})}\underbrace{\big{(}y_{n}-\mu_{\theta}(g_{n})\big{)}^{2}}_{\text{squared
error}}+\log\sigma_{\theta}^{2}(g_{n})+\underbrace{\log
2\pi}_{\text{constant}}\bigg{)}\,.$ (3)

The last term in (3) is constant since it does not depend on $\mu_{\theta}(g)$
or $\sigma_{\theta}^{2}(g)$ and can be ignored for the purpose of training the
model. Notice how for constant variance (homoscedastic uncertainty) this is
equivalent to minimising the standard mean squared error (MSE) loss. Notice
also how the predicted uncertainty acts as learned loss attenuation by letting
examples with high predicted uncertainty have smaller impact on the total
loss, while the $\log\sigma_{\theta}^{2}$ term discourages large uncertainties
[9].

In practice, we found that training directly with NLL loss can be unstable
because of interactions between the mean and variance output in the loss
function. To mitigate this, we initially train the mean output of the network
before introducing the variance terms by interpolating from MSE to NLL loss:

$\mathcal{L}(\theta)=\lambda\,\text{MSE}(\theta)+(1-\lambda)\,\text{NLL}(\theta)\,,$
(4)

where $\lambda$ is set to 1 for a number of warmup steps and then decreased
linearly from 1 to 0 over a number of interpolation steps. The resulting loss
function is quite natural since the NLL loss includes the squared error term
(see eq. 3) and as a result we found that model training becomes more stable
and robust to outliers in the training data. Additional measures exist to
promote the stability of training variance networks [21, 22, 23], but we found
the method above to be sufficient in this case.

###  2.4 Ensemble mixture

To produce the ensemble predictive distribution $p_{*}(y|g)$ and capture
epistemic uncertainty, we follow the approach of [11] and make an ensemble
approximation by combining the predictions of the $M$ individual models as a
uniformly-weighted mixture of normal distributions:

$p_{*}(y|g)=\frac{1}{M}\sum_{m=1}^{M}p_{\theta_{m}}(y|g)\,,$ (5)

whose mean $\mu_{*}(g)$ and variance $\sigma_{*}^{2}(g)$ are given by the
following expressions:

$\displaystyle\mu_{*}(g)$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\mu_{\theta_{m}}(g)\,,$ (6)
$\displaystyle\sigma_{*}^{2}(g)$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\big{(}\sigma_{\theta_{m}}^{2}(g)+\mu_{\theta_{m}}^{2}(g)\big{)}-\mu_{*}^{2}(g)$
(7)
$\displaystyle=\underbrace{\frac{1}{M}\sum_{m=1}^{M}\sigma_{\theta_{m}}^{2}(g)}_{\text{aleatoric
uncertainty}}+\underbrace{\frac{1}{M}\sum_{m}\mu_{\theta_{m}}^{2}(g)-\mu_{*}^{2}(g)}_{\text{epistemic
uncertainty}}\,.$ (8)

The variance of the ensemble predictive distribution represents the total
predicted uncertainty and can be decomposed into aleatoric and epistemic
uncertainty as shown in eq. 8 above.

###  2.5 Uncertainty calibration

Training with NLL loss can result in overfitting of the uncertainty to the
training data once a good fit for the mean is found [21]. This can lead to
overconfident predictions as we expect the model to have higher error, and
thus higher uncertainty, when presented with unseen data. On the other hand,
using an ensemble of models can reduce the overall error, and should thus lead
to lower uncertainty. This is not reflected in the ensemble variance (eq. 8)
which is strictly higher than the average of the individual variances.
Furthermore, there is nothing in the training procedure which ensures that the
epistemic uncertainty fits the error distribution. Thus there is a need to re-
calibrate the uncertainty to unseen data. This can be achieved by applying a
calibration function, that maps the predictive distribution to a calibrated
distribution.

Several approaches to re-calibrate regression models have been proposed in the
literature [24, 25, 26]. A simple, yet robust, method is to simply scale the
predicted uncertainty estimates by a scaling factor $s_{n}^{2}$ optimised to
fit the error distribution on unseen data [25]. Uncertainty scaling has the
advantage that it does not modify the mean prediction and the calibrated
predictive distribution remains a normal distribution:

$p_{s^{2}}(y_{n}|g_{n})=\mathcal{N}\big{(}\mu_{*}(g_{n}),s_{n}^{2}\sigma_{*}^{2}(g_{n})\big{)}\,.$
(9)

The solution proposed in [25] is to scale all uncertainty estimates by a
constant scaling factor optimised to minimise the NLL on a held out
calibration dataset, but it was also suggested to use non-linear methods. We
found that with our method better calibration results were achieved by
applying a non-linear function to scale the uncertainty estimates.
Specifically, to obtain the scaled uncertainty estimates we apply an isotonic
regression model111Specifically, we use the implementation of isotonic
regression available from the scikit-learn Python package [27]:
sklearn.isotonic.IsotonicRegression. $f_{\phi}(\cdot)$ to fit the empirical
squared errors $(y_{n}-\mu_{\theta}(g_{n}))^{2}$ on a held out calibration
dataset:

$s_{n}^{2}\sigma_{*}^{2}(g_{n})=f_{\phi}(\sigma_{*}^{2}(g_{n}))\,.$ (10)

The isotonic regression results in a monotonic increasing scaling function and
thus has the desired property of being non-linear while maintaining the
overall ordering of the uncertainty estimates.

###  2.6 Evaluating the uncertainty calibration

Evaluating uncertainty estimates in regression tasks is not straight forward
as the true uncertainties are generally unknown [24, 25, 26]. In the
classification domain, a classifier is said to be calibrated, if the predicted
class probability corresponds to the probability that the prediction is
correct. In other words, the classifier is expected to correctly predict its
error. In a similar way, calibration for regression can be defined as
correctly predicting the expected error for each data point [25].

The NLL is the standard metric for evaluating the quality of probabilistic
models by measuring the probability of observing the data given the predicted
distribution. In regression, the NLL depends both on the predicted mean and
variance, but it is also useful to evaluate the predicted uncertainty on its
own. To evaluate the uncertainty calibration of a model, we can thus compare
the predicted uncertainties to the observed errors on a test dataset. To do
this, we sort examples by their predicted uncertainty, divide them into $K$
bins and in each bin $k$ compute the predicted root mean variance (RMV) and
the empirical root mean squared error (RMSE). The resulting values can be used
to produce a reliability diagram where RMSE is plotted as a function of RMV
[25]. For a well calibrated model, the plotted points should lie close to the
diagonal line corresponding to the identity function. The calibration error
can then be summarized by the expected normalized calibration error (ENCE):

$\text{ENCE}=\frac{1}{K}\sum_{k=1}^{K}\frac{|\text{RMV}_{k}-\text{RMSE}_{k}|}{\text{RMV}_{k}}\,.$
(11)

This quantity is analogous to the expected calibration error (ECE) often used
in classification [25].

Additionally, it is useful to assess the dispersion of the uncertainty
estimates by computing the coefficient of variation (CV) [25]:

$\text{CV}=\frac{\sqrt{\tfrac{1}{N}\sum_{n=1}^{N}(\sigma_{*}(g_{n})-\mu_{\sigma_{*}})^{2}}}{\mu_{\sigma_{*}}}\,,$
(12)

where $\sigma_{*}(g_{n})$ is the predicted standard deviation (uncertainty) of
instance $n$, $\mu_{\sigma_{*}}=\frac{1}{N}\sum_{n=1}^{N}\sigma_{*}(g_{n})$ is
the mean predicted standard deviation and $N$ in this case iterates over the
test dataset. A high CV indicates large dispersion (heteroscedasticity) of the
uncertainty estimates and thus a high input dependence, whereas a CV of zero
indicates constant (homoscedastic) uncertainty and consequently the
uncertainty estimates are uninformative.

As a secondary diagnostic tool we can compare the predicted quantiles to the
observed quantiles averaged over the entire test dataset [24]. Visualising
these quantities in a calibration plot reveals if the predictive distribution
fits the error distribution on average. For a well calibrated model, the
plotted line should lie close to the diagonal. The quantile calibration can be
summarised by the sum of squared errors (SE) between the predicted and
observed quantiles.

##  3 Experiments and results

###  3.1 Datasets

In our experiments we consider two publicly available datasets: QM9 [12],
which is a widely used benchmark for machine learning predictions of molecular
properties, and the more recent PC9 [13], that contains a more diverse set of
molecules, but with the same general constraints as QM9.

The QM9 dataset consists of 133,885 small organic molecules in equilibrium
state with up to 9 heavy atoms (C, O, N, F) besides hydrogen. For each
molecule, the dataset contains several quantum chemical properties calculated
at the B3LYP/6-31G(2df,p) level of theory including total energy $U_{0}$,
which incorporates the vibrational zero point energy (ZPE) [12]. We
additionally compute the total energy without the ZPE, $E=U_{0}-\text{ZPE}$,
to enable comparison with PC9, that does not include $U_{0}$. The PC9 dataset
[13] consists of 99,234 molecules extracted from the PubChem database [28] by
applying the constraints of QM9 mentioned above and was found to represent a
more diverse set of molecules than QM9. PC9 includes properties calculated at
the B3LYP/6-31G(d) level of theory including total energy $E$. Structures that
appear in both datasets were identified by comparing International Chemical
Identifiers (InChI) [29] (see supplementary material A for details). We found
that 21,777 molecules from QM9 are also in PC9 and 21,619 molecules from PC9
are also in QM9.

In line with previous work, we actually consider the atomisation energies (the
energy remaining after subtracting the energies of the constituent atoms) in
our experiments, rather than the total energies. Thus in subsequent sections,
$U_{0}$ and $E$ will be used to refer to the respective atomisation energies.

###  3.2 Experimental setup

To evaluate the proposed method, we performed experiments of predicting
atomisation energies on the QM9 and PC9 datasets. In each experiment, we
trained a deep ensemble of $M=5$ message passing neural network models
extended to predict uncertainty as described in Section 2. The models were
trained individually using the same hyperparameters and data splits, but with
random parameter initialisation and random shuffling of the training data to
induce model diversity. Following previous work [10], the networks were
constructed with $T=3$ interaction steps, a cutoff distance of 5.0 Å for
constructing the molecular graphs, and an embedding size of 256 nodes in each
layer. We used the PyTorch implementation of the AdamW optimizer [30] with an
initial learning rate of 0.0001, an exponential decay learning rate scheduler,
and a weight decay coefficient of 0.01. Each model was trained for up to
3,000,000 gradient steps with 1,000,000 warmup steps training only with MSE
loss ($\lambda=1$) and loss interpolated linearly from MSE to NLL on the
subsequent 1,000,000 steps (see eq. 4). The validation set was used for early
stopping with NLL criterion and was also used as calibration set for fitting
the calibration function $f_{\phi}$ as described in Section 2.5.

Dataset | Error (eV) | Calibration
---|---|---
Train | Test | Target | MAE | RMSE | NLL | ENCE | CV
QM9 | QM9 | $U_{0}$ | 0.0094 | 0.0313 | -3.1593 | 0.0484 | 1.8939
QM9 | QM9 | $E$ | 0.0101 | 0.0342 | -3.0759 | 0.0720 | 1.7293
PC9 | PC9 | $E$ | 0.0199 | 0.0844 | -2.5956 | 0.0650 | 2.2011
QM9 | PC9 | $E$ | 0.4192 | 0.7410 | 0.8107 | 0.0220 | 0.6294
PC9 | QM9 | $E$ | 0.1165 | 0.1737 | -0.5366 | 0.0312 | 0.5597
Table 1:  Test results of ensemble models trained to predict atomisation
energy properties on the QM9 and PC9 datasets. Mean absolute error (MAE) and
root mean squared error (RMSE) are presented in electron volt (eV). The
uncertainty calibration in each experiment is summarised by the mean negative
log likelihood (NLL), expected normalised calibration error (ENCE), and
coefficient of variation (CV) of the uncertainty estimates.

###  3.3 Prediction of $U_{0}$ on QM9 with random split

In this first experiment, we trained an ensemble to predict the total energy
$U_{0}$ property of the QM9 dataset. Following previous work [2, 10, 20], we
randomly split the data into a training set of 110,000 molecules, a validation
set of 10,000 molecules, and a test set consisting of the remaining 13,885
molecules. The test set results are presented in the first row of Table 1. The
ensemble achieved a MAE = 0.0094 eV which is comparable to previous work using
a similar model [10] (MAE = 0.0105 eV).

Figure 1 shows the trade off between error and ensemble size of up to 10
models on the validation set. As expected, using a larger ensemble clearly
improves the error. In this case, reasonably low errors can be achieved with
an ensemble of $M=5$ models and not much is gained by increasing the size of
the ensemble beyond that, so we choose to use ensembles of this size
throughout our experiments. Similar behavior was observed for all experiments.

(a)

(b)

Figure 1: Trade off between error and ensemble size evaluated on the QM9
validation set when predicting the $U_{0}$ property: (a) mean absolute error
(MAE) in eV and (b) mean negative log likelihood (NLL). The models are sorted
by NLL in increasing order (best first). Reasonably low errors can be achieved
with an ensemble size of $M=5$ models.

The ensemble uncertainty estimates were re-calibrated by fitting an isotonic
regression calibration function (see Section 2.5) on the validation set and
applying it on the test set leading to an average scaling factor of 0.2965 (SD
= 0.5346) on the test set. Even though each individual model is expected to
have increased error when presented with unseen data, the ensemble
approximation significantly improved the overall error in this case resulting
in a calibration function that effectively shrinks the uncertainty of the
predictive distribution. Uncertainty calibration diagrams are presented in
Figure 2 and uncertainty evaluation metrics are summarised in the first row of
Table 1. The reliability diagram shows that in general the model assigns
higher uncertainty to instances with higher error and overall the model is
well calibrated since the predicted uncertainties correspond closely to the
observed errors on average, placing the points close to the diagonal. The
rightmost bin, representing predictions with the highest uncertainty
estimates, includes instances with relatively large errors, placing this point
far from the rest. However, the model correctly assigns high uncertainty to
these instances, thereby identifying them as problematic. The diagram also
reveals that epistemic uncertainty is relatively low in most cases, indicating
a high level of agreement among the individual models of the ensemble, and
consequently the aleatoric uncertainty is responsible for the majority of the
total uncertainty in these cases. The last bin is an exception showing high
epistemic uncertainty, corresponding to a high level of disagreement among the
individual models, indicating these molecules are out of distribution and
therefore the predictions are also more likely to have high error.
Additionally, the uncertainty estimates have a CV (eq. 12) well above 0,
indicating that the uncertainties assigned by the model have large dispersion
and thus are highly input dependent.

(a)

(b)

Figure 2: Evaluation of uncertainty on the QM9 test set when predicting the
$U_{0}$ property. The reliability diagram (a) shows observed root mean squared
error (RMSE) as a function of predicted uncertainty measured by root mean
variance (RMV) computed in bins. The error is measured by expected normalized
calibration error (ENCE). The calibration plot (b) compares predicted
quantiles and observed quantiles averaged over the test data. The error is
measured by the sum of squared errors (SE).

Learning curves for this experiment are presented in Figure 3 showing test set
metrics as a function of the amount of training data when predicting $U_{0}$
on QM9. As expected, the error metrics decrease with more training data.
Interestingly, good calibration in terms of the ENCE was obtained with
relatively small training datasets and does not vary significantly when adding
more data, while the dispersion of uncertainty estimates measured by the CV
clearly increases with the amount of training data, making the uncertainty
estimates more input dependent and thus more informative.

(a)

(b)

Figure 3: Learning curves showing test set metrics as a function of training
set size on the QM9 dataset when predicting $U_{0}$. (a) The mean absolute
error (MAE) and root mean squared error (RMSE) improve with more data as
expected. (b) The calibration in terms of expected normalised calibration
error (ENCE) does not vary significantly, while the dispersion of uncertainty
estimates measured by the coefficient of variation (CV) increases with the
amount of training data.

###  3.4 Prediction of $E$ on QM9 with random split

Complementary to the first experiment, we also trained an ensemble to predict
the total energy $E$ of the QM9 dataset using the same data split. The test
set results are presented in the second row of Table 1. The ensemble model
achieved a MAE = 0.0101 eV, which is a little higher than when predicting
$U_{0}$, suggesting predicting $E$ is harder. A similar finding was reported
in [13] using a SchNet [20] model.

The uncertainty estimates were likewise re-calibrated by fitting an isotonic
regression calibration function on the validation set and applying it on the
test set leading to an average scaling factor of 0.3116 (SD = 0.3966) on the
test set, effectively shrinking the predictive distribution similarly to the
first experiment. Also in this case we found that the model succeeds at
assigning uncertainty estimates that correlates with the error and the model
is well calibrated. Both ENCE and CV results are similar to the first
experiment. Uncertainty calibration figures for this experiment are included
in the supplementary material B1.

###  3.5 Prediction of $E$ on PC9 with random split

Next, we trained an ensemble to predict the total energy $E$ of the more
diverse PC9 dataset. The data was split randomly into a training set of 80,000
molecules, a validation set of 10,000 molecules, and a test set consisting of
the remaining 9,234 molecules. The test set results are presented in the third
row of Table 1. The ensemble model achieved a MAE = 0.0199 eV which is
approximately twice as high as when predicting $E$ on QM9. We attribute this
increase in error to PC9 representing a more diverse set of molecules, making
the task more difficult, and additionally to the smaller amount of available
training data. A similar increase in error between QM9 and PC9 was reported in
[13] using a SchNet [20] model.

The uncertainty estimates were re-calibrated by fitting a calibration function
on the validation set and applying it on the test set leading to an average
scaling factor of 0.2938 (SD = 0.6449) on the test set, effectively shrinking
the uncertainty of the predictive distribution similarly to the two previous
experiments. This model also succeeds in assigning uncertainty estimates that
correlates with the observed error and the model is well calibrated. As in the
two previous experiments, some instances among the predictions with the
highest uncertainty have relatively large errors. The ENCE is similar to the
two previous experiments while the CV is higher, indicating more dispersion of
the uncertainty estimates corresponding to the overall increase in error.
Uncertainty calibration figures for this experiment are included in the
supplementary material B2.

###  3.6 Generalisation from QM9 to PC9

In this experiment we examine the effect of applying an ensemble of models
trained on QM9 to the more diverse set of molecules included in PC9 and
especially how it affects the uncertainty estimates as we anticipate larger
errors. When constructing the data splits, we utilize the fact that the
datasets are overlapping, using the 112,108 molecules that are unique to QM9
as the training set and the 21,777 structures from QM9 that are also present
in PC9 as the validation set. Then we compute a linear correction on the
21,619 structures from PC9 that are present in QM9 to account for the
different level of theory used to calculate the energy properties. Following
[13], the linear correction was performed by fitting a Huber regression model
(coefficient = 1.0038, intercept = 1.1428) on the predicted and observed
energies. Finally, we use the remaining 77,615 molecules exclusive to PC9 as
the the test set and apply the linear correction to the predictions. The test
set results are presented in the fourth row of Table 1. The ensemble model
achieved a MAE = 0.4192 eV, which is comparable to the findings reported in
[13] using a SchNet [20] model. The relatively high error is caused primarily
by out of distribution instances, and indicates that the model has problems
generalising under domain shift, and secondly by the different level of theory
used to calculate the energies in the different datasets, which was shown to
produce large errors (see Figure 3 in [13]).

The uncertainty estimates were re-calibrated by fitting an isotonic regression
calibration function on the validation set and applying it on the test set
leading to an average scaling factor of 135.0313 (SD = 31.5991) on the test
set, reflecting the large increase in error when extrapolating to a more
diverse dataset containing out of distribution molecules. Interestingly, the
uncertainty estimates are still well calibrated as shown by a low ENCE, while
the dispersion of the uncertainty estimates is low indicated by the relatively
low CV. Uncertainty calibration figures for this experiment are included in
the supplementary material B3. The calibration plot in Figure B3 shows that
the predicted quantiles do not fit the observed quantiles averaged over the
dataset in this experiment. This is primarily because the errors are not
normally distributed in this particular case [13].

###  3.7 Generalisation from PC9 to QM9

Now going in the opposite direction, in this last experiment we examine the
effect of applying an ensemble of models trained on PC9 to the less diverse
set of molecules in QM9. Analogous to the previous experiment, we use the
77,615 molecules that are unique to PC9 as the training set and the 21,619
structures from PC9 that are also present in QM9 as the validation set.
Similarly to the previous experiment, we compute a linear correction on the
21,777 structures from QM9 that are also present in PC9 by fitting a Huber
regression model (coefficient = 0.9994, intercept = −0.6830) on the predicted
and observed energies. Finally, we use the remaining 112,108 molecules
exclusive to QM9 as the test set. The test set results are presented in the
fifth and final row of Table 1. The ensemble achieved a MAE = 0.1165 eV, which
is comparable to the findings in [13] using a SchNet [20] model. While high
compared to the experiment of predicting $E$ on QM9 above, the error is
significantly lower than the previous experiment of training on QM9 and
testing on PC9 as might be expected when going from a more diverse dataset to
an overlapping and less diverse dataset. Some of the error may be attributed
to the different level of theory used to calculate the energies in the
datasets.

The uncertainty estimates were re-calibrated by fitting a calibration function
on the validation set and applying it on the test set leading to an average
scaling factor of 6.2404 (SD = 1.4109) on the test set, which like the error
is also significantly lower than the previous experiment. Similarly to the
previous experiment, the uncertainty remains well calibrated shown by low ENCE
and the dispersion of uncertainty estimates is also relatively low in this
case. Uncertainty calibration figures for this experiment are included in the
supplementary material B4.

##  4 Discussion

In our first three experiments, we trained ensembles to predict atomisation
energies on the QM9 and PC9 benchmark datasets, respectively, with random data
splits. The results of predicting total energy $U_{0}$ on QM9 are comparable
to previous work by the authors using the same base model [10], meaning we did
not loose accuracy with our approach. We saw a small increase in the error
when predicting $E$ on QM9 which is consistent with results reported in [13].
The error when predicting $E$ on the more diverse PC9 dataset was almost twice
as high compared to QM9, which is also consistent with results reported in
[13], indicating that the increase in chemical diversity makes the predictions
task harder. In all three random split experiments, the proposed method
produced well calibrated uncertainty estimates characterised by highly
correlated average uncertainties and errors, as shown in the reliability
diagram in Figure 2 and additionally in the supplementary material B, and
summarized by low ENCE values (see Table 1). The reliability diagrams further
show that for the test examples with high error the epistemic uncertainty is
high relative to the aleatoric uncertainty, indicating high variance among the
predictions of the individual models in the ensemble. This means that the
ensemble model is good at identifying instances that are out of distribution
and therefore have high expected error, and exemplifies why it is useful to be
able to distinguish between epistemic and aleatoric uncertainty in the
predictions. The dispersion of uncertainty estimates, measured by the CV, was
high in all three random split experiments, indicating the predicted
uncertainty estimates are highly input dependent and thereby informative.

In the fourth experiment, we aimed to generalise from QM9 to the more diverse
PC9 dataset by training on QM9 and testing on molecules exclusive to PC9. The
analysis of the PC9 structures presented in [13] showed that some molecules
included in PC9 are chemically different from molecules in QM9, making this a
difficult out-of-distribution prediction task. Additionally, the properties of
the datasets where computed at different levels of theory (B3LYP/6-31G(2df,p)
in QM9 and B3LYP/6-31G(d) in PC9), which we accounted for with a linear
correction, following [13]. The error we observed in this experiment was quite
high, but comparable to what is reported in [13]. However, the uncertainty
estimates of our model were still well calibrated, meaning the model is
capable of correctly identifying the high error instances and assigning them
high uncertainty, which means the out of distribution cases can be detected.
In the fifth and final experiment, we went in the opposite direction and
trained on PC9 to predict the molecules exclusively in QM9. This should be an
easier task, since QM9 is similar to but less diverse than PC9. As expected,
the error we observed is significantly lower than in the previous experiment
and comparable to what was reported in [13]. Our method again produced
calibrated uncertainty estimates, thus on average estimating the expected
error of the predictions. Interestingly, the two generalisation experiments
resulted in the lowest ENCE of all our experiments (see Table 1) despite
having the highest errors. Additionally, in the learning curve experiment
presented in Figure 3 we also observed that good calibration was achieved even
for small training set sizes where the error is relatively high. As such there
does not appear to be a direct correlation between error and calibration. We
find this to be an interesting topic for further research.

In the proposed method, the effectiveness of the ensemble approximation,
including the quality of the epistemic uncertainty estimates, relies on
training a diverse set of models to ensure variance of predictions beyond the
training data distribution. The method proposed in this work depends on random
initialisation of network parameters and random shuffling of the training data
to induce model diversity, but other more deliberate methods exist.
Bootstrapping, i.e. re-sampling the training set with replacement, is a
popular technique for inducing diversity in ensemble models, but some evidence
suggests that this method is less appropriate for deep models as they
typically perform better with more training data [11]. We tried to apply
bootstrapping in our experiments, but did not observe any improvements in
terms of error or calibration, so we left it out for simplicity. Another more
recent approach to induce diversity is to use randomized prior functions [31],
which we consider an interesting direction for future work.

A major advantage of the proposed method is the possibility to distinguish
between epistemic and aleatoric uncertainty in the predictions. Capturing
epistemic uncertainty may be especially important in the area of atomic
structure analysis. In some other areas, where a lot of data is available, it
may be sufficient to model aleatoric uncertainty, which cannot be reduced with
more data, and reduce epistemic uncertainty with large amounts of training
data. Chemical space, however, is so vast that it is not feasible to gather
enough training data to cover the entire domain [4, 5]. Thus, identifying
cases beyond the training data distribution where the model is not expected to
perform well is more critical. In particular, distinguishing between epistemic
and aleatoric uncertainty can be utilised in a screening system for atomic
structures in the following way: If the epistemic uncertainty of a prediction
is low, the aleatoric uncertainty indicates the expected error. If, on the
other hand, the epistemic uncertainty is high, there is a high level of
disagreement in the ensemble and therefore low confidence in the prediction,
and the system can automatically fall back to a more accurate method such as
DFT [32]. The specific thresholds for decision making can be tuned depending
on the data, application and computational resources available.

##  5 Conclusion

In this work we have explored a complete framework for obtaining well
calibrated uncertainty estimates for molecular property prediction by using a
deep ensemble of message passing neural networks and by re-calibrating the
uncertainty estimates to unseen data. Our experiments on two publicly
available benchmark datasets show that the method is able to produce well
calibrated uncertainty estimates in and out of the training data distribution.
A major advantage of this approach is that the uncertainty estimates can be
decomposed into epistemic and aleatoric uncertainty, which allows for better
interpretability and provides useful information for decision making.
Additionally, the proposed method does not depend on the particular
architecture of the neural network model, and can thus easily be adapted to
use new models based on the specific application and as research in model
development continues to advance.

##  6 Acknowledgements

The authors acknowledge support from the Novo Nordisk Foundation (SURE,
NNF19OC0057822) and the European Union’s Horizon 2020 research and innovation
program under grant agreement No 957189 (BIG-MAP) and No 957213
(BATTERY2030PLUS).

## References

  * [1] Pavlo O Dral.  Quantum chemistry in the age of machine learning.  J. Phys. Chem. Lett., 11(6):2336–2347, March 2020.
  * [2] 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, volume 70 of Proceedings of Machine Learning Research, pages 1263–1272. PMLR, 2017.
  * [3] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q. Weinberger.  On calibration of modern neural networks.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 1321–1330. JMLR.org, 2017.
  * [4] Rafael Gómez-Bombarelli, Jennifer N. Wei, David Duvenaud, José Miguel Hernández-Lobato, Benjamín Sánchez-Lengeling, Dennis Sheberla, Jorge Aguilera-Iparraguirre, Timothy D. Hirzel, Ryan P. Adams, and Alán Aspuru-Guzik.  Automatic chemical design using a data-driven continuous representation of molecules.  ACS Central Science, 4(2):268–276, 2018.  PMID: 29532027.
  * [5] Pavel G Polishchuk, Timur I Madzhidov, and Alexandre Varnek.  Estimation of the size of drug-like chemical space based on gdb-17 data.  Journal of computer-aided molecular design, 27(8):675–679, 2013\.
  * [6] Rocco Peter Fornari, Murat Mesta, Johan Hjelm, Tejs Vegge, and Piotr de Silva.  Molecular engineering strategies for symmetric aqueous organic redox flow batteries.  ACS Materials Letters, 2(3):239–246, 2020.
  * [7] Felix T. Bölle, Nicolai R. Mathiesen, Alexander J. Nielsen, Tejs Vegge, Juan Maria Garcia-Lastra, and Ivano E. Castelli.  Autonomous discovery of materials for intercalation electrodes.  Batteries & Supercaps, 3(6):488–498, 2020.
  * [8] Eyke Hüllermeier and Willem Waegeman.  Aleatoric and epistemic uncertainty in machine learning: an introduction to concepts and methods.  Machine Learning, 110(3):457–506, Mar 2021.
  * [9] Alex Kendall and Yarin Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 5580–5590, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [10] Peter Bjørn Jørgensen, Karsten Wedel Jacobsen, and Mikkel Nørgaard Schmidt.  Neural message passing with edge updates for predicting properties of molecules and materials.  In Machine Learning for Molecules and Materials, Neural Information Processing Systems workshop, 2018.
  * [11] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  Advances in neural information processing systems, 30:6402–6413, 2017.
  * [12] Raghunathan Ramakrishnan, Pavlo O Dral, Matthias Rupp, and O Anatole Von Lilienfeld.  Quantum chemistry structures and properties of 134 kilo molecules.  Scientific data, 1(1):1–7, 2014.
  * [13] Marta Glavatskikh, Jules Leguy, Gilles Hunault, Thomas Cauchy, and Benoit Da Mota.  Dataset’s chemical diversity limits the generalizability of machine learning predictions.  Journal of Cheminformatics, 11, 11 2019.
  * [14] Gabriele Scalia, Colin A. Grambow, Barbara Pernici, Yi-Pei Li, and William H. Green.  Evaluating scalable uncertainty estimation methods for deep learning-based molecular property prediction.  Journal of Chemical Information and Modeling, 60(6):2697–2717, 2020\.  PMID: 32243154.
  * [15] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In Maria Florina Balcan and Kilian Q. Weinberger, editors, Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1050–1059, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * [16] Kevin Tran, Willie Neiswanger, Junwoong Yoon, Qingyang Zhang, Eric Xing, and Zachary W Ulissi.  Methods for comparing uncertainty quantifications for material property predictions.  Mach. Learn.: Sci. Technol., 1(2):025006, May 2020.
  * [17] AkshatKumar Nigam, Robert Pollice, Matthew F. D. Hurley, Riley J. Hickman, Matteo Aldeghi, Naruki Yoshikawa, Seyone Chithrananda, Vincent A. Voelz, and Alán Aspuru-Guzik.  Assigning confidence to molecular property prediction.  CoRR, abs/2102.11439, 2021.
  * [18] Félix Musil, Michael J Willatt, Mikhail A Langovoy, and Michele Ceriotti.  Fast and accurate uncertainty estimation in chemical machine learning.  J. Chem. Theory Comput., 15(2):906–915, February 2019.
  * [19] Giulio Imbalzano, Yongbin Zhuang, Venkat Kapil, Kevin Rossi, Edgar A Engel, Federico Grasselli, and Michele Ceriotti.  Uncertainty estimation for molecular dynamics and sampling.  J. Chem. Phys., 154(7):074102, February 2021.
  * [20] Kristof Schütt, Pieter-Jan Kindermans, Huziel Enoc Sauceda Felix, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller.  Schnet: A continuous-filter convolutional neural network for modeling quantum interactions.  Advances in neural information processing systems, 30:991–1001, 2017.
  * [21] D. A. Nix and A. S. Weigend.  Estimating the mean and variance of the target probability distribution.  In Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94), volume 1, pages 55–60 vol.1, 1994.
  * [22] Alexander Amini, Wilko Schwarting, Ava Soleimany, and Daniela Rus.  Deep evidential regression.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 14927–14937. Curran Associates, Inc., 2020.
  * [23] Nicki Skafte Detlefsen, Martin Jørgensen, and Søren Hauberg.  Reliable training and estimation of variance networks.  In Proceedings of 33rd Conference on Neural Information Processing Systems, 2019.  33rd Conference on Neural Information Processing Systems, NeurIPS 2019 ; Conference date: 08-12-2019 Through 14-12-2019.
  * [24] Volodymyr Kuleshov, Nathan Fenner, and Stefano Ermon.  Accurate uncertainties for deep learning using calibrated regression.  In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2796–2804. PMLR, 10–15 Jul 2018\.
  * [25] Dan Levi, Liran Gispan, Niv Giladi, and Ethan Fetaya.  Evaluating and calibrating uncertainty prediction in regression tasks.  arXiv, 2020.
  * [26] Hao Song, Tom Diethe, Meelis Kull, and Peter Flach.  Distribution calibration for regression.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 5897–5906. PMLR, 09–15 Jun 2019.
  * [27] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [28] Sunghwan Kim, Jie Chen, Tiejun Cheng, Asta Gindulyte, Jia He, Siqian He, Qingliang Li, Benjamin Shoemaker, Paul Thiessen, Bo Yu, Leonid Zaslavsky, Jian Zhang, and Evan Bolton.  Pubchem in 2021: New data content and improved web interfaces.  Nucleic Acids Research, 49, 11 2020.
  * [29] Alan McNaught.  The iupac international chemical identifier.  Chemistry international, pages 12–14, 2006.
  * [30] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization, 2019.
  * [31] Ian Osband, John Aslanides, and Albin Cassirer.  Randomized prior functions for deep reinforcement learning.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31, pages 8617–8629. Curran Associates, Inc., 2018.
  * [32] Arghya Bhowmik, Ivano E. Castelli, Juan Maria Garcia-Lastra, Peter Bjørn Jørgensen, Ole Winther, and Tejs Vegge.  A perspective on inverse design of battery interphases using multi-scale modelling, experiments and generative deep learning.  Energy Storage Materials, 21:446–456, 2019.

## Supplementary material:

Calibrated Uncertainty for Molecular Property Prediction
using Ensembles of Message Passing Neural Networks

##  Appendix A InChi comparison details

The overlapping set of structures that appear in both the QM9 [12] and PC9
[13] datasets were identified by comparing International Chemical Identifiers
(InChI) strings [29]. The InChI strings for both datasets were computed using
the Open Babel command line tool (obabel v. 3.1.0):

    $ obabel [input_file.xyz] -o inchi -xr -O [output_file.inchi]

or similarly for multiple files:

    $ for f in *.xyz;
    > do obabel $f -o inchi -xr -O ../inchi/${f:0:-3}inchi;
    > done

Then the InChi strings were truncated as to not differentiate between
stereoisomers (structures with the same chemical formula and connectivity).
Specifically, the `b }, {\verb `t , `m }, and {\verb `s layers of the InChi
strings were removed. When comparing the truncated InChi strings of the two
datasets, we found that that 21,777 molecules from QM9 are also in PC9 and
21,619 molecules from PC9 are also in QM9. The numbers are not identical since
QM9 and PC9 contains a different amount of duplicate truncated InChi strings,
so a structure from one dataset can appear multiple times in the other
dataset.

In [13] it was reported that 18,357 structures from PC9 also belong to QM9,
based on comparing InChi strings computed with the Open Babel software. We
were not able to reproduce this number using any combination of InChi layers,
so we instead used the method and result described above in this section.

##  Appendix B Additional results

(a)

(b)

Figure B1: Evaluation of uncertainty on the QM9 test set when predicting $E$.
Reliability diagram (a) and calibration plot (b).

(a)

(b)

Figure B2: Evaluation of uncertainty on the PC9 test set when predicting $E$.
Reliability diagram (a) and calibration plot (b).

(a)

(b)

Figure B3: Evaluation of uncertainty when training on QM9 and testing on PC9.
Reliability diagram (a) and calibration plot (b).

(a)

(b)

Figure B4: Evaluation of uncertainty when training on PC9 and testing on QM9.
Reliability diagram (a) and calibration plot (b).