#  _Ab initio_ computation of the longitudinal response function in 40Ca

J. E. Sobczyk  Institut für Kernphysik and PRISMA+ Cluster of Excellence,
Johannes Gutenberg-Universität, 55128 Mainz, Germany  B. Acharya  Institut für
Kernphysik and PRISMA+ Cluster of Excellence, Johannes Gutenberg-Universität,
55128 Mainz, Germany  S. Bacca  Institut für Kernphysik and PRISMA+ Cluster of
Excellence, Johannes Gutenberg-Universität, 55128 Mainz, Germany  Helmholtz-
Institut Mainz, Johannes Gutenberg-Universität Mainz, D-55099 Mainz, Germany
G. Hagen  Physics Division, Oak Ridge National Laboratory, Oak Ridge, TN
37831, USA  Department of Physics and Astronomy, University of Tennessee,
Knoxville, TN 37996, USA

###### Abstract

We present a consistent _ab initio_ computation of the longitudinal response
function $R_{L}$ in 40Ca using the coupled-cluster and Lorentz integral
transform methods starting from chiral nucleon-nucleon and three-nucleon
interactions. We validate our approach by comparing our results for $R_{L}$ in
4He and the Coulomb sum rule in 40Ca against experimental data and other
calculations. For $R_{L}$ in 40Ca we obtain a very good agreement with
experiment in the quasi-elastic peak up to intermediate momentum transfers,
and we find that final state interactions are essential for an accurate
description of the data. This work presents a milestone towards _ab initio_
computations of neutrino-nucleus cross sections relevant for experimental
long-baseline neutrino programs.

Understanding a wide variety of nuclear phenomena in terms of constituent
nucleons is a major ongoing initiative in nuclear theory Hergert (2020).
Theoretical predictions that start from the forces among nucleons and their
interactions with external probes as described by chiral effective field
theory are arguably the doorway to connect experimental observations with the
underlying fundamental theory of quantum chromo-dynamics van Kolck (1994);
Bedaque and van Kolck (2002); Epelbaum  _et al._ (2009); Machleidt and Entem
(2011). This approach is key to interpret existing data, provide guidance for
future experiments, and support interdisciplinary efforts at the interface
with nuclear physics, such as neutrino physics Alvarez-Ruso  _et al._ (2018).

Neutrino oscillation experiments aim at addressing some of the biggest
unanswered questions in physics by measuring the charge conjugation-parity
violating phase in the lepton sector of the Standard Model of particle
physics. For the current neutrino oscillation experiments the systematic
errors are at the order of $\sim 10\%$ and largely influenced by considerable
cross-section uncertainties. Next generation experiments set their precision
goal much higher. The T2HK Abe  _et al._ (2015) and DUNE Acciarri  _et al._
(2015) experiments aim at achieving much smaller statistical fluctuations,
comparable with present systematic errors. It is therefore crucial to control
the systematics, whose major part comes from the limited precision of
theoretical modeling of neutrino-nucleus cross sections. Furthermore, the
exposure needed to achieve a desired sensitivity also depends on the ability
of reducing systematic errors. The models which are presently in use,
particularly the ones implemented in the Monte Carlo event generators, should
be benchmarked with the predictions given by _ab initio_ models of nuclear
dynamics for relevant nuclei such as 12C, 16O and 40Ar.

Due to recent developments of accurate quantum many-body methods with
controlled approximations, ever-increasing computing power, and advancements
in the description of nuclear interactions and electroweak currents, we are
now entering an era where the _ab initio_ description of lepton-nucleus
scattering is becoming possible. The Green’s Function Monte Carlo (GFMC)
method was used to calculate nuclear responses of 4He and 12C Lovato  _et al._
(2014, 2015), and was recently able to make direct comparisons with the
neutrino-nucleus experimental cross sections Lovato  _et al._ (2018, 2020).
Using the same dynamical ingredients, other simplified methods are being
developed to reduce the computational load and address the quasi-elastic peak
Pastore  _et al._ (2020). In another set of studies, the lepton-nucleus
scattering cross sections of 4He, 16O and 40Ar were obtained using spectral
functions calculated in the self-consistent Green’s function method with
final-state interactions included using mean-field potentials Rocco and
Barbieri (2018); Barbieri  _et al._ (2019).

In this Letter, we lay out the tools for an _ab initio_ method that accurately
accounts for final state interactions, consistently with the treatment of
initial state interactions, and demonstrate its advantages by comparing to
available longitudinal electron scattering data for 40Ca. We base our approach
on the coupled–cluster (CC) method Coester (1958); Coester and Kümmel (1960);
Kümmel  _et al._ (1978); Mihaila and Heisenberg (2000); Dean and Hjorth-Jensen
(2004); Włoch  _et al._ (2005); Hagen  _et al._ (2008, 2010); Binder  _et al._
(2014); Hagen  _et al._ (2014), which stands out as one of the most suitable
and promising methods for calculations involving medium-mass and heavy nuclei
due to the polynomial scaling of its computational cost with the mass number
$A$. Initially applied to closed-shell nuclei (see Ref. Hagen  _et al._ (2014)
for a review), it has since been extended to accurately describe doubly open-
shell neighbors such as 40Ar Liu  _et al._ (2019); Payne  _et al._ (2019), and
more recently starting from an axially deformed reference state entire isotope
chains Novario  _et al._ (2020); Koszorús  _et al._ (2021). Combining CC with
the Lorentz integral transform method Efros  _et al._ (1994, 2007), the LIT-CC
approach extends the reach of this theory to processes involving excitation of
bound nuclear states to the continuum. Originally applied to calculate low-
energy nuclear dipole responses Bacca  _et al._ (2013, 2014), recently it was
extended to compute the Coulomb sum rule for 4He and 16O Sobczyk  _et al._
(2020). By devising a method to project out the spurious center-of-mass (CoM)
excitations, Ref. Sobczyk  _et al._ (2020) has also tackled the major
technical challenge of removing CoM contaminations in calculations utilizing
translationally non-invariant nuclear electroweak operators. These
developments open the door to go beyond the sum rule calculations and gain
deeper insights into the dynamics of the nucleus by computing the nuclear
response functions. With the goal of eventually applying the theory to
neutrino-nucleus scattering, where experimental data are scarce or imprecise,
we first benchmark our results for inelastic electron-scattering by comparing
them with existing data for 40Ca.

The inclusive cross section of this process can be expressed in terms of two
response functions: the longitudinal, $R_{L}(\omega,q)$, and the transverse,
$R_{T}(\omega,q)$, where $\omega$ is the energy transferred from the electron
to the nucleus. These are induced by the charge and the current operator,
respectively, and can be experimentally disentangled using the so-called
Rosenbluth separation. We study the longitudinal response in this work and
defer the transverse response, which receives large two-nucleon
electromagnetic current contributions Lovato  _et al._ (2020), to a future
work. Formally, the longitudinal response function can be defined as

$R_{L}(\omega,q)\\!=\\!\int\\!\\!\\!\\!\\!\\!\\!\sum_{f}|\\!\left\langle\Psi_{f}|\rho(q)|\Psi_{0}\right\rangle\\!|^{2}\,\delta\\!\\!\left(\\!E_{f}+\frac{q^{2}}{2M}-E_{0}-\omega\\!\right),$
(1)

where $M$ is the mass of the target nucleus, and $|\Psi_{0/f}\rangle$ and
$E_{0/f}$ respectively denote the initial/final-state nuclear wave functions
and energies, which we compute using nucleon-nucleon and three-nucleon forces
from chiral effective field theory. In order to estimate the sensitivity of
our results on the employed Hamiltonian we use two different chiral
interactions, namely NNLOsat Ekström  _et al._ (2015) and $\Delta$NNLOGO(450)
Jiang  _et al._ (2020). These interactions are both given at next-to-next-to-
leading order in the chiral expansion and employ a regulator cutoff of $450$
$\mathrm{MeV}/c$, but they differ in that $\Delta$NNLOGO(450) includes
explicit $\Delta$-isobars in its construction while NNLOsat does not. These
interactions are well suited for our study of 4He and 40Ca as they have been
shown to provide an accurate description of radii and binding energies of
light and medium-mass nuclei nuclei, and the saturation point of symmetric
nuclear matter Ekström  _et al._ (2015); Jiang  _et al._ (2020).

The charge density operator considered in this work is

${\rho}(q)=\frac{e}{2}\sum_{i=1}^{A}\,\left(G_{E}^{S}(Q^{2})+\tau_{i}^{3}\,G_{E}^{V}(Q^{2})\right)\exp{(i{\bf
q}\cdot{\bf r}_{i})}\,,$ (2)

where $e$ is the proton charge, while ${\bf r}_{i}$ and $\tau_{i}^{3}$ are the
coordinate and the third isospin component of nucleon $i$. We use the
parametrization of Ref. Kelly (2004) for the nucleon isoscalar/isovector
electric form factors, $G_{E}^{S/V}(Q^{2})$. The Darwin-Foldy and the spin-
orbit relativistic corrections, as well as the two-nucleon current
contributions, are not included in Eq. (2) since we strive for consistency
between the power-counting and truncation in the chiral expansions of the
current and the interactions. Specifically, corrections to Eq. (2) are at
least four orders higher in the chiral expansion when the inverse of the
nucleon mass is counted as two chiral orders Krebs  _et al._ (2019), which is
beyond the order at which the interactions we use are truncated.

The sum over $\Psi_{f}$ in Eq. (1) poses a serious computational challenge,
since it involves an integration over the continuum states, when $\omega$ is
above the particle emission threshold $\omega_{th}$. To overcome this issue,
we use the LIT method, where through the application of a Lorentzian-kernel
transform

${\cal L}_{L}(\sigma,q)=\frac{\sigma_{I}}{\pi}\int
d\omega\frac{R_{L}(\omega,q)}{(\omega-\sigma_{R})^{2}+\sigma_{I}^{2}}=\langle\widetilde{\Psi}_{\sigma,q}^{\rho}|\widetilde{\Psi}_{\sigma,q}^{\rho}\rangle$
(3)

with $\sigma_{I}\neq 0$, one reduces the problem to solving

$(H-E_{0}-\sigma)|\widetilde{\Psi}_{\sigma,q}^{\rho}\rangle={\rho}(q)|{\Psi_{0}}\rangle\,,$
(4)

where $H$ denotes the nuclear Hamiltonian. Effectively,
$\widetilde{\Psi}_{\sigma,q}^{\rho}$ is the solution of a bound-state
“Schrödinger-like” equation with a source term, which can be solved also in
coupled-cluster theory.

The CC method allows for the inclusion of many-body correlations as a
controlled expansion by writing the nuclear wave function as
$|\Psi\rangle=e^{T}|\Phi_{0}\rangle$. Here $|\Phi_{0}\rangle$ is a suitably
chosen reference state, and $T=T_{1}+T_{2}+\ldots$ is a linear expansion in
particle-hole excitations typically truncated at some low excitation rank. In
this work we truncate $T=T_{1}+T_{2}$ which is known as the coupled-cluster
singles and doubles (CCSD) method. Inserting the CCSD wave function into the
many-body Schrödinger equation and projecting from the left with $e^{-T}$, it
is seen that the reference state $|\Phi_{0}\rangle$ is the ground-state of the
similarity transformed normal-ordered Hamiltonian
$\overline{H}_{N}=e^{-T}H_{N}e^{T}$. In the LIT-CC formulation one has to
employ the equation-of-motion coupled-cluster technique (EOM-CC) Stanton and
Bartlett (1993) with a source term (see r.h.s. of Eq. (4)) and the similarity
transformed normal-ordered operator $\overline{\Theta}_{N}\equiv
e^{-T}\Theta_{N}e^{T}$ Miorelli  _et al._ (2018). Here, $\Theta$ are the
rank-$J$ multipoles of the electromagnetic charge operator given by Eq. (2).
To obtain the LIT, we perform EOM-CC calculations for each multipole
$[\rho(q)]^{J}$, and perform the sum over all multipoles at the end (see also
Ref. Acharya and Bacca (2020)).

The response function $R_{L}(\omega,q)$ for a given value of $q$ is then
obtained by inverting the integral transform from Eq. (3). To perform the
inversions, which require the solution of an ill-posed problem, we perform the
expansion
$R_{L}(\omega)=\sum_{i}^{N}c_{i}\omega^{n_{0}}e^{-\frac{\omega}{\beta i}}$ and
seek for stable solutions by varying the non-linear parameter $\beta$ (as well
as $n_{0}$) in a certain range. The inversion procedure involves the
determination of the coefficients $c_{i}$ of the $N$ basis functions by a
least-squares fit Efros  _et al._ (2007). We impose $R(\omega)$ to be zero for
$\omega\leq\omega_{th}$, using the values we obtain for a given nuclear
Hamiltonian in the CCSD approximation. We estimate the uncertainty associated
with the inversion procedure by inverting LITs with three different values of
$\sigma_{I}=5,10$ and 20 MeV and by varying $N$ from 6 to 9.

In all our results we employ a model space consisting of 15 major oscillator
shells ($e_{\rm max}=2n+l=14$) with an additional cut on the matrix elements
of the three-nucleon force given by $e_{3{\rm
max}}=2n_{1}+l_{1}+2n_{2}+l_{2}+2n_{3}+l_{3}\leq 16$. We checked that we can
reach a satisfactory convergence of $\mathcal{L}_{L}$ in terms of the single-
particle model space size $e_{\rm max}$. The latter can be tested, e.g., by
studying the residual dependence on the underlying harmonic oscillator
frequency $\hslash\Omega$. In particular, for LITs with $\sigma=20$ MeV we
estimate the convergence in the quasi-elastic peak to be at the 2$\%$ level
for $q\leq 350$ MeV/c and of $4\%$ for $q\geq 400$ MeV/c, by varying
$\hslash\Omega$ in the range 18 to 22 MeV.

Figure 1: Longitudinal response function for 4He at $q=300$ MeV/c. HH results
taken from Ref. Barnea  _et al._ (2001), GFMC results from Ref. Rocco  _et
al._ (2018), and experimental data from Ref. Carlson  _et al._ (2002).

Benchmark on the ${}^{4}\mathrm{He}$ nucleus— We begin by presenting our
results for $R_{L}$ in the case of 4He at $q=300$ MeV/c. In Fig. 1, we show
calculations performed with the NNLOsat interaction in the CCSD scheme for an
underlying harmonic oscillator frequency of $\hslash\Omega=16$ MeV. Here the
small band reflects only the uncertainty associated with the LIT inversion.
For comparison, we also show calculations performed with the hyperspherical
harmonics method (HH) Bacca  _et al._ (2009) using the AV18+UIX potential and
Green Function Monte Carlo (GFMC) Rocco  _et al._ (2018) calculations that
used the AV18+IL7 potential. We obtain very good agreement with the
experimental data as well as with other theoretical calculations. This
comparison corroborates our method and further validates the protocol we
developed in Ref. Sobczyk  _et al._ (2020) to remove center of mass
contamination.

Figure 2: 40Ca results for Coulomb sum rule for N2LOsat and $\hslash\omega=22$
MeV compared with CVMC results of Ref. Lonardoni  _et al._ (2017) and
experimental data taken from Ref. Williamson  _et al._ (1997).

Figure 3: Longitudinal response of 40Ca for $q=300,\,350,\,400$ MeV/c for
NNLOsat and $\Delta$NNLOGO(450) potentials. For $q=200$ MeV/c the strength of
excited states was quenched by factor of 2 for better visibility. Experimental
data taken from Williamson  _et al._ (1997).

Benchmark on the ${}^{40}\mathrm{Ca}$ nucleus — Following the same steps as in
Ref. Sobczyk  _et al._ (2020), we calculate the Coulomb sum rule for 40Ca
using the NNLOsat interaction. We observe that the CoM contamination is
negligible for $q>200$ MeV/c, and is overall much smaller than in the
previously considered cases of 4He and 16O Sobczyk  _et al._ (2020). In Fig. 2
we compare it to the cluster variational Monte Carlo (CVMC) results from Ref.
Lonardoni  _et al._ (2017) which used the AV18+UIX potential and included
Darwin-Foldy and spin-orbit corrections. Results are compatible at low-$q$ due
to the larger uncertainty in the CVMC curve, and show the same increasing
trend for $q>100$ MeV/c with small differences. We have verified that the
difference at $q=500$ MeV/c is mainly due to relativistic effects which we
omitted in order to be consistent with the chiral order we work at. Most
importantly, both theoretical predictions are in agreement with experimental
data Williamson  _et al._ (1997) in the range between 300 and 375 MeV/c and
are higher than the data above $q=400$ MeV/c, likely because experimental data
are obtained by integrating $R_{L}$ up to a finite $\omega$, and not up to
infinity as is done in the theoretical calculations. We consider this a
successful benchmark of our method and point out that only a mild Hamiltonian
dependence is observed.

The ${}^{40}\mathrm{Ca}$ longitudinal response function — We now turn to our
_ab initio_ calculation of $R_{L}$ in 40Ca where the full final state
interaction is considered. We choose 40Ca because we can compare our
calculations with existing data, and it is also a stepping stone for coupled-
cluster computations of neutrino scattering on 40Ar. For both NNLOsat and
$\Delta$NNLOGO(450) we perform computations of $R_{L}$ at the momentum
transfers $q=200$, 300, 350 and 400 MeV/c. In CCSD, the obtained ground-state
energies $E_{0}$ (proton separation energies $\omega_{th}$) are $300.1$ (6.32)
MeV and $322.12$ (6.12) MeV for the NNLOsat and the $\Delta$NNLOGO(450)
potential, respectively.

First, we find two bound excited $J^{\pi}=3^{-}$, $5^{-}$ states lying
respectively at 4.5(3.8) MeV and 4.7(4.0) MeV with the
NNLOsat($\Delta$NNLOGO(450)) interactions, which are in reasonable agreement
with experimental data at $3.7$ MeV $(J^{\pi}=3^{-})$ and at $4.5$ MeV
$(J^{\pi}=5^{-})$. We plot their strengths as a line in Fig. 3, and we observe
that it decreases with $q$. Second, for the continuum response we show a band
that reflects the uncertainty associated with the LIT inversion and the model
space, as we vary the harmonic oscillator frequency $\hslash\Omega$ from 18 to
20 and 22 MeV. As can be seen in Fig. 3, for each momentum transfer we observe
a mild dependence on the interaction, the latter being stronger at $q=200$
MeV/c. Comparing to the available experimental data from Ref. Williamson  _et
al._ (1997), we find a generally very good agreement, which is best for
$q=300$ MeV/c. At $q=400$ MeV/c, we see a quenching of the quasi-elastic peak
and an enhancement in the tail with respect to experiment. We speculate that
this could potentially be explained by relativistic boost effects Rocco  _et
al._ (2018) or by the fact that, especially at high $q$ and high $\omega$, we
are reaching the limits of applicability of chiral effective field theory set
by the regulator cutoff 450 MeV/c.

Finally, to quantify the effect of the final state interaction, we will
contrast the LIT-CC results with those of the simple plane wave impulse
approximation (PWIA). The point-proton longitudinal response function is
obtained in PWIA assuming one outgoing free proton with mass $m$ and a
spectator (A-1)-system with mass $M_{s}$,

$R_{L}^{\mathrm{PWIA}}(\omega,q)=\\!\\!\int\\!\\!d{\bf p}\,n({\bf
p})\,\delta\left(\omega-\frac{(\\!{\bf p}+{\bf q})^{2}}{2m}-\frac{\bf
p^{2}}{2M_{s}}-\omega_{th}\\!\right),$ (5)

and then augmented with nucleon electric form factors. Here $n({\bf p})$
represents the proton momentum distribution calculated from coupled-cluster
theory using the NNLOsat interaction, where CoM corrections are found to be
negligible com . Unlike the LIT-CC results, the PWIA curves shown in Fig. 3
are in poor agreement with the data: $(i)$ they miss the quasi-elastic peak
position by up to 20 MeV, $(ii)$ they overestimate considerably the
quasi–elastic peak size by up to 40$\%$ and $(iii)$ and they do not fully
account for the asymmetric shape of the response. The differences between the
LIT-CC and the PWIA results are very strong at lower $\omega$, where we
observe that even for the highest momentum transfers here considered $q=400$
MeV/c, we describe the experimental data very well. This highlights the
importance of consistently including the final state interaction.

In order to provide a prediction for future measurements as opposed to a sole
postdiction of existing data, we have calculated also the $q=200$ MeV/c
kinematics, where no data exist yet. While this low-$q$ range may be less
important for neutrino physics, this is where we have the largest uncertainty
band (range of low-$q$ and low-$\omega$). New precise data could provide
important tests of the _ab initio_ nuclear structure theory. An experimental
program in this direction is presently under development in Mainz Dor .

Conclusions— We performed an _ab initio_ calculation of the longitudinal
response function of 40Ca and obtained very good agreement with existing data.
Our results are a proof of principle that the LIT-CC method is suitable to
deliver responses for lepton-nucleus scattering at the momentum transfers
relevant for neutrino oscillation experiments. Consequently, we extended the
reach of consistent _ab initio_ calculations of electromagnetic responses at
intermediate momentum transfers into a region of medium-mass nuclei, which
until now was limited to systems with $A\leq 12$.

Our framework allows for quantification of uncertainties stemming from
truncations of model space, chiral effective-field-theory, and coupled-cluster
expansions. In this work, we estimated errors that arise from the inversion
procedure, and studied the dependencies on the model space and the nuclear
Hamiltonian. Our quantified uncertainties does not yet include effects of
missing higher-order excitations in the coupled-cluster expansion or terms in
the chiral effective field theory interactions and currents. A thorough
analysis of all theory uncertainties entering lepton-nucleus cross sections is
part of our future plans.

###### Acknowledgements.

We thank Nir Barnea and Thomas Papenbrock for useful comments and discussions.
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through
the Collaborative Research Center [The Low-Energy Frontier of the Standard
Model (SFB 1044)], and through the Cluster of Excellence “Precision Physics,
Fundamental Interactions, and Structure of Matter” (PRISMA+ EXC 2118/1) funded
by the DFG within the German Excellence Strategy (Project ID 39083149), by the
Office of Nuclear Physics, U.S. Department of Energy, under grants desc0018223
(NUCLEI SciDAC-4 collaboration) and by the Field Work Proposal ERKBP72 at Oak
Ridge National Laboratory (ORNL). Computer time was provided by the Innovative
and Novel Computational Impact on Theory and Experiment (INCITE) program and
by the supercomputer Mogon at Johannes Gutenberg-Universität Mainz. This
research used resources of the Oak Ridge Leadership Computing Facility located
at ORNL, which is supported by the Office of Science of the Department of
Energy under Contract No. DE-AC05-00OR22725.

## References

  * Hergert (2020) Heiko Hergert, “A guided tour of ab initio nuclear many-body theory,” Frontiers in Physics 8, 379 (2020), and references therein.
  * van Kolck (1994) U. van Kolck, “Few-nucleon forces from chiral Lagrangians,” Phys. Rev. C 49, 2932–2941 (1994).
  * Bedaque and van Kolck (2002) P. F. Bedaque and U. van Kolck, “Effective field theory for few-nucleon systems,” Annual Review of Nuclear and Particle Science 52, 339–396 (2002), nucl-th/0203055 .
  * Epelbaum  _et al._ (2009) E. Epelbaum, H.-W. Hammer, and Ulf-G. Meißner, “Modern theory of nuclear forces,” Rev. Mod. Phys. 81, 1773–1825 (2009).
  * Machleidt and Entem (2011) R. Machleidt and D.R. Entem, “Chiral effective field theory and nuclear forces,” Physics Reports 503, 1 – 75 (2011).
  * Alvarez-Ruso  _et al._ (2018) L. Alvarez-Ruso _et al._ (NuSTEC), “NuSTEC White Paper: Status and challenges of neutrino–nucleus scattering,” Prog. Part. Nucl. Phys. 100, 1–68 (2018), arXiv:1706.03621 [hep-ph] .
  * Abe  _et al._ (2015) K. Abe _et al._ (Hyper-Kamiokande Proto-Collaboration), “Physics potential of a long-baseline neutrino oscillation experiment using a J-PARC neutrino beam and Hyper-Kamiokande,” PTEP 2015, 053C02 (2015).
  * Acciarri  _et al._ (2015) R. Acciarri _et al._ (DUNE), “Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE),” (2015), arXiv:1512.06148 [physics.ins-det] .
  * Lovato  _et al._ (2014) A. Lovato, S. Gandolfi, J. Carlson, Steven C. Pieper, and R. Schiavilla, “Neutral weak current two-body contributions in inclusive scattering from ${}^{12}\mathrm{C}$,” Phys. Rev. Lett. 112, 182502 (2014).
  * Lovato  _et al._ (2015) A. Lovato, S. Gandolfi, J. Carlson, Steven C. Pieper, and R. Schiavilla, “Electromagnetic and neutral-weak response functions of 4He and 12C,” Phys. Rev. C 91, 062501 (2015), arXiv:1501.01981 [nucl-th] .
  * Lovato  _et al._ (2018) A. Lovato, S. Gandolfi, J. Carlson, Ewing Lusk, Steven C. Pieper, and R. Schiavilla, “Quantum Monte Carlo calculation of neutral-current $\nu-^{12}C$ inclusive quasielastic scattering,” Phys. Rev. C 97, 022502 (2018), arXiv:1711.02047 [nucl-th] .
  * Lovato  _et al._ (2020) A. Lovato, J. Carlson, S. Gandolfi, N. Rocco, and R. Schiavilla, “Ab initio study of $({\nu}_{\ell},{\ell}^{-})$ and $({\overline{\nu}}_{\ell},{\ell}^{+})$ inclusive scattering in ${}^{12}\mathrm{C}$: Confronting the miniboone and t2k ccqe data,” Phys. Rev. X 10, 031068 (2020).
  * Pastore  _et al._ (2020) S. Pastore, J. Carlson, S. Gandolfi, R. Schiavilla, and R. B. Wiringa, “Quasielastic lepton scattering and back-to-back nucleons in the short-time approximation,” Phys. Rev. C 101, 044612 (2020).
  * Rocco and Barbieri (2018) N. Rocco and C. Barbieri, “Inclusive electron-nucleus cross section within the Self Consistent Green’s Function approach,” Phys. Rev. C 98, 025501 (2018), arXiv:1803.00825 [nucl-th] .
  * Barbieri  _et al._ (2019) C. Barbieri, N. Rocco, and V. Somà, “Lepton scattering from 40Ar and Ti in the quasielastic peak region,” (2019), arXiv:1907.01122 [nucl-th] .
  * Coester (1958) F. Coester, “Bound states of a many-particle system,” Nuclear Physics 7, 421 – 424 (1958).
  * Coester and Kümmel (1960) F. Coester and H. Kümmel, “Short-range correlations in nuclear wave functions,” Nuclear Physics 17, 477 – 485 (1960).
  * Kümmel  _et al._ (1978) H. Kümmel, K. H. Lührmann, and J. G. Zabolitzky, “Many-fermion theory in expS- (or coupled cluster) form,” Physics Reports 36, 1 – 63 (1978).
  * Mihaila and Heisenberg (2000) B. Mihaila and J. H. Heisenberg, “Microscopic Calculation of the Inclusive Electron Scattering Structure Function in ${}^{16}O$,” Phys. Rev. Lett. 84, 1403–1406 (2000).
  * Dean and Hjorth-Jensen (2004) D. J. Dean and M. Hjorth-Jensen, “Coupled-cluster approach to nuclear physics,” Phys. Rev. C 69, 054320 (2004).
  * Włoch  _et al._ (2005) M. Włoch, D. J. Dean, J. R. Gour, M. Hjorth-Jensen, K. Kowalski, T. Papenbrock, and P. Piecuch, “Ab-Initio coupled-cluster study of ${}^{16}\mathrm{O}$,” Phys. Rev. Lett. 94, 212501 (2005).
  * Hagen  _et al._ (2008) G. Hagen, T. Papenbrock, D. J. Dean, and M. Hjorth-Jensen, “Medium-mass nuclei from chiral nucleon-nucleon interactions,” Phys. Rev. Lett. 101, 092502 (2008).
  * Hagen  _et al._ (2010) G. Hagen, T. Papenbrock, D. J. Dean, and M. Hjorth-Jensen, “Ab initio coupled-cluster approach to nuclear structure with modern nucleon-nucleon interactions,” Phys. Rev. C 82, 034330 (2010).
  * Binder  _et al._ (2014) Sven Binder, Joachim Langhammer, Angelo Calci, and Robert Roth, “Ab initio path to heavy nuclei,” Phys. Lett. B 736, 119 – 123 (2014).
  * Hagen  _et al._ (2014) G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, “Coupled-cluster computations of atomic nuclei,” Rep. Prog. Phys. 77, 096302 (2014).
  * Liu  _et al._ (2019) H. N. Liu, A. Obertelli, P. Doornenbal, C. A. Bertulani, G. Hagen, J. D. Holt, G. R. Jansen, T. D. Morris, A. Schwenk, R. Stroberg, N. Achouri, H. Baba, F. Browne, D. Calvet, F. Château, S. Chen, N. Chiga, A. Corsi, M. L. Cortés, A. Delbart, J.-M. Gheller, A. Giganon, A. Gillibert, C. Hilaire, T. Isobe, T. Kobayashi, Y. Kubota, V. Lapoux, T. Motobayashi, I. Murray, H. Otsu, V. Panin, N. Paul, W. Rodriguez, H. Sakurai, M. Sasano, D. Steppenbeck, L. Stuhl, Y. L. Sun, Y. Togano, T. Uesaka, K. Wimmer, K. Yoneda, O. Aktas, T. Aumann, L. X. Chung, F. Flavigny, S. Franchoo, I. Gašparić, R.-B. Gerst, J. Gibelin, K. I. Hahn, D. Kim, T. Koiwai, Y. Kondo, P. Koseoglou, J. Lee, C. Lehr, B. D. Linh, T. Lokotko, M. MacCormick, K. Moschner, T. Nakamura, S. Y. Park, D. Rossi, E. Sahin, D. Sohler, P.-A. Söderström, S. Takeuchi, H. Törnqvist, V. Vaquero, V. Wagner, S. Wang, V. Werner, X. Xu, H. Yamada, D. Yan, Z. Yang, M. Yasuda, and L. Zanetti, “How robust is the $n=34$ subshell closure? first spectroscopy of ${}^{52}\mathrm{Ar}$,” Phys. Rev. Lett. 122, 072502 (2019).
  * Payne  _et al._ (2019) C. G. Payne, S. Bacca, G. Hagen, W. Jiang, and T. Papenbrock, “Coherent elastic neutrino-nucleus scattering on 40Ar from first principles,” Phys. Rev. C 100, 061304 (2019), arXiv:1908.09739 [nucl-th] .
  * Novario  _et al._ (2020) S. J. Novario, G. Hagen, G. R. Jansen, and T. Papenbrock, “Charge radii of exotic neon and magnesium isotopes,” Phys. Rev. C 102, 051303 (2020).
  * Koszorús  _et al._ (2021) Á. Koszorús, X. F. Yang, W. G. Jiang, S. J. Novario, S. W. Bai, J. Billowes, C. L. Binnersley, M. L. Bissell, T. E. Cocolios, B. S. Cooper, R. P. de Groote, A. Ekström, K. T. Flanagan, C. Forssén, S. Franchoo, R. F. Garcia Ruiz, F. P. Gustafsson, G. Hagen, G. R. Jansen, A. Kanellakopoulos, M. Kortelainen, W. Nazarewicz, G. Neyens, T. Papenbrock, P. G. Reinhard, C. M. Ricketts, B. K. Sahoo, A. R. Vernon, and S. G. Wilkins, “Charge radii of exotic potassium isotopes challenge nuclear theory and the magic character of n = 32,” Nature Physics (2021), 10.1038/s41567-020-01136-5.
  * Efros  _et al._ (1994) Victor D. Efros, Winfred Leidemann, and Giuseppina Orlandini, “Response functions from integral transforms with a lorentz kernel,” Phys. Lett. B 338, 130 – 133 (1994).
  * Efros  _et al._ (2007) V D Efros, W Leidemann, G Orlandini, and N Barnea, “The lorentz integral transform (lit) method and its applications to perturbation-induced reactions,” Journal of Physics G: Nuclear and Particle Physics 34, R459 (2007).
  * Bacca  _et al._ (2013) S. Bacca, N. Barnea, G. Hagen, G. Orlandini, and T. Papenbrock, “First principles description of the giant dipole resonance in ${}^{16}\mathbf{O}$,” Phys. Rev. Lett. 111, 122502 (2013).
  * Bacca  _et al._ (2014) S. Bacca, N. Barnea, G. Hagen, M. Miorelli, G. Orlandini, and T. Papenbrock, “Giant and pigmy dipole resonances in 4He, 16,22O, and 40Ca from chiral nucleon-nucleon interactions,” ArXiv e-prints (2014), arXiv:1410.2258 [nucl-th] .
  * Sobczyk  _et al._ (2020) J. E. Sobczyk, B. Acharya, S. Bacca, and G. Hagen, “Coulomb sum rule for 4He and 16O from coupled-cluster theory,” Phys. Rev. C 102, 064312 (2020), arXiv:2009.01761 [nucl-th] .
  * Ekström  _et al._ (2015) A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, “Accurate nuclear radii and binding energies from a chiral interaction,” Phys. Rev. C 91, 051301 (2015), arXiv:1502.04682 [nucl-th] .
  * Jiang  _et al._ (2020) W. G. Jiang, A. Ekström, C. Forssén, G. Hagen, G. R. Jansen, and T. Papenbrock, “Accurate bulk properties of nuclei from $a=2$ to $\infty$ from potentials with $\mathrm{\Delta}$ isobars,” Phys. Rev. C 102, 054301 (2020).
  * Kelly (2004) J.J. Kelly, “Simple parametrization of nucleon form factors,” Phys. Rev. C 70, 068202 (2004).
  * Krebs  _et al._ (2019) H. Krebs, E. Epelbaum, and U. G. Meißner, “Nuclear Electromagnetic Currents to Fourth Order in Chiral Effective Field Theory,” Few Body Syst. 60, 31 (2019), arXiv:1902.06839 [nucl-th] .
  * Stanton and Bartlett (1993) John F. Stanton and Rodney J. Bartlett, “The equation of motion coupled-cluster method. a systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties,” J. Chem. Phys. 98, 7029–7039 (1993).
  * Miorelli  _et al._ (2018) M. Miorelli, S. Bacca, G. Hagen, and T. Papenbrock, “Computing the dipole polarizability of ${}^{48}\mathrm{Ca}$ with increased precision,” Phys. Rev. C 98, 014324 (2018).
  * Acharya and Bacca (2020) Bijaya Acharya and Sonia Bacca, “Neutrino-deuteron scattering: Uncertainty quantification and new ${L}_{1,A}$ constraints,” Phys. Rev. C 101, 015505 (2020).
  * Barnea  _et al._ (2001) Nir Barnea, Winfried Leidemann, and Giuseppina Orlandini, “State-dependent effective interaction for the hyperspherical formalism with noncentral forces,” Nuclear Physics A 693, 565 – 578 (2001).
  * Rocco  _et al._ (2018) Noemi Rocco, Winfried Leidemann, Alessandro Lovato, and Giuseppina Orlandini, “Relativistic effects in ab-initio electron-nucleus scattering,” Phys. Rev. C 97, 055501 (2018), arXiv:1801.07111 [nucl-th] .
  * Carlson  _et al._ (2002) J. Carlson, J. Jourdan, R. Schiavilla, and I. Sick, “Longitudinal and transverse quasielastic response functions of light nuclei,” Phys. Rev. C 65, 024002 (2002).
  * Bacca  _et al._ (2009) Sonia Bacca, Nir Barnea, Winfried Leidemann, and Giuseppina Orlandini, “Role of the final-state interaction and three-body force on the longitudinal response function of ${}^{4}\mathrm{He}$,” Phys. Rev. Lett. 102, 162501 (2009).
  * Lonardoni  _et al._ (2017) D. Lonardoni, A. Lovato, Steven C. Pieper, and R. B. Wiringa, “Variational calculation of the ground state of closed-shell nuclei up to $a=40$,” Phys. Rev. C 96, 024326 (2017).
  * Williamson  _et al._ (1997) C. F. Williamson _et al._ , “Quasielastic electron scattering from Ca-40,” Phys. Rev. C 56, 3152–3172 (1997).
  * (48) J. E. Sobczyk, B. Acharya, S. Bacca, G. Hagen, and T. Papenbrock, in preparation.
  * (49) L. Doria and M. Mihovilovic, private communication.


# Geometry-aware Transformer for
molecular property prediction

Bumju Kwak
Recommendation Team
Kakao Corporation
Gyeonggi, Republic of Korea
<EMAIL_ADDRESS>
&Jeonghee Jo
Institute of New Media and Communications
Seoul National University
Seoul 08826, Republic of Korea
<EMAIL_ADDRESS>
Byunghan Lee
Department of Electronic Engineering
Seoul National University of Science and Technology
Seoul 01811, Republic of Korea
<EMAIL_ADDRESS>
&Sungroh Yoon
Department of Electrical and Computer Engineering
Seoul National University
Seoul 08826, Republic of Korea
<EMAIL_ADDRESS>

###### Abstract

Recently, message passing neural networks (MPNNs) have achieved remarkable
performances for various chemical tasks. However, the meesage passing scheme
can only cover a local field, and cannot capture long-range relationship
between atoms. This behavior does not accord with chemical theory, which is
one of the fundamental limitations of the MPNNs. In this work, we propose a
novel Transformer-based framework for representation learning on molecules,
named Geometry-aware Transformer (GeoT). GeoT learns sequential
representations of molecule graphs based on global attentions. By virtue of
the Transformer-based architecture and additional refinements, the smoothing
problem can be significantly alleviated compared to graph convolution-based
architectures. GeoT achieved comparable performances with other recent MPNN-
based models. We also verified the intuition behind GeoT with respect to the
reasonability and interpretability, through attention visualization.

##  1 Introduction

Quantum mechanical (QM) calculations have been used for development of
chemicals in various fields, such as drugs, catalysts, or novel industrial
materials. Density functional theory (DFT) is the most widely used
computational modeling method used for modern QM research. However, DFT
requires tremendous computations to predict properties for a small molecule.

For this reason, several machine learning-based methods have been developed to
replace DFT calculation with more cost-effective alternatives [1, 2, 3, 4] .
Motivated by recent advances in deep learning, many recent models [5, 6, 7, 8,
9, 10, 11] adopted deep learning-based methods to predict various molecular
properties including energy and forces.

Message passing neural network (MPNN) [12] is one of the most widely used
frameworks for graph neural networks (GNNs) in recent days. In these kinds of
works, a molecule is described as a graph, and its atoms are represented as
discrete nodes scattered in Euclidean space. These models do not consider any
bond type information. Instead, they construct edges between all node pairs
located within the fixed cutoff distance from each other. Then, the message
passing layers do a graph convolution with each atom and its local
substructures (messages) which are controlled by radial basis filters. Many
atomistic models [6, 12, 7, 13] based on MPNN showed competitive performances
in various QM tasks.

Although promising performances on existing MPNNs, there are several
fundamental limitations in the localized view. In theory, every atom pair in a
molecule has their own interaction defined by their charges and distance. In
other words, atom-atom interactions exist over a broad range in the molecule.
Therefore, all atom-atom interactions should be considered properly in
molecular property prediction models regardless of their distances. However,
the localized convolutions of MPNN have conflict with this conformational
behavior of molecules, since MPNN with a finite cutoff value assumes that
messages are transferred within a restricted region. It means that MPNN cannot
identify any interatomic relationship between atom pairs located outside the
cutoff distance from each other. Consequently, MPNNs fail to capture long-
range interatomic relationships, which is not a desired property of molecular
representation.

Furthermore, the MPNN models are vulnerable to over-smoothing problem [14],
which is caused by inherent Laplacian-smoothing effects of message passing.
The atom features cannot be distinguished from each other, because the
message-passing updates an atom feature by mixing its neighbors in every step.
Recent studies [14, 15, 16] reported that the over-smoothing restricts the
performance of deep GNNs in various tasks including molecular dataset.

To address the above issues, we developed a novel geometry-aware Transformer
(GeoT) consisting of softmax-free attentions, for reasonable structure
representation combined with a global view. We removed the positional encoding
and introduced radial basis functions (RBFs) to reflect a complete molecular
structure and ensuring the permutation invariance of a molecular sequence. We
adopted several modifications and refinements to the original self-attention,
for enhancing the performance and preventing from the over-smoothing problem.

We summarize our contributions as below.

  * •

We propose a Transformer-based architecture GeoTransformer for molecular
graphs, incorporating radial basis embeddings of molecular conformations.

  * •

GeoT is simple yet effective, by capturing global three dimensional structure
with less over-smoothing problem. In addition, GeoT provides an intuitive
interpretation of atom-atom interactions via the attention values with a
global view.

  * •

GeoT achieved competitive performances in various types of benchmarks.
Furthermore, we adopted several model refinement tricks to GeoT and analyzed
effects thereof.

##  2 Related Works

###  2.1 Graph neural networks for molecules

The molecular representation as a graph has been widely used from early
stages. These methods regard a molecule as a graph consisting of discrete atom
nodes with bonds as edges. BPNN [3] is the first neural network model which
calculates energy from molecular graphs. The first work using graph
convolution on molecules [5] trained each atom representation with localized
filters. GDML [17, 18, 19] used conservative force fields of each atom, and
MPNN [12] proposed message passing framework for local update of each atom
representation in the network.

SchNet [6], PhysNet [7] and several MPNN-based models constructed localized
messages centered on atoms with continuous filters based on interatomic
distances. More recent works such as DimeNet [10], DimeNet++ [13], GemNet
[20], and SphereNet [21] compute angles to describe molecular conformations
more concretely. Some works like TFN [22], Cormorant [8], and
SE(3)-Transformer [23] adopt group equivariances in their models to preserve
roto-translational symmetry. These previous models need a cutoff value in
common, which restricts the receptive field of an atom-wise convolution. In
addition, these previous methods depends on tremendous computations for
geometric embedding for embedding explicit equivariances.

###  2.2 Self-attention and Transformer on molecular graphs

The self-attention has been widely adopted for learning sequential data,
because it can capture the relations between two components from the sequence.
Transformer [24] is a type of encoder-decoder architecture with a stack of
self-attention layers, to efficiently handle sequential data such as sentences
regardless of their lengths. We also applied AttnScale [25] and other
strategies for effective implementation of the attention mechanism on graphs.

For molecular property prediction tasks, Transformer-like structures have been
used for sequential representation of molecules such as SMILES [26], rather
than graph-structured form. MAT [27] and GRAT [9] used a self-attention on
molecular graphs, whereas the distance representations between atom pairs of
their model were incomplete to capture long-range interactions. Recently, [28,
29, 30], which are another Transformer-based models with a global view of
attention. However, they only count hops between two nodes rather than real
geometric distances.

##  3 Preliminary

###  3.1 Notation

We denote a molecule as a set of $N$ atoms
$\\{a_{i}\\}_{i\in\\{1,\ldots,N\\}}$ with coordinates
$\\{\mathbf{r_{i}}\\}\in\mathbb{R}^{3}$. We calculate the Euclidean distance
$r_{ij}=\lVert\mathbf{r_{i}}-\mathbf{r_{j}}\rVert_{2}$ between two atom
$a_{i}$ and $a_{j}$. Our model is a stack of $L$ layers, which transforms
input $\bm{X}^{l}$ into output $\bm{X}^{l+1}$, especially the first layer
input $\bm{X}^{0}$ is given from embedding layer which only depends on atomic
number.

###  3.2 Transformer encoder

Figure 1: Overall architecture of Geometry-aware Transformer

####  3.2.1 Multi-head self attention (MSA)

For given set of vectors $\\{\bm{X}^{l}_{i}\\}_{i\in\\{1,...,N\\}}$, MSA is
calculated by incorporating vectors using weight called self-attention, which
represents interaction between $\bm{X}^{l}_{i}$ and $\bm{X}^{l}_{j}$. Self-
attention first transforms input feature $\bm{X}$ into three types of tensors,
query $\bm{Q}$, key $\bm{K}$, and value $\bm{V}$, respectively. Each matrices
are obtained by multiplication of trainable weights $W_{Q},W_{K},W_{V}$.

$\bm{Q}=W_{Q}\bm{X},\;\;\bm{K}=W_{K}\bm{X},\;\;\bm{V}=W_{V}\bm{X}$ (1)

With $\bm{Q},\bm{K},and\bm{V}$, self-attention yields weighted sum of
$\bm{V}$, which weights are given by scaled dot-product between $\bm{Q}$ and
$\bm{K}$ scaled by vector dimension $d_{m}$ .

$\mathrm{Attention}(\bm{Q},\bm{K},\bm{V})=\mathrm{softmax}(\bm{QK}^{T}/\sqrt{d_{m}})\cdot\bm{V}$
(2)

MSA [24] was proposed to improve model performance instead of using single
self-attention only. For MSA, $\\{\bm{Q},\bm{K},\bm{V}\\}$ are split into $h$
matrices $\\{\bm{Q}_{i},\bm{K}_{i},\bm{V}_{i}\\}_{i\in\\{1,...,h\\}}$ with
dimension $d_{m}/h$. Self-attention is then applied to
$\\{\bm{Q}_{i},\bm{K}_{i},\bm{V}_{i}\\}_{i\in\\{1,...,h\\}}$, followed by
concatenation.

$\begin{split}\bm{H}_{i}&=\mathrm{Attention}(\bm{Q}_{i},\bm{K}_{i},\bm{V}_{i})\;\;\;\;i\in\\{1,...,h\\}\\\
\mathrm{MSA}(\bm{X})&=\mathrm{Concat}(\bm{H}_{1},...,\bm{H}_{h})\\\
\end{split}$ (3)

####  3.2.2 Transformer encoder layer

In the Transformer encoder [24] layer, the output from
$\mathrm{MSA}(\bm{X}^{l})$ added to original input $\mathbf{X}^{l}$ by
residual connection accompanied by layer normalization [31]. After that, the
last output $\tilde{\mathbf{X}}^{l}$ fed into feed forward network
($\mathrm{FFN}$), consisting of two linear transformations with an ELU [32]
activation. Finally, another layer normalization is applied to the output from
FFN added by $\tilde{\mathbf{X}}^{l}$.

$\begin{split}\tilde{\mathbf{X}}^{l}&=\mathrm{LayerNorm}(\mathrm{MSA}(\bm{X}^{l})+\bm{X}^{l})\\\
\bm{X}^{l+1}&=\mathrm{LayerNorm}(\mathrm{FFN}(\tilde{\mathbf{X}}^{l})+\tilde{\mathbf{X}}^{l})\end{split}$
(4)

##  4 Geometry-aware Transformer

We propose GeoT architecture to solve chemical property prediction task with
capturing global perspective of a molecule structure. GeoT is built based on
consecutive Transformer encoders, with modification for capturing geometric
information in Euclidean space.

Geometry-aware Transformers create distance embeddings expanded by multiple
RBFs. The embeddings from RBFs are merged with query $\mathbf{Q}$ and key
$\mathbf{K}$, which yields geometry-aware attention. GeoT is achieved by
calculating weighted sum of the values $\bm{V}$ with weights given from
geometry-aware attention. At the end of layer, the readout layer produces
predicted graph labels.

###  4.1 Two-body kernel

GeoT adopt distance information by two-body kernel
$\mathbf{\Lambda}:[1,...,N]\times[1,...,N]\rightarrow\mathbb{R}$ with its
matrix representation $\mathbf{\Lambda}$. To construct $\mathbf{\Lambda}$, we
define multiple RBFs $\mathbf{g}(r):\mathbb{R}\rightarrow\mathbb{R}^{d_{m}}$
with distances between atom pairs. Many previous studies suggested different
types of radial basis functions for these distances [6, 7, 13, 8]. We selected
multiple different Gaussian basis functions for RBFs, proposed in SchNet [6]
followed by concatenation, where $\gamma=10$ adjusts the scales of distance
representations and $\delta=0.1\textup{~{}\AA}$ is multiplier.

$\mathbf{g}(r)_{k}=exp(-\gamma||r-\delta k||^{2})\;\;$ (5)

With $\mathbf{g}(r)_{k}$, we created two-body kernel ${\Lambda}^{n}$ depending
distance $r$ with two-layer multilayer perceptron (MLP) $f^{n}$ as follows.

${\Lambda}^{n}(r)=f^{n}(\mathbf{g}(r))\\\ $ (6)

Since atomic interaction not only depends on distance between two atoms but
their atomic numbers, it is natural to combine atomic number information with
distance information $\mathbf{g}(r)$ to describe two-body interaction between
atoms. For this reason, we added vector representation $Z^{n}$ of two atom
$a_{i}$ and $a_{j}$ symmetrically to $\mathbf{g}(r)$. The two-body kernels are
created and trained in each layer individually.

$\begin{split}\mathbf{g}_{\mathrm{emb}}(r_{ij},a_{i},a_{j})&=\mathbf{g}(r_{ij})\oplus(Z^{n}(a_{i})+Z^{n}(a_{j}))),Z^{n}:\mathbb{N}\rightarrow\mathbb{R}^{d_{\mathrm{emb},2}}\\\
\bm{\Lambda}^{n}_{ij}&=f^{n}_{\mathrm{emb}}(\mathbf{g}_{\mathrm{emb}}(r_{ij},a_{i},a_{j}))\\\
\end{split}$ (7)

###  4.2 Geometry-aware self-attention

Note that we did not use the standard self-attention in our model. Our
Geometry-aware self-attention is different from the standard self-attention in
two points: 1) inserting two-body kernel matrix $\bm{\Lambda}$ and 2) the
absence of softmax operation. In details, we replaced the original attention
$Attn(\bm{Q},\bm{K})=(\bm{QK}^{T})/\sqrt{d_{m}}$ in Eq. 3 with geometry-aware
attention
$Attn_{G}(\bm{Q},\bm{K},\bm{\Lambda})=(\bm{QK}^{T}\bm{\Lambda})/\sqrt{d_{m}}$,
by introducing RBF from Eq. 7 to capture interactions between atoms.

Figure 2: Comparison between the original Transformer and Geometry-aware
Transformer. The key differences of two architectures are that the Geometry-
aware Transformer 1) adopted the radial basis functions to represent distances
of atom pairs, and 2) removed softmax function existed in the original
Transformer.

$\begin{split}\mathrm{GeoAttention}(\bm{Q},\bm{K},\bm{V},\bm{\Lambda})=(\bm{QK}^{T}\bm{\Lambda})/\sqrt{d_{m}}\cdot\bm{V}\end{split}$
(8)

The corresponding MSA is defined as follows.

$\begin{split}\bm{H}^{G}_{i}&=\mathrm{GeoAttention}(\bm{Q}_{i},\bm{K}_{i},\bm{V}_{i},\bm{\Lambda})\;\;\;\;i\in\\{1,...,h\\}\\\
\mathrm{MSA}_{G}(\bm{X})&=\mathrm{Concat}(\bm{H}^{G}_{1},...,\bm{H}^{G}_{h})\\\
\end{split}$ (9)

###  4.3 Readout layer

We introduced a pooling layer at the final step, to obtain molecular
properties by aggregating atom features. We applied sum-pooling and a linear
layer to get the final scalar-valued predictions.

$\hat{y}=W_{\mathrm{pool}}\bm{X^{L}}+b,\;\;$ (10)

###  4.4 Model refinement

Despite Transformer is more powerful in oversmoothing phenomenon, such works
reported that the original Transformer architecture still suffers by over-
smoothing problem [25, 33], resulting deeper models fails to converge.
Therefore, we introduced several strategies to GeoT layers for more robust and
stable convergence.

Li et al. [14] found that a graph convolution is a kind of a Laplacian
smoothing on graphs. It means that a GCN plays a role in repeated smoothing
which may cause undesirable and harmful effects on graph tasks, especially
when the network depth is deeper. In recent days, many studies [14, 15, 16]
argued that this is problematic and must be overcome, for both enhanced
performance in various tasks and richer representation learning on graphs.
Several studies pointed out that the over-smoothing problem is not restricted
in GNN architectures, but also exists in variants of Transformer architectures
[25, 34, 33], because the self-attention is also a type of the global
aggregating operation.

####  4.4.1 AttnScale

For this reason, we adopted several strategies to minimize smoothing effect in
Transformer. First, we applied AttnScale [25], which is developed to reduce
over-smoothing effect in deep Transformer with a simple strategy.

Transformer networks with AttnScale accentuate more meaningful pattern by
amplifying high-frequency signal of self-attention weights.

In specifically, AttnScale layer aims to amplify high-frequency filter aspect
of attention map in order to diminish inherent low-frequency filter property
[25]. Since low-frequency component of attention map as
$\bm{A}^{n}_{DC}=\frac{1}{N}\bm{1}\bm{1}^{T}$, remnant considered as high-
frequency term of attention map, where $\bm{A}^{n}_{DC}$ means the direct
component (DC) of an attention. AttnScale layer which scale high-frequency
term by $1+w_{a}$ which is greater than one, achieves attention map preserving
more HF environment of input.

$\begin{split}\bm{A}&=\mathrm{Attn}(\bm{Q},\bm{K})\\\
\bm{A}_{DC}&=\frac{1}{N}\bm{1}\bm{1}^{T}\\\ \bm{A}_{HF}&=\bm{A}-\bm{A}_{DC}\\\
\mathrm{Attn}^{\prime}(\bm{Q},\bm{K})&=\bm{A}_{DC}+(1+w_{a})\bm{A}_{HF}\\\
\end{split}$ (11)

In this paper, we modified AttnScale to consider our unnormalized attention
map without softmax layer. Especially, we replace original low-frequency
component $\bm{A}^{n}_{DC}$ with $\bm{A}^{\prime
n}_{DC}=\frac{1}{N}\sum_{k}A^{n}_{ik}$ to preserve the equivalence of original
$\bm{A}^{n}_{DC}$ with channel-wise average of attention map.

$\begin{split}\bm{A}^{n}&=\mathrm{Attn}_{G}(\bm{Q},\bm{K},\bm{\Lambda})\\\
\bm{A}^{\prime n}_{DC}&=\frac{1}{N}\sum_{k}A^{n}_{ik}.\\\ \bm{A}^{\prime
n}_{HF}&=\bm{A}^{n}-\bm{A}^{\prime n}_{DC}\\\
\mathrm{Attn}^{\prime}_{G}(\bm{Q},\bm{K},\bm{\Lambda})&=\bm{A}^{\prime
n}_{DC}+(1+w^{n}_{a})\bm{A}^{\prime n}_{HF}\\\
\mathrm{GeoAttention}^{\prime}(\bm{Q},\bm{K},\bm{V},\bm{\Lambda})&=\mathrm{Attn}^{\prime}_{G}(\bm{Q},\bm{K},\bm{\Lambda})\cdot\bm{V}\\\
\end{split}$ (12)

####  4.4.2 Attention with Parallel MLP

Our second strategy modifies layer normalization scheme. We were inspired from
the results of [33], which analyze the nature of MSA. Following from the above
paper, we modified the network structure outside the MSA. Original Transformer
encoder implements two successive layer normalizations for MSA with residual
connection, and feed forward network (FFN) with residual connections,
respectively, as below.

$\begin{split}\widetilde{\bm{X}^{l}}&=\mathrm{LayerNorm}(\mathrm{MSA}_{G}(\bm{X}^{l})+\bm{X}^{l})\\\
\bm{X}^{l}&=\mathrm{LayerNorm}(\mathrm{FFN}(\widetilde{\bm{X}}^{l}_{i})+\widetilde{\bm{X}}^{l}_{i})\\\
\end{split}$ (13)

Alternatively, we modified the above equations by implementing only one layer
normalization which takes MSA, FFN, and residual connection simultaneously. It
can be formulated as below.

$\begin{split}\bm{X}^{l+1}_{i}=\mathrm{LayerNorm}(\mathrm{MSA}_{G}(\bm{X}^{l}_{i})+\mathrm{FFN}(\bm{X}^{l}_{i})+\bm{X}^{l}_{i})\\\
\end{split}$ (14)

##  5 Experiments

###  5.1 Dataset

We evaluated our model on two benchmark datasets QM9 [35, 36] MD17 [17, 18,
19], and OC20 [37]. QM9 consists of 134k small organic molecules consist of
carbon, hydrogen, nitrogen, oxygen, and fluorine in equilibrium state. All
molecular conformations and their corresponding properties were created by
computational simulation based on DFT calculation. Dataset contains each
molecule’s stable three-dimensional coordinate and 12 scalar quantum chemical
properties consist of geometric, energetic, electronic, and thermodynamic
properties of given molecular structure. The task of the datasets are
predicting properties from molecular conformation. Since the dataset contains
comprehensive chemicals with high consistence [36], many molecular property
prediction tasks were evaluated on this dataset.

The original MD17 dataset [17] consists of 10 small molecules. For each
molecule, dataset contains energy and forces for different geometry with more
than 10k structures, based on DFT calculation. Instead of using original
dataset we used dataset consist of 1,000 training sample [19] with more
precise calculation than original MD17 dataset. The task of the datasets are
predicting forces and energy of given conformation, from other conformations
of specific molecule.

Open Catalyst 2020 (OC20) is an open-source for learning catalysis dynamical
property with chemical configuration. Similar to the benchmarks mentioned
above, OC20 data consists of a three-dimensional configurations with atomic
number of each adsorbate onto surfaces in the initial and DFT relaxation
states. Nevertheless, OC20 is a more advanced task since the relaxed state
between adsorbates and surfaces needs to be predicted the initial state
thereof. OC20 includes three different types of tasks: structure to energy and
forces (S2EF), initial structure to relaxed structure (IS2RS) and initial
structure to relaxed energy (IS2RE) tasks. Each task is subdivided into
several tasks according to dataset size. We focused on the IS2RE tasks on 10k
and full sized-datasets, which needs relatively less computations.

###  5.2 Implementation

In this paper, we stacked up eight Geometry-aware Transformer layers ($L=8$)
for QM9 prediction, and four ($L=4$) for MD17 and OC20 dataset. We used 300
($n_{basis}=300$) different Gaussian functions with $\mu_{k}=0.1\cdot k$ to
embed atom-atom distances. For MLPs in RBF, We set the hidden dimension to 64
($d_{rbf}=64$). All embedded representations ($d_{e}$), model dimensions
($d_{m}$) and all hidden dimensions ($d_{h}$) of the GeoT layers were selected
between (256, 512, 1,024). with (4, 16, 64) heads ($h=4,16,64$). We used swish
activation function [38] for RBF and ELU activation function [32] for feed
forward network and readout phase.

###  5.3 Training

For QM9 training, we removed 3k data which reported as structures are unstable
following previous works for fair comparison [6, 7, 13, 8, 9]. Each label was
trained individually. Mean absolute error (MAE) was used for evaluation
according to the guidelines. Adam [39] was used as the optimizer with MAE loss
and the batch size was 32 in all experiments. The learning rate was set to
0.0002 with linear warmup 3,000 steps. The learning rate decreased at 0.95 for
every 200k step. We applied the early stopping method by evaluating every 10k
step of training. The number of training epochs was up to 300.

To calculate forces, we assumed that the predicted energy as a function of
positions $E(\\{\mathbf{r}_{i}\\})$. Based on the relation between force and
potential energy of atoms, we calculated forces by differentiating energies
with atom positions as previous works [18, 6].

$F_{j}=\frac{\partial}{\partial\mathbf{r}_{j}}E(\\{\mathbf{r}_{i}\\})$ (15)

In the MD17 dataset, both the energy and forces of molecular conformations are
provided. To utilize both features with Eq. 15, we trained our model with the
modified loss function with additional force terms as Eq. 16 [18, 6].

$L(E,\hat{E})=|E-\hat{E}|+c\times\sum_{j=i}^{n}|F_{j}-\frac{\partial}{\partial\mathbf{r}_{j}}\hat{E}(\\{\mathbf{r}_{i}\\})|$
(16)

Where $c$ is a weight coefficient for the loss on forces. We set $c$ as 1,000
in our experiments.

For training OC20 dataset, the trained model was evaluated on the in-
distribution subset for validation, following the official guideline [37].
Note that we did not implement the periodic boundary condition (PBC) trick,
whereas all other methods did. PBC refers the approximation technique for
analysis in large system with small repeated pattern. With PBC method,
multiple replica of a small unique pattern are arranged to represent a regular
system. By virtue of this, broad regular surface molecules can be represented
with a much smaller-sized unique cell. It significantly lowers the
computations and ensures consistent representations of different atoms in same
surface. We are preparing further experimental results for the next version of
our manuscript.

##  6 Result

Table 1: Mean absolute error on QM9

Target | Unit | SchNet | PhysNet | GRAT | Cormorant | DimeNet++ | GeoT
---|---|---|---|---|---|---|---
$\mu$ | D | 0.033 | 0.0529 | 0.03898 | 0.038 | 0.0297 | 0.0297∗
$\alpha$ | ${a_{0}}^{3}$ | 0.235 | 0.0615 | 0.07219 | 0.085 | 0.0435 | 0.05269
$\epsilon_{\mathrm{HOMO}}$ | eV | 0.041 | 0.0329 | 0.02228 | 0.034 | 0.0246 | 0.0250∗
$\epsilon_{\mathrm{LUMO}}$ | eV | 0.034 | 0.0247 | 0.02053 | 0.038 | 0.0195 | 0.0202∗
$\Delta_{\epsilon}$ | eV | 0.063 | 0.0425 | 0.035 | 0.061 | 0.0326 | 0.04392
<$R^{2}$> | ${a_{0}}^{2}$ | 0.073 | 0.765 | 0.76681 | 0.961 | 0.331 | 0.3008
$zpve$ | eV | 0.0017 | 0.0014 | 0.00208 | 0.0020 | 0.00121 | 0.00173∗
$U_{0}$ | eV | 0.014 | 0.0082 | 0.0705 | 0.022 | 0.00632 | 0.0111
$U$ | eV | 0.019 | 0.0083 | 0.02825 | 0.021 | 0.0062 | 0.0117
$H$ | eV | 0.014 | 0.0084 | 0.02549 | 0.023 | 0.0065 | 0.0113
$G$ | eV | 0.014 | 0.0094 | 0.02505 | 0.020 | 0.0075 | 0.0117
$C_{v}$ | cal/mol K | 0.033 | 0.028 | 0.0326 | 0.026 | 0.023 | 0.0276

  * •

∗ $\mu$, $\epsilon_{\mathrm{HOMO}}$, $\epsilon_{\mathrm{LUMO}}$, and $zpve$
were evaluated on $L=16$, which is different from other targets ($L=8$).

###  6.1 Model performance

We compared our result with five previous studies. Four studies [6, 8, 7, 13]
are MPNN-based studies, and [9] is the only Transformer-based study. Note that
all methods but SchNet [6] uses both distances and angles between atoms. Our
model outperformed previous models which are only used distance between atoms.
However, the performances from our model did not exceed those of DimeNet++
[13], which adopted angular information. We achieved comparable performances
with DimeNet++ among six targets and outperformed for one. We achieved better
performances than those from GRAT, which also adopted the Transformer
architecture for QM9, except three targets.

Table 2: Basis functions for ablation study Type | Basis function |
---|---|---
Linear basis | $a+b\cdot d_{ij}$ |
Gaussian basis | exp($-\gamma||d_{ij}-d_{k}||^{2}$) |
Bessel basis | $\sqrt{\frac{2}{c}}\frac{\sin(\frac{n\pi}{c}d)}{d}$ |
Table 3: Ablation study for basis functions in Table 2 in QM9 Target | Unit | Gaussian basis | Linear basis | Bessel basis
---|---|---|---|---
$\mu$ | D | 0.0299 | 0.664 | 0.0876
$\alpha$ | ${a_{0}}^{3}$ | 0.05269 | 0.752 | 0.15
$\epsilon_{\mathrm{HOMO}}$ | eV | 0.02675 | 0.0659 | 0.056
$\epsilon_{\mathrm{LUMO}}$ | eV | 0.02145 | 0.122 | 0.0326
$\Delta_{\epsilon}$ | eV | 0.04392 | 0.158 | 0.0575
<$R^{2}$> | ${a_{0}}^{2}$ | 0.3008 | 55.4 | 0.3262
$zpve$ | eV | 0.00174 | 0.00845 | 0.00255
$U_{0}$ | eV | 0.0111 | 0.2026 | 0.0147
$U$ | eV | 0.0117 | 51.99 | 0.0275
$H$ | eV | 0.0113 | 0.952 | 0.06
$G$ | eV | 0.0117 | 0.574 | 0.013
$C_{v}$ | cal/mol K | 0.0276 | 0.195 | 0.06
Table 4: Force Prediction MAE in MD17

kcal/mol/A | sGDML | DimeNet | GemNet-T | GeoT | GeoT+A | GeoT+B | GeoT+C | GeoT+B+C | GeoT+A+B+C
---|---|---|---|---|---|---|---|---|---
Aspirin | 0.6803 | 0.4981 | 0.2191 | 0.7081 | 0.7147 | 0.7657 | 0.606 | 0.741 | 0.721
Benzene[9] | - | 0.1868 | 0.1453 | - | 0.10879 | 0.109 | 0.1132 | 0.1206 | 0.0978
Ethanol | 0.3298 | 0.2306 | 0.0853 | 0.09494 | 0.093508 | 0.08967 | 0.0952 | 0.07741 | 0.08599
Malonaldehyde | 0.4105 | 0.3828 | 0.1545 | 0.13985 | - | 0.1408 | 0.14215 | 0.1476 | 0.151
Naphthalene | 0.1107 | 0.2145 | 0.0553 | 0.2198 | 0.2264 | 0.1947 | 0.2321 | 0.1949 | 0.1924
SalicylicAcid | 0.2790 | 0.3736 | 0.1268 | 0.2869 | 0.289 | 0.2849 | 0.2758 | 0.2755 | 0.2667
Toluene | 0.1407 | 0.2168 | 0.0600 | 0.1541 | 0.1538 | 0.1307 | 0.1642 | 0.1403 | 0.147
Uracil | 0.2398 | 0.3021 | 0.0969 | 0.1471 | 0.1499 | 0.1341 | 0.1631 | 0.13097 | 0.1285

  * •

A means that the use of atom-aware two-body kernel, which is described in Eq.
7), otherwise, Eq. 6 was used instead. B and C are the use of Paralleled MLP
and the AttnScale, respectively.

We executed ablation studies on various types of radial basis functions. We
used other two types of radial bases: Bessel basis functions, adopted in
DimeNet [10, 13] and simple linear basis functions. The formula is described
in 2. We also examined for linear basis function, which is most simple case of
RBF. As shown in Table 3, we found that Gaussian basis functions used in our
model significantly improved the model performances.

The performance from MD17 is described in Table 4. We compared our results
with previously reported methods [18, 10, 20]. GeoT achieved the best
performance in three molecules, whereas GemNet-T outperformed in the rest five
molecules. Note that DimeNet uses angle and GemNet-T uses torsional angle
features in a molecule. We did another ablation study with two types of model
refinements mentioned above and the existence of atom embeddings. The model
refinements improved the model performance in six out of eight types of
molecules, except for Aspirin and Malonaldehyde.

Table 5 and 6 showed the performances of OC20 IS2RE task of 10k and full
datasets. The baseline results of [6, 10, 40, 13] were retrieved from the open
leaderboard. We figured out that implementing PBC clearly improves the model
performance in SchNet. Without PBC, SchNet does not exceed our model in both
of the tasks. The model refinements did not effective in 10k task. We remain
the issue with open question.

Table 5: In-distribution validation energy MAE for OC20 IS2RE 10k task

model | MAE
---|---
SchNet* | 1.059
SchNet (without PBC) | 1.077
DimeNet* | 1.0117
CGCNN* | 0.9881
DimeNet++* | 0.8837
GeoT | 0.9631
GeoT + Atom aware kernel | 0.996
GeoT + ParallelMLP | 0.966
GeoT + AttnScale | 1.018
GeoT + Attnscale + ParallelMLP | 1.0516
GeoT + Attnscale + ParallelMLP + Atom aware kernel | 0.9786

  * •

The results with * were retrieved from [37]. Note that we did not implement
PBC conditions into GeoT, while other methods adopted it.

Table 6: In-distribution validation energy MAE for OC20 IS2RE full task model | MAE
---|---
SchNet | 0.6458
SchNet (without PBC) | 0.729
GeoT (without PBC) | 0.7002

(a) Short-range local attention

(b) Long-range attention

(c) Dispersed attention

(d) Localized attention

Figure 3: The four sample pairs of the attention map of QM9 model. Each pair
is retrieved from the Second (left) and Seventh GeoT layer (right),
respectively. The red asterisk is the arbitrarily selected query atom of a
given molecule. The attention maps varies according to molecule types. (a) and
(b) are the examples of short-range and long-range attention from query atom,
respectively. (c) and (d) are the distributed and localized attention,
respectively. All model is trained on HOMO.

(a) Salicylic acid

(b) Uracil

Figure 4: The analysis of the model refinements and atom embedding effects in
MD17. Note that GeoT showed better performances accompanied with AttnScale,
paralleld MLP, and atom embeddings in two-body kernel, in Salicylic acid and
Uracil predictions. In both of the cases, the left parts were extracted from
the model without refinements while the right parts from the model with all
refinements. We found that attentions become more specified and the prediction
performance was improved, with the additional strategies.

(a) Aspirin

(b) Toluene

Figure 5: The analysis of the attention maps in each layer of GeoT. In both
cases, the attention values of right parts (the fourth layer) are more
specified and localized than those of left parts (the first layer).

###  6.2 Attention map visualization

We visualized attention maps from the trained model to analyze the behavior of
our attention mechanism. First, we will describe how we extract attention
values. Next, we will discuss two subjects in this section, the aspects of
attentions and the effects of model refinements.

Different from the previous work [41], we extracted
$||\mathrm{Attn}^{\prime}_{G}(\bm{Q},\bm{K},\bm{\Lambda})||_{2}$ of several
molecules from the trained model. We picked a query atom $i$ of each molecule,
and visualized attention values with another atom $j$. Query atoms are marked
with red asterisk, and attention norms of key atom $j$ from the query atom $i$
are colored with blue shades.

The attention values in GeoT are differ from the original Transformer in two
points. First, GeoT updates each atom representation with unnormalized
attention scale, given by Eq. 8. Second, $Attn^{\prime}_{G}$ includes a novel
feature $\bm{\Lambda}$, which endows a geometric feature with the attention.
By capturing unnormalized geometric attenion values,
$||Attn^{\prime}_{G}||_{2}$ provide us an intuitive interpretation that the
intensity of one atom $i$ attend another atom $j$ with their distance
$d_{ij}$.

Especially, we focused on two aspects of $||Attn_{G}^{\prime}||_{2}$, the
$range$ and the $dispersity$. In GeoT, an atom can attend any other atom
without cutoff, so attention values exist all atoms in a molecule. Note that
the single bond between two carbons are about 1.54 Å, and the diameter of
benzene rings are approximately 3 Å. Since existing message-passing model
restricts inside cutoff boundary to 5Å, they cannot cover neighboring two
successive benzene rings or three carbon-carbon bonds far from an atom $i$.
Meanwhile, attention values should also be distinguished. It means that the
model should focus on more interactive atom pairs. It is not contrary to
cutoff-free attention: attention values might range globally but can be more
concentrated around the query atom if no informative atoms exist in far away.

We found that the attention values become more significant as they locate
closer to the query atom (Figure 3(a), 3(b)). At the same time, these values
are more evenly dispersed around in the case of aromatic rings (Figure 3(c),
3(d)). In consideration of the characteristics of aromaticity, these pattern
suggest that GeoT learnt some concepts of Chemistry.

We did a further analysis on the effect of model refinement tricks described
in 4.4. In Figure 4, we visualized salicylic acid and uracil molecule in MD17
dataset, respectively. The left part of each pair shows the attention map from
GeoT trained without any model refinements, whereas the right part contains
AttnScale, paralleled MLP and atom embedding. Applying model refinements
clearly effective on both of the molecules (Table 4). In general, attention
values seems to be more distributed with model refinements, and capture the
aromatic ring structure. Regardless of whether model refinements are adopted
or not, attention values become more specified in deeper layers, which might
be indicate suppression of over-smoothing.

Finally, we analyzed the attention maps of each layer inside the trained GeoT.
We visualize two examples in MD17: Aspirin and Toluene in Figure 5(a) and
5(b), respectively. The left and right parts of each pair are the attention
map extracted from the first and fourth layer, respectively. We found that the
attention map from latter layer is slightly more specified and localized in
both of the cases.

##  7 Discussion

Based on the above results and visualizations, we discuss several issues on
the behavior of GeoT with multiple aspects.

Q1. Is Transformer better than message-passing for molecular graph
representation?
The most significant difference is that the ability of learning representation
between any two atoms in a molecule. In MPNN, messages are only passed to
closely located atoms limited to cutoff values. However, the attention
mechanism (both of the original self-attention and GeoAttention) acts globally
in a molecule so that any two atom configuration can contribute to the whole
molecule representation learning.

Second, the deep message-passing models often suffer from the over-smoothing
problem, while Transformer-based networks showed outstanding performances in
various tasks with deeper structures.

Q2. Is GeoAttention better than the original self-attention for molecular
graph representation?
We emphasize that GeoAttention includes a radial-basis feature
$\mathbf{\Lambda}$ for atom pair distances and removes a softmax function.
First, $\mathbf{\Lambda}$ provides all pairwise atom-atom distances with
multiple bases. We assumed that if all atom-atom pairs with their distances
from all combinations are available, then the model can recover the global
view of its three-dimensional configurations. Second, the scale of
GeoAttention is not normalized so that it is not related with the number of
atoms.

Note that we do not use either message passing framework or graph convolution
in our method. Therefore, our proposed model can capture the global
perspective of a molecular conformation without any local restriction. In
addition, our method is relatively free from over-smoothing problems
regardless of the model depth with additional refinements.

Q3. What insights can we derive from the attention map based on GeoAttention?
The attention visualization provides deeper insights than the prediction
performance only, because we can investigate each atom and their relations in
each molecule. The attention maps showed that the nature of atom-atom
interaction and the aromaticity are properly reflected. Further, it can
suggest unknown chemical interactions in complex system, which can be a
significant clue of synthetic chemistry researches.

##  8 Conclusion

In this study, we introduced geometry-aware attention which can capture
interactions regardless of distances. We added geometric information by
adopting self-attentions of the distance representations between atoms with
the RBF. We evaluated our model on QM9 and MD17. Our model outperformed
previous works based on distances and comparable with state-of-the-art models
which utilized additional angular information. We also studied on the
candidates of the RBF, and showed Gaussian basis function performed
significantly better than other functions. We also displayed attention
weights, which supported our assumption that model can attend atom pairs
located far from each other. We argue that our model achieved that long-range
interactions between any intramolecular atom pairs, which cannot be available
for previous methods using the restricted convolutional filters. We will
explore additional geometric features to improve our model performance. In the
future, we will explore more extensive experiments on OC20, including PBC and
more variable tasks.

## References

  * [1] Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller, and O Anatole Von Lilienfeld.  Fast and accurate modeling of molecular atomization energies with machine learning.  Physical review letters, 108(5):058301, 2012.
  * [2] Sönke Lorenz, Axel Groß, and Matthias Scheffler.  Representing high-dimensional potential-energy surfaces for reactions at surfaces by neural networks.  Chemical Physics Letters, 395(4-6):210–215, 2004.
  * [3] Jörg Behler and Michele Parrinello.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  Physical review letters, 98(14):146401, 2007.
  * [4] 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.  Physical review letters, 104(13):136403, 2010.
  * [5] David Duvenaud, Dougal Maclaurin, Jorge Aguilera-Iparraguirre, Rafael Gómez-Bombarelli, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P Adams.  Convolutional networks on graphs for learning molecular fingerprints.  arXiv preprint arXiv:1509.09292, 2015.
  * [6] Kristof T Schütt, Pieter-Jan Kindermans, Huziel E Sauceda, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller.  Schnet: A continuous-filter convolutional neural network for modeling quantum interactions.  arXiv preprint arXiv:1706.08566, 2017.
  * [7] Oliver T Unke and Markus Meuwly.  Physnet: A neural network for predicting energies, forces, dipole moments, and partial charges.  Journal of chemical theory and computation, 15(6):3678–3693, 2019\.
  * [8] Brandon Anderson, Truong-Son Hy, and Risi Kondor.  Cormorant: Covariant molecular neural networks.  arXiv preprint arXiv:1906.04015, 2019.
  * [9] Sanghyun Yoo, Young-Seok Kim, Kang Hyun Lee, Kuhwan Jeong, Junhwi Choi, Hoshik Lee, and Young Sang Choi.  Graph-aware transformer: Is attention all graphs need?  arXiv preprint arXiv:2006.05213, 2020.
  * [10] Johannes Klicpera, Janek Groß, and Stephan Günnemann.  Directional message passing for molecular graphs.  arXiv preprint arXiv:2003.03123, 2020.
  * [11] Yoni Choukroun and Lior Wolf.  Geometric transformer for end-to-end molecule properties prediction.  arXiv preprint arXiv:2110.13721, 2021.
  * [12] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In International Conference on Machine Learning, pages 1263–1272. PMLR, 2017.
  * [13] Johannes Klicpera, Shankari Giri, Johannes T. Margraf, and Stephan Günnemann.  Fast and uncertainty-aware directional message passing for non-equilibrium molecules.  In NeurIPS-W, 2020.
  * [14] Qimai Li, Zhichao Han, and Xiao-Ming Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
  * [15] Deli Chen, Yankai Lin, Wei Li, Peng Li, Jie Zhou, and Xu Sun.  Measuring and relieving the over-smoothing problem for graph neural networks from the topological view.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 3438–3445, 2020.
  * [16] Wenbing Huang, Yu Rong, Tingyang Xu, Fuchun Sun, and Junzhou Huang.  Tackling over-smoothing for general graph convolutional networks.  arXiv preprint arXiv:2008.09864, 2020.
  * [17] Stefan Chmiela, Alexandre Tkatchenko, Huziel E Sauceda, Igor Poltavsky, Kristof T Schütt, and Klaus-Robert Müller.  Machine learning of accurate energy-conserving molecular force fields.  Science advances, 3(5):e1603015, 2017.
  * [18] Stefan Chmiela, Huziel E Sauceda, Igor Poltavsky, Klaus-Robert Müller, and Alexandre Tkatchenko.  sgdml: Constructing accurate and data efficient molecular force fields using machine learning.  Computer Physics Communications, 240:38–45, 2019.
  * [19] Stefan Chmiela, Huziel E Sauceda, Klaus-Robert Müller, and Alexandre Tkatchenko.  Towards exact molecular dynamics simulations with machine-learned force fields.  Nature communications, 9(1):1–10, 2018.
  * [20] Johannes Gasteiger, Florian Becker, and Stephan Günnemann.  Gemnet: Universal directional graph neural networks for molecules.  Advances in Neural Information Processing Systems, 34:6790–6802, 2021.
  * [21] Yi Liu, Limei Wang, Meng Liu, Yuchao Lin, Xuan Zhang, Bora Oztekin, and Shuiwang Ji.  Spherical message passing for 3d molecular graphs.  In International Conference on Learning Representations, 2022.
  * [22] Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley.  Tensor field networks: Rotation-and translation-equivariant neural networks for 3d point clouds.  arXiv preprint arXiv:1802.08219, 2018.
  * [23] Fabian B Fuchs, Daniel E Worrall, Volker Fischer, and Max Welling.  Se (3)-transformers: 3d roto-translation equivariant attention networks.  arXiv preprint arXiv:2006.10503, 2020.
  * [24] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  arXiv preprint arXiv:1706.03762, 2017.
  * [25] Peihao Wang, Wenqing Zheng, Tianlong Chen, and Zhangyang Wang.  Anti-oversmoothing in deep vision transformers via the fourier domain analysis: From theory to practice.  arXiv preprint arXiv:2203.05962, 2022.
  * [26] Sheng Wang, Yuzhi Guo, Yuhong Wang, Hongmao Sun, and Junzhou Huang.  Smiles-bert: large scale unsupervised pre-training for molecular property prediction.  In Proceedings of the 10th ACM international conference on bioinformatics, computational biology and health informatics, pages 429–436, 2019.
  * [27] Łukasz Maziarka, Tomasz Danel, Sławomir Mucha, Krzysztof Rataj, Jacek Tabor, and Stanisław Jastrzębski.  Molecule attention transformer.  arXiv preprint arXiv:2002.08264, 2020.
  * [28] Yu Rong, Yatao Bian, Tingyang Xu, Weiyang Xie, Ying Wei, Wenbing Huang, and Junzhou Huang.  Self-supervised graph transformer on large-scale molecular data.  Advances in Neural Information Processing Systems, 33:12559–12571, 2020.
  * [29] Chengxuan Ying, Tianle Cai, Shengjie Luo, Shuxin Zheng, Guolin Ke, Di He, Yanming Shen, and Tie-Yan Liu.  Do transformers really perform badly for graph representation?  Advances in Neural Information Processing Systems, 34:28877–28888, 2021.
  * [30] Dexiong Chen, Leslie O’Bray, and Karsten Borgwardt.  Structure-aware transformer for graph representation learning.  In International Conference on Machine Learning, pages 3469–3489. PMLR, 2022.
  * [31] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.  Layer normalization.  arXiv preprint arXiv:1607.06450, 2016.
  * [32] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  arXiv preprint arXiv:1511.07289, 2015.
  * [33] Namuk Park and Songkuk Kim.  How do vision transformers work?  arXiv preprint arXiv:2202.06709, 2022.
  * [34] Han Shi, Jiahui Gao, Hang Xu, Xiaodan Liang, Zhenguo Li, Lingpeng Kong, Stephen Lee, and James T Kwok.  Revisiting over-smoothing in bert from the perspective of graph.  International Conference on Learning Representations, 2022.
  * [35] Lars Ruddigkeit, Ruud Van Deursen, Lorenz C Blum, and Jean-Louis Reymond.  Enumeration of 166 billion organic small molecules in the chemical universe database gdb-17.  Journal of chemical information and modeling, 52(11):2864–2875, 2012.
  * [36] Raghunathan Ramakrishnan, Pavlo O Dral, Matthias Rupp, and O Anatole von Lilienfeld.  Quantum chemistry structures and properties of 134 kilo molecules.  Scientific Data, 1, 2014.
  * [37] Lowik Chanussot*, Abhishek Das*, Siddharth Goyal*, Thibaut Lavril*, Muhammed Shuaibi*, Morgane Riviere, Kevin Tran, Javier Heras-Domingo, Caleb Ho, Weihua Hu, Aini Palizhati, Anuroop Sriram, Brandon Wood, Junwoong Yoon, Devi Parikh, C. Lawrence Zitnick, and Zachary Ulissi.  Open catalyst 2020 (oc20) dataset and community challenges.  ACS Catalysis, 2021.
  * [38] Prajit Ramachandran, Barret Zoph, and Quoc V Le.  Searching for activation functions.  arXiv preprint arXiv:1710.05941, 2017.
  * [39] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [40] Tian Xie and Jeffrey C Grossman.  Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties.  Physical review letters, 120(14):145301, 2018.
  * [41] Jesse Vig.  A multiscale visualization of attention in the transformer model.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics: System Demonstrations, pages 37–42, Florence, Italy, July 2019. Association for Computational Linguistics.