# Global performance of covariant density functional theory in description of
charge radii and related indicators.

A.V.Afanasjev  U.C.Perera  P. Ring  Department of Physics and Astronomy,
Mississippi State University, MS 39762  Fakultät für Physik, Technische
Universität München, D-85748 Garching, Germany

(Day Month Year (Insert date of submission))

Afanasjev, Perera and RingGlobal performance of CDFT in description of charge
radii and related indicators

1, 1, 2

1 2

A short review of existing efforts to understand charge radii and related
indicators on a global scale within the covariant density functional theory
(CDFT) is presented. Using major classes of covariant energy density
functionals (CEDFs), the global accuracy of the description of experimental
absolute and differential charge radii within the CDFT framework has been
established. This assessment is supplemented by an evaluation of theoretical
statistical and systematic uncertainties in the description of charge radii.
New results on the accuracy of the description of differential charge radii in
deformed actinides and light superheavy nuclei are presented and the role of
octupole deformation in their reproduction is evaluated. Novel mechanisms
leading to odd-even staggering in charge radii are discussed. Finally, we
analyze the role of self-consistency effects in an accurate description of
differential charge radii.

Covariant density functional theory, charge radii, differential charge radii,
theoretical uncertainties, self-consistency effects

##  1 Introduction

The charge radii are among the most fundamental properties of atomic nuclei.
Essential information on the saturation density of symmetric nuclear matter is
imprinted into them. They also depend on the properties of nuclear forces and
nuclear many-body dynamics. Charge radii have been in the focus of
experimental and theoretical studies for a long time and, in the last decade,
there was a significant increase in the interest to them due to the
availability of new facilities and new techniques. Absolute charge radii are
determined via the studies of muonic spectra and electronic scattering
experiments [1, 2]. However, such experiments are impossible in radioactive
elements. Laser spectroscopy is a more flexible tool in that respect and
allows a significantly higher precision of the measurements of the changes of
charge radii within the isotopic chain [1, 2].

The CDFT [3] is one of the modern theoretical tools which has been used with
reasonable success to the description of charge radii and related indicators.
For example, the first-ever successful description of the kink in charge radii
of the lead isotopes has been achieved in the CDFT framework in Ref. [4] and
in the 90ies it followed by a study of differential charge radii in spherical
even-even nuclei of the Ca, Sn, and Pb isotopic chains in Ref. [5]. However,
only during the last seven years did new CDFT results appear, dedicated to the
study of charge radii, related indicators, and other connected issues. The
goal of this paper is to present a short review of different facets of physics
of charge radii in CDFT and to provide new clarifying results. All the quoted
and new calculations are performed within the framework of relativistic
Hartree-Bogolibov (RHB) approach with separable pairing interaction of Ref.
[6]. The technical details of the RHB approach and numerical details of the
calculations are presented in Ref. [7]. The calculations are carried out with
the CEDFs (such as NL3* [8], DD-ME2 [9], DD-ME$\delta$ [10], DD-PC1 [11], PC-
PK1 [12] and NL5(C), NL5(D) and NL5(E) [13]) representing the major classes of
the CDFT models.

##  2 Physical observables

We focus here on physical observables, which provide information on charge
radii and related indicators. The charge radii were calculated from the
corresponding point proton radii as (see Ref. [7])

$r_{ch}=\sqrt{<r^{2}>_{p}+0.64}\,\,\,\,{\rm fm}$ (1)

where $<r^{2}>_{p}$ stands for proton mean-square point radius and the factor
0.64 accounts for the finite-size effects of the proton.

In addition, two differential indicators are commonly used to facilitate the
quantitative comparison of the experimental results with those from
theoretical calculations. One of them is the differential mean-square charge
radius

$\displaystyle\delta\left<r^{2}\right>_{p}^{N,N^{\prime}}=r^{2}_{ch}(N)-r^{2}_{ch}(N^{\prime})$
(2)

where $N^{\prime}$ is the neutron number of the reference nucleus. Another is
the three-point indicator

$\displaystyle\Delta\langle
r^{2}\rangle^{(3)}(N)=\frac{1}{2}\left[r_{ch}^{2}(N-1)+r_{ch}^{2}(N+1)-2r_{ch}^{2}(N)\right]$
(3)

which quantifies OES in charge radii.

##  3 Global performance of the CDFT in the description of charge radii

Figure 1: The difference between measured and calculated charge radii $r_{ch}$
for the CEDF DD-PC1. Two-proton and two-neutron drip lines obtained with this
functional are shown by solid black lines. The experimental data are taken
from Ref. [1]. Based on figure from Ref. [14].

The global performance of CDFT in the description of charge radii for major
classes of CEDFs has been studied in Refs. [7, 14, 13], and Fig. 1 presents an
example of a comparison between theory and experiment. The following rms
deviations $\Delta r_{ch}^{rms}$ between calculated and experimental charge
radii have been reported: $\Delta r_{ch}^{rms}=0.0283$ fm for CEDF NL3*,
$\Delta r_{ch}^{rms}=0.0230$ fm for DD-ME2, $\Delta r_{ch}^{rms}=0.0329$ fm
for DD-ME$\delta$*, $\Delta r_{ch}^{rms}=0.0253$ fm for DD-PC1, $\Delta
r_{ch}^{rms}=0.0284$ fm for NL5(C), $\Delta r_{ch}^{rms}=0.0277$ fm for
NL5(D), and $\Delta r_{ch}^{rms}=0.0288$ fm for NL5(E). Note that the DD-PC1
functional has been adjusted only to nuclear binding energies while other
functionals include the information on charge radii in their fitting
protocols. The average rms deviations of calculated $r_{ch}$ from experimental
ones for this group of functionals are at the level of 0.028 fm. Considering
that the average experimental rms charge radius in the nuclear chart is around
4.8 fm (see, for example, Fig. 23 in Ref. [7] and Figs. 2-4 in Ref. [1]), this
amounts to a high average precision of 0.58% in the prediction of charge
radii.

Fig. 1 reveals several regions on the nuclear chart which exhibit enhanced
(with respect of average trend) differences $(r_{ch}^{exp}-r_{ch}^{th})$
between experimental and calculated charge radii. These are (i) proton-rich
nuclei around Hg ($Z=80$), (ii) proton-rich nuclei around Sr ($Z=38$) located
below the $N=50$ shell closure and (iii) light and very light nuclei in the
region with $Z<20$. The detailed analysis of the nuclei in the regions (i) and
(ii) performed in Refs. [15, 16, 17, 18, 19] clearly indicates that shape
coexistence and the correlations beyond mean-field are the sources of these
discrepancies. Note that the latter are neglected in our studies. The
importance of the correlations beyond mean-field is also expected to be
increased in light and very light nuclei, the potential energy surfaces of
which are very frequently soft in collective coordinates. The fact that the
accuracy of the description of the masses declines with decreasing mass also
points in that direction (see discussion in Sec. V of Ref. [7]). Thus, the
inclusion of the correlations beyond mean-field in the regions (i-iii) may
improve the description of charge radii and binding energies. The accuracy of
the description of charge radii also deteriorates in the actinides (see Fig.
1). However, this is most likely due to the extrapolations used in Ref. [1]
for the definition of absolute charge radii since they have not been measured
experimentally in the $Z>83$ nuclei (uranium is an exception).

##  4 Statistical and systematic uncertainties in the description of charge
radii

Because of underlying physical approximations in modeling of the nuclear many-
body problem it is important to estimate statistical and systematic
theoretical uncertainties in the predictions of charge radii. This is
especially important when one deals with the extrapolations beyond the known
regions, as, for example, in particle number or deformation, since
experimental data which act as a substitute of an exact solution are not
available there. Such estimates are also required for the evaluation of the
predictive power of the models and the robustness of their predictions.

Systematic theoretical uncertainties emerge from underlying theoretical
approximations about the form of the energy density functional. In the DFT
framework, there are two major sources of these approximations, namely, the
range of interaction and the form of the density dependence of the effective
interaction [20, 21, 3]. In reality, only the lower limit of these
uncertainties is possible to estimate (see Refs. [7, 22, 13]) and even then,
it contains a degree of arbitrariness related to the choice of the set of the
functionals. In our works, it is estimated via the spread of theoretical
predictions within a set of selected functionals (NL3*, DD-ME2, DD-PC1, DD-
ME$\delta$) as

$\Delta r_{ch}(Z,N)=|r_{ch}^{max}(Z,N)-r_{ch}^{min}(Z,N)|$ (4)

where $r_{ch}^{max}(Z,N)$ and $r_{ch}^{min}(Z,N)$ are the largest and smallest
values of the charge radii obtained with the four employed CEDFs for the
$(Z,N)$ nucleus [7]. The charge radii spreads are shown in Fig. 2. For the
nuclei shown in Fig. 1 they are on the level of 0.01 fm for the majority of
medium mass and heavy nuclei but become larger for $Z\leq 40$ nuclei where
they exceed 0.02 fm in most of the nuclei. The latter feature is most likely
due to softer potential energy surfaces of such nuclei and thus stronger
dependence of the results on the details of the underlying single-particle
structure. As a result, in some nuclei the systematic theoretical
uncertainties are either comparable or larger than the differences between
experimental and calculated charge radii (compare Figs. 2 and 1). This fact
has to be kept in mind when performing a detailed comparison between theory
and experiment.

Figure 2: Systematic uncertainties in the description of charge radii. From
Ref. [7].

Figure 3:  The propagation of statistical errors in charge radii with neutron
number for the Ca (Z = 20), Ni (Z = 28), Sn (Z = 50), and Pb (Z = 82) isotopic
chains. The RHB calculations have been performed at spherical shape with the
NL5(C) functional for all even-even nuclei located between the two-proton and
two-neutron drip lines (see Refs. [13] for details). Open symbols are used to
indicate the nuclei whose experimental data have been used in the fitting
protocol of this functional. From Ref. [13].

Statistical errors $\sigma(r_{ch})$, which apply only to a given functional,
are due to the details of the fitting protocol, such as the choice of
experimental data and the selection of adopted errors (see Refs. [22, 13]).
They are shown in Fig. 3 for the NL5(C) functional for several isotopic chains
of spherical nuclei. However, similar magnitude of $\sigma(r_{ch})$ is
expected both for other functionals and for deformed nuclei [13]. One can see
that, on average, they are by a factor of approximately 6 smaller than global
rms deviations between theory and experiment. They are also substantially
smaller than systematic theoretical uncertainties for charge radii. Note that
the $\sigma(r_{ch})$ values show some fluctuations as a function of neutron
number due to the underlying shell structure; they become especially
pronounced in the Ni isotopes. These fluctuations are expected to wash out in
deformed nuclei because of a more equal distribution of deformed single-
particle states originating from different spherical subshells. On average,
the $\sigma(r_{ch})$ values shown in Fig. 3 display a very modest increase on
going from the two-proton to the two-neutron drip line. This is in contrast to
the results obtained in the Skyrme DFT calculations with the UNEDF0 functional
(see Ref. [23]) which for the nuclei near the proton drip line shows
$\sigma(r_{ch})$ values comparable to those obtained in the CDFT(NL5(C))
calculations, but are by a factor of 17-33 larger at the neutron drip line
(see Ref. [13]). Based on these considerations, one can conclude that
statistical errors represent only a relatively small part of the total
theoretical error in the prediction of charge radii and in many cases they can
be neglected.

##  5 Differential charge radii: the examples of the Sn/Gd region and the
region of actinides and light superheavy nuclei

The differential charge radius $\delta\left<r^{2}\right>^{N,N^{\prime}}$ (see
Eq.(2)) is an indicator which is very sensitive to the shell closures,
underlying single-particle structure, pairing interaction, and occupation
pattern of the single-particle states and their evolution with neutron number
(see Refs. [15, 24] and references quoted therein). Their systematic global
investigation has been performed within the CDFT framework in Ref. [15] and
theoretical results were compared with experimental data which crosses the
neutron shell closures at N = 28, 50, 82, and 126. An example of such a
comparison is presented in Fig. 4. Absolute differential radii of different
isotopic chains and their relative properties are reasonably well described in
model calculations in the cases when the mean-field approximation is
justified. The observed clusterization of differential charge radii of
different isotopic chains (see Ref. [25]) is well described above the $N=50$
and $N=126$ shell closures (see Ref. [15]). However, a similar clusterization
above the $N=28$ and $N=82$ shell closures is not that well reproduced due to
deficiencies in the underlying single-particle structure.

Figure 4: The experimental and calculated differential charge radii
$\delta\left<r^{2}\right>^{N,N^{\prime}}$ of the indicated isotopic chains.
The experimental data are taken from Ref. [1]. The calculations are performed
with the indicated CEDFs. Theoretical results are shown by lines with symbols
and by only lines for spherical and deformed ground states, respectively.

Figure 5: The experimental and calculated differential charge radii
$\delta\left<r^{2}\right>^{N,N^{\prime}}$ of the indicated isotopic chains.
The results of reflection symmetric and reflection asymmetric (with octupole
deformation) RHB calculations are shown by solid and dashed lines,
respectively. Note that the latter are based on the results obtained in Ref.
[26] and are shown only for the nuclei which have non-zero octupole
deformation $\beta_{3}$. The arrows indicate the reference nuclei in the
isotopic chains. The experimental data are taken from Refs. [1, 27, 28].

In the present paper, we focus on differential charge radii of deformed
actinides and light superheavy nuclei which were not covered in Ref. [15]. The
aim is to understand the accuracy of the description of the
$\delta\left<r^{2}\right>^{N,N^{\prime}}$ values in high-$Z$ systems and the
impact of octupole deformation on differential charge radii. Note that due to
the difficulties of the measurements of absolute charge radii in the $Z>83$
nuclei, the experimental data covers only differential charge radii [the only
exception is uranium]. The experimental and calculated differential charge
radii of the Th, U, Pu, Cm, and No isotopic chains are compared in Fig. 5.
Note that only even-even nuclei are considered in the calculations. The format
of the figure is similar to that of Fig. 2 of Ref. [28] in which the results
of Skyrme DFT calculations with UNEDF1 and SV-min functionals are compared
with experimental data. Note that the RHB calculations indicate the presence
of octupole deformation in low-$N$ Th and U isotopes (see Ref. [26]) but, in a
given nucleus, it somewhat depends on the employed CEDF and on the type of the
DFT [26, 30]. The accounting of octupole deformation leads to an increase of
the $\delta\left<r^{2}\right>^{N,N^{\prime}}$ value by approximately $\approx
0.1$ fm2 and leads to an improvement in the description of differential charge
radii in the Th and U isotopes.

On average, the accuracy of the description of differential charge radii
obtained with covariant (DD-ME2, DD-PC1, PC-PK1, and NL3*) and non-
relativistic energy density functionals (UNEDF1 and SV-min) is similar
(compare Fig. 5 in the present paper with Fig. 2 in Ref. [28]). Among the
employed CEDFs, the best agreement with experiment is obtained for the Th
isotopes by DD-ME2, DD-PC1 and NL3*, for the U isotopes by DD-PC1 and NL3*,
for the Pu isotopes by DD-ME2 and DD-PC1, and for the Cm isotopes by PC-PK1
and NL3* (see Fig. 5). The results of the calculations for the No isotopes are
the same with all employed functionals111Note that the analysis of the No
isotopes performed in Ref. [29] suggests that the measurements of at least two
atomic transitions would allow to disentangle the contributions of the changes
in deformation and charge radius into field isotopic shifts.. The local
deviations between theory and experiment for
$\delta\left<r^{2}\right>^{N,N^{\prime}}$ are most likely due to the
differences between experimental and calculated energies of the deformed
quasiparticle states (see Refs. [31, 32, 33]) which leads to some
dissimilarity of their contributions to charge radii. The same is true when
the results of the calculations with two functionals are compared (see Sect.
IV of Ref. [15] for a discussion of the impact of the underlying single-
particle structure and occupation pattern of the single-particle states on
differential charge radii of the spherical Pb nuclei).

##  6 New mechanisms leading to odd-even staggering in charge radii

The description of odd-even staggering (OES) in charge radii has been a
challenge for a long time. In the most commonly accepted scenario, it is due
to the gradient terms both in the surface part of the energy density
functional and in the pairing interaction (see Refs. [34, 35, 36]). However,
this has recently been challenged in Ref. [37] (see also Ref. [24] for a more
detailed study) where it has been shown that the scattering of the occupation
of different single-particle states between even-even and odd nuclei in the
isotopic chain can lead to OES in charge radii. This scattering is due to the
fact that particle-vibration coupling (PVC) leads to a change of the relative
order of the single-particle states in odd-$N$ nuclei as compared with
even-$N$ ones in which PVC plays a minor role.

The analysis of Refs. [37, 24, 15] also clearly indicates the importance of
taking into account both the structure of the ground state in odd-$A$ nuclei
and the fact that the DFT calculations do not reproduce the structure of the
ground states in odd-$A$ nuclei in more than half of the cases (see Refs. [32,
38, 31, 39]). These facts have been ignored in the absolute majority of the
studies of OES in charge radii. The importance of the accounting of these
facts is illustrated in Fig. 6 in which OES in charge radii of the Sn isotopes
obtained with two procedures (LES and EGS) in the RHB calculations of odd-$N$
nuclei are compared. In the LES procedure, the lowest in energy calculated
configuration is used. In the EGS procedure, the configuration with the spin
and parity of the blocked state corresponding to those of the experimental
ground state is employed, although it is not necessarily the lowest in energy.
The results of the RHB calculations with the LES procedure significantly
underestimate the magnitude of experimental OES and occasionally provide a
wrong phase of the OES (see Fig. 6(b)). The use of the EGS procedure
significantly improves the description of both the phases and the magnitude of
OES (see Fig. 6(a)). However, even in that case, the magnitude of OES is
underestimated on average by a factor of approximately two.

Figure 6: OES in charge radii of the Sn isotopes. The experimental data are
taken from Ref. [1]. The RHB results with the NL3* CEDF are shown by blue
lines with triangles for the LES and EGS procedures in odd-$N$ nuclei. The
results corrected for the fragmentation of the single-particle content of the
dominant single-particle state in odd-$N$ nuclei within the schematic model
discussed in the text are shown by red lines with squares for a representative
spectroscopic factor $S=0.8$.

The remaining discrepancies between theory and experiment suggest that the
energy considerations discussed above are not sufficient and special attention
has to be paid to the wavefunction of the ground state in odd-$N$ nuclei.
There is a substantial difference between the ground states in even-even and
odd-$A$ nuclei with stiff parabola-like potential energy surfaces (such as the
Sn isotopes). Indeed, the impact of the correlations beyond mean-field on the
wavefunction of the ground state is relatively small in even-even nuclei, but
it is substantial in odd-$A$ ones. The latter is due to the coupling of the
single-particle motion with phonons (particle-vibration coupling) which leads
to a substantial fragmentation of the wavefunction; this feature is confirmed
in numerous experiments (see Refs. [40, 38, 39] and references quoted
therein). This reduces the single-particle content of the unpaired neutron
ground state in odd-$A$ nuclei to the value characterized by the spectroscopic
factor $S$ ($S\leq 1.0$). As a result, the pull provided by the odd neutron on
the proton radii (see Refs. [41, 15] and Sec. 7) is reduced as compared with
the mean-field value. Assuming that at the mean-field level the increase of
charge radii on going from even to odd neutron number in the even $Z$ isotopic
chain is $\Delta(r_{ch})_{MF}^{qp}=r_{ch}(N+1)-r_{ch}(N)$ [here $N$ is even],
the accounting of the fragmentation of the wavefunction in odd-$N$ nuclei will
reduce this increase to $\Delta(r_{ch})_{frag}=S\Delta(r_{ch})_{MF}^{qp}$.
Using this schematic model (see Ref. [15] for more details), one can show that
the fragmentation of the structure of the unpaired neutron in odd-$N$ nuclei
due to particle-vibration coupling leads to an increase of the magnitudes of
OES in charge radii and correct phases of the OES both in the LES and EGS
procedures for typical values of $S$ (see Fig. 6). Both of these factors act
in the direction of improving the agreement with experiment.

##  7 The role of self-consistency effects

As reviewed in the introduction of Ref. [15], the absolute majority of the
studies of differential charge radii has been performed either in DFT or in ab
initio approaches. However, it is well known that these types of models have
some deficiencies in the description of spectroscopic properties related to
the energies of the single-particle states and their wave functions [32, 42,
31, 33, 39]. In contrast, the spherical shell models with empirical
interactions provide a better description of experimental spectroscopic data
in spherical nuclei located in the vicinity of doubly magic nuclei and
microscopic+macroscopic (mic+mac) models based on phenomenological potentials
such as the Woods-Saxon one do the same in the region of deformed nuclei.
However, we are not aware of successful investigations of differential charge
radii within these types of models and naturally the question of the
underlying reasons emerges.

Table 1: Proton single-particle rms radii $r_{p,i}=\sqrt{\left<r^{2}\right>_{i}}$ [in $fm$] (columns 4 and 5) and single-particle energies $e_{i}$ [in MeV] (columns 2 and 3) of the indicated single-particle ($i$-th) in the nuclei 208Pb and 220Pb. The pairing is neglected in the calculations and they employ NL3* functional. The column 6 shows the change of proton single-particle radii $\delta r_{p,i}=r_{p,i}(^{220}{\rm{Pb}})-r_{p,i}(^{208}\rm{Pb})$ between the nuclei 208Pb and 220Pb. Note that the results for the single-particle states with principal quantum number $n=1$ are shown in bold.  subshell |  $e_{i}$ (208Pb) |  $e_{i}$ (220Pb) |  $r_{p,i}$ (208Pb) |  $r_{p,i}$ (220Pb) | $\delta r_{p,i}$
---|---|---|---|---|---
1 | 2 | 3 | 4 | 5 | 6
$1s_{1/2}$ | -48.886 | -48.151 | 4.067759 | 4.277595 | 0.209836
$1p_{3/2}$ | -43.193 | -43.155 | 4.667454 | 4.889208 | 0.221754
$1p_{1/2}$ | -42.510 | -42.620 | 4.575735 | 4.805227 | 0.229492
$1d_{5/2}$ | -36.097 | -36.853 | 5.11159 | 5.324115 | 0.212525
$1d_{3/2}$ | -34.533 | -35.636 | 4.972567 | 5.187897 | 0.21533
$2s_{1/2}$ | -30.844 | -32.215 | 4.458931 | 4.544538 | 0.085607
$1f_{7/2}$ | -28.042 | -29.573 | 5.487019 | 5.682079 | 0.19506
$1f_{5/2}$ | -25.265 | -27.366 | 5.325364 | 5.512845 | 0.187481
$2p_{3/2}$ | -20.890 | -22.819 | 5.000532 | 5.032818 | 0.032286
$2p_{1/2}$ | -19.839 | -21.843 | 5.004416 | 5.03297 | 0.028554
$1g_{9/2}$ | -19.369 | -21.587 | 5.82501 | 5.998839 | 0.173829
$1g_{7/2}$ | -15.175 | -18.149 | 5.673305 | 5.82303 | 0.149725
$2d_{5/2}$ | -11.154 | -13.413 | 5.53793 | 5.542713 | 0.004783
$1h_{11/2}$ | -10.335 | -13.120 | 6.14104 | 6.291571 | 0.150531
$2d_{3/2}$ | -9.537 | -11.909 | 5.572499 | 5.577024 | 0.004525
$3s_{1/2}$ | -8.428 | -10.725 | 5.495602 | 5.482959 | -0.012643

In order to answer this question, we consider single-particle rms radii of the
proton orbitals occupied in 208Pb and 220Pb as obtained in the calculations
without pairing (see Table 1). These nuclei differ by 12 neutrons which
completely occupy the $\nu 1i_{11/2}$ subshell222We consider a full neutron
subshell in order to magnify the impact of the occupation of neutron states on
the proton single-particle rms radii. (see discussion in Ref. [15]). One can
see333Similar impact of the occupation of the $1i_{11/2}$ neutron subshell on
charge radii of proton subshells is seen also in the Skyrme DFT calculations
of Ref. [41]. that the pull provided by neutrons in this subshell on rms radii
of proton single-particle states depends on their principal quantum number
$n$. This pull is quantified by $\delta r_{p,i}$ (see Table 1) and it is the
largest for the $n=1$ orbitals (especially those located at the bottom of the
potential). It is reduced drastically for the $n=2$ and 3 proton orbitals.
Note that the proton rms point radius entering into Eq. (1) is given by

$\left<r^{2}\right>_{p}=\sum_{i}(2j_{i}+1)r_{p,i}^{2}$ (5)

where the sum runs over all occupied spherical subshells and $(2j_{i}+1)$ is
the degeneracy of the $i$-th subshell. As a result, the differential charge
radius between two nuclei (see Eq. (2)) can be rewritten as

$\delta\left<r^{2}\right>^{N,N^{\prime}}_{p}=\sum_{i}(2j_{i}+1)(r_{p,i}^{2}(N)-r_{p,i}^{2}(N^{\prime}))$
(6)

One can see that it depends on the pull on the rms radii of all occupied
proton subshells provided by the occupation of neutron states. Thus the
introduction of the core (as it is done in the spherical shell model) leads to
a neglect of the pull provided by the extra neutron(s) on the proton single-
particle states forming the core. This introduces uncontrollable errors in the
calculations of differential charge radii and thus severely limits the
applicability of spherical shell model to the description of this observable.

In a similar fashion, the lack of self-consistency effects is expected to
affect the description of differential charge radii in the mic+mac model. This
is because the addition of neutron does not affect the proton subsystem in a
self-consistent manner on the level of single-particle subshells (as in Table
1). For example, in the Woods-Saxon potential it affect the total radius of
nucleus $R$ only via mass dependence $R=1.2A^{1/3}$: this means that the
occupation of the neutron $1i_{11/2}$ and $2g_{9/2}$ subshells in the $N>126$
Pb isotopes will lead to the same differential charge radii contrary to the
results of self-consistent calculations (see Refs. [4, 5, 37, 24].) Moreover,
there is a lack of self-consistency in the definition of the radial properties
of the density distributions in the macroscopic (liquid drop) and microscopic
(single-particle potential) parts of the mic+mac model. To our knowledge, this
aspect of the problem has not been studied in details. However, the physical
observables similar to $\delta\left<r^{2}\right>_{p}^{N,N^{\prime}}$, namely,
relative charge quadrupole moments of superdeformed bands are affected by the
lack of self-consistency between microscopic and macroscopic parts [43].

##  8 Conclusions

Recent applications of covariant density functional theory to the analysis and
description of charge radii and related indicators have been reviewed. It
covers global accuracy of the description of experimental absolute and
differential charge radii and the evaluation of related theoretical
statistical and systematic uncertainties. This theory provides reasonably good
description of these physical observables already at the mean field level. The
discrepancies between theory and experiment have two sources: the deficiencies
in the description of the energies of the single-particle states and the
neglect of the correlations beyond mean field. For example, the former affects
the reproduction of the evolution of differential charge radii in isotopic
chains which cross the $N=28$ and $N=82$ neutron shell closures. The later
shows itself in two areas: (i) in the description of charge radii in shape
coexistent nuclei or in the nuclei with soft potential energy surfaces
[accurate description of such systems requires the methods such as generator
coordinate method or 5-dimensional collective Hamiltonian] and (ii) in the
description of odd-even staggering in charge radii since the energies and the
wave functions of one-quasiparticle states in odd-$A$ nuclei are affected
[strongly in spherical nuclei] by particle-vibration coupling.

In addition, new calculations reveal that differential charge radii in
deformed actinides and light superheavy nuclei are well described in the CDFT
calculations. However, it is important to take into account octupole
deformation in low-$N$ Th and U nuclei. Finally, we analyzed the importance of
self-consistency effects on the accurate description of differential charge.
This analysis reveals the challenging situation. The CDFT as well as the
Skyrme and Gogny DFTs take these effects into account on the level of single-
particle states but suffer from some deficiencies in the description of their
energies. In contrast, the models which are specifically fitted for a very
accurate description of the single-particle energies suffer from the neglect
of these self-consistency effects either due to introduction of the closed
core (as in the most of spherical shell models) or due to the split of the
model space into macroscopic and microscopic part (as in the mic+mac models).

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Nuclear Physics under Award No. DE-SC0013037. PR
acknowledges partial support from the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany Excellence Strategy
EXC-2094-390783311, ORIGINS.

## References

  * [1] I. Angeli and K. P. Marinova, At. Data Nucl. Data Tables 99, 69 (2013).
  * [2] P. Campbell, I. Moore, and M. Pearson, Prog. Part. Nucl. Phys. 86, 127 (2016).
  * [3] D. Vretenar, A. V. Afanasjev, G. A. Lalazissis, and P. Ring, Phys. Rep. 409, 101 (2005).
  * [4] M. M. Sharma, G. A. Lalazissis, and P. Ring, Phys. Lett. B 317, 9 (1993).
  * [5] P.-G. Reinhard and H. Flocard, Nucl. Phys. A 584, 467 (1995).
  * [6] Y. Tian, Z. Y. Ma, and P. Ring, Phys. Lett. B 676, 44 (2009).
  * [7] S. E. Agbemava and A. V. Afanasjev and D. Ray and P. Ring, Phys. Rev. C 89, 054320 (2014).
  * [8] G. A. Lalazissis, S. Karatzikos, R. Fossion, D. Peña Arteaga, A. V. Afanasjev and P. Ring, Phys. Lett. B671, 36 (2009).
  * [9] G. A. Lalazissis, T. Nikšić, D. Vretenar and P. Ring, Phys. Rev. C 71, 024312 (2005).
  * [10] X. Roca-Maza, X. Viñas, M. Centelles, P. Ring and P. Schuck, Phys. Rev. C 84, 055309 (2011).
  * [11] T. Nikšić, D. Vretenar and P. Ring, Phys. Rev. C 78, 034318 (2008).
  * [12] P. W. Zhao, Z. P. Li, J. M. Yao and J. Meng, Phys. Rev. C 82, 054319 (2010).
  * [13] S. E. Agbemava, A. V. Afanasjev, and A. Taninah, Phys. Rev. C 99, 014318 (2019).
  * [14] A. V. Afanasjev and S. E. Agbemava, Phys. Rev. C 93, 054310 (2016).
  * [15] U. C. Perera, A. V. Afanasjev, and P. Ring, subm. to Phys. Rev. C (2021), see also nuclear theory arXiv:2108.02245
  * [16] S. Frauendorf and V. V. Pashkevich, Phys. Lett. B 55, 365 (1975).
  * [17] J.-P. Delaroche, M. Girod, J. Libert, H. Goutte, S. Hilaire, S. Peru, N. Pillet, and G. F. Bertsch, Phys. Rev. C 81, 014303 (2010).
  * [18] E. W. Otten, Nuclear radii and moments of unstable isotopes, in Treatise on Heavy Ion Science: Volume 8: Nuclei Far From Stability, edited by D. A. Bromley (Springer, Boston, MA, (1989), pp. 517638.
  * [19] S. Sels et al, Phys. Rev. C 99, 044306 (2019).
  * [20] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * [21] J. Boguta and R. Bodmer, Nucl. Phys. A 292, 413 (1977).
  * [22] J. Dobaczewski, W. Nazarewicz, and P.-G. Reinhard, J. Phys. G 41, 074001 (2014).
  * [23] Y. Gao, J. Dobaczewski, M. Kortelainen, J. Toivanen, and D. Tarpanov, Phys. Rev. C 87, 034324 (2013).
  * [24] T. D. Goodacre et al, Phys. Rev. C, in press
  * [25] R. F. Garcia-Ruiz and A. Vernon, Eur. Phys. J. A 56, 136 (2020).
  * [26] S. E. Agbemava and A. V. Afanasjev and P. Ring, Phys. Rev. C 93, 044303 (2016).
  * [27] A. Voss, V. Sonnenschein, P. Campbell, B. Cheal, T. Kron, I. D. Moore, I. Pohjalainen, S. Raeder, N. Trautmann, and K. Wendt, Phys. Rev. A 95, 032506 (2017).
  * [28] S. Raeder et al., Phys. Rev. Lett. 120, 232503 (2018).
  * [29] Saleh O. Allehabi, V. A. Dzuba, V. V. Flambaum, A. V. Afanasjev and S. E. Agbemava, Phys. Rev. C 102, 024326 (2020).
  * [30] Yunchen Cao, S.E. Agbemava, A.V. Afanasjev, W. Nazarewicz and E. Olsen, Phys. Rev. C 102, 024311 (2020).
  * [31] A. V. Afanasjev and S. Shawaqfeh, Phys. Lett. B 706, 177 (2011).
  * [32] L. Bonneau, P. Quentin and P. Möller, Phys. Rev. C 76, 024320 (2007).
  * [33] J. Dobaczewski, A. V. Afanasjev, M. Bender, L. M. Robledo and Yue Shi, Nucl. Phys. A 944, 388 (2015).
  * [34] S. A. Fayans, S. V. Tolokonnikov, E. L. Trykov, and D. Zawischa, Phys. Lett. B 338, 1994 (1).
  * [35] S. A. Fayans, S. V. Tolokonnikov, E. L. Trykov, and D. Zawischa, Nucl. Phys. A 676, 49 (2000).
  * [36] P.-G. Reinhard and W. Nazarewicz, Phys. Rev. C 95, 064328 (2017).
  * [37] T. D. Goodacre et al, Phys. Rev. Let. 126, 032502 (2021).
  * [38] E. V. Litvinova and A. V. Afanasjev, Phys. Rev. C 84, 014305 (2011).
  * [39] A. V. Afanasjev and E. Litvinova, Phys. Rev. C 92, 044317 (2015).
  * [40] E. Litvinova and P. Ring, Phys. Rev. C 73, 044328 (2006).
  * [41] P. M. Goddard, P. D. Stevensson, and A. Rios, Phys. Rev. Lett. 110, 032503 (2013).
  * [42] G. Coló, H. Sagawa, and P. F. Bortignon, Phys. Rev. C 82, 064307 (2010).
  * [43] L. B. Karlsson, I. Ragnarsson and S. Åberg, Nucl. Phys. A 639, 654 (1998).


# DirectMultiStep: Direct Route Generation for Multi-Step Retrosynthesis

Yu Shee
Yale University
<EMAIL_ADDRESS>
Haote Li*
Yale University
<EMAIL_ADDRESS>
Anton Morgunov*
Yale University
<EMAIL_ADDRESS>
Victor Batista
Yale University
<EMAIL_ADDRESS>
Equal contribution. Listing order is random.

###### Abstract

Traditional computer-aided synthesis planning (CASP) methods rely on iterative
single-step predictions, leading to exponential search space growth that
limits efficiency and scalability. We introduce a transformer-based model that
directly generates multi-step synthetic routes as a single string by
conditionally predicting each molecule based on all preceding ones. The model
accommodates specific conditions such as the desired number of steps and
starting materials, outperforming state-of-the-art methods on the PaRoutes
dataset with a 2.2x improvement in Top-1 accuracy on the n1 test set and a
3.3x improvement on the n5 test set. It also successfully predicts routes for
FDA-approved drugs not included in the training data, showcasing its
generalization capabilities. While the current suboptimal diversity of the
training set may impact performance on less common reaction types, our
approach presents a promising direction towards fully automated retrosynthetic
planning.

##  1 Introduction

Finding the most efficient route to a desired chemical compound is a daily
challenge for synthetic organic chemists in both fundamental research and drug
discovery. Route efficiency is determined by various factors, some of which
can be objectively assessed, such as overall yield (not all chemical reactions
have 100% conversion rate) and enantiomeric excess (in case of chiral
compounds), where higher values are always preferred. Other factors, such as
atom efficiency (minimization of byproducts) and availability (cost) of
starting materials, are more case-dependent. A used chemical reactant is
considered waste unless it can be repurposed as a reactant in a different
process. Similarly, the choice of starting materials depends on factors such
as budget, logistics, and the availability of specific equipment. It’s worth
noting that many commercially available compounds can be synthesized from
other commercially available compounds, adding another layer of complexity to
the decision-making process.

A significant contribution to the development of algorithmic frameworks for
identifying synthetic routes leading to any desired compound (subject of the
1990 Nobel Prize in Chemistry) was made by Elias James Corey. Corey’s
framework, now known as retrosynthetic analysis, begins with identifying atoms
that would serve as reaction centers. Disconnecting bonds between these
centers results in the formation of hypothetical fragments (called synthons)
from which a precursor molecule can be created. This mapping from fragments to
actual molecules is one-to-many because there is usually more than one
functional group that could participate in any given type of reaction.
Importantly, meticulous application of Corey’s framework (i.e., systematically
breaking small subsets of bonds) will eventually lead to commercially
available starting materials. The algorithmic nature of this process allowed
Corey to envision automating these rules to create Computer-Aided Synthesis
Planning (CASP) as early as 1969 [1].

Recent advancements in data science and machine learning (ML) methods have led
to a surge of interest in developing CASP methods [2, 3]. The vast majority of
existing methods [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35] are designed
to automate single-step retrosynthetic (SSR) analysis, i.e. predicting a list
of compounds from which a target product could be made in one step. A full
multi-step route can be determined by iteratively applying SSR methods to the
identified precursors until a termination condition is satisfied (e.g.,
identifying reactants in the database of commercially available compounds).
Notably, because each SSR method call creates a list of candidates, iterative
application of these methods generates an exponentially growing search space.
Graph traversal algorithms such as Monte Carlo Tree Search (MCTS), Depth-First
Proof Number (DFPN) search, and A-star (Retro*) have been applied to
efficiently traverse this exponential search space [36, 37, 38].

However, evaluation of performance of these SSR methods on full route
prediction has been limited by the scarcity of open-source datasets containing
valid multi-step routes. The creation of PaRoutes [39], an open-source dataset
containing 450k multi-step routes extracted from the United States Patent and
Trademark Office (USPTO), marked a major development. Notably, even state-of-
the-art SSR models [6, 40] combined with advanced search algorithms (MCTS,
DFPN, Retro*) correctly identify multi-step routes (Top-1 accuracy) for only
17% and 10% of the target compounds in n1 and n5, respectively. This
performance can be rationalized by recognizing that retrosynthesis is
inherently a multi-step problem: the optimal choice of reaction to make a
compound depends on subsequent steps in the synthesis. For example, a common
pattern in multi-step routes includes: 1) protection of certain functional
groups, 2) the desired transformation of unprotected groups, and 3)
deprotection of the protected groups. An SSR method applied to the reactant of
reaction 2 may output numerous candidate precursors; however, knowing that
protective groups are removed in subsequent steps dramatically changes the
probability distribution over those precursors.

In this work, we propose a novel approach for direct prediction of multi-step
routes, bypassing the need for single-step models and sophisticated
exponential graph traversal algorithms. Our DirectMultiStep model demonstrates
state-of-the-art performance on n1 and n5 evaluation sets and shows
generalizability by successfully predicting routes for FDA-approved drugs
absent from the training set. We also discuss the model’s limitations, which
may include the need to split longer routes into sub-routes, provide the
structure of the starting material, or supply a precursor to the target
compound when dealing with transformations not well-represented in the
training dataset. Despite these limitations, we believe our results are
promising enough to warrant further attention to developing retrosynthetic
methods that directly predict multi-step routes.

##  2 Methods

###  2.1 Definitions and multi-step route representation

This subsection defines key terms related to our multi-step retrosynthetic
methodology. The "target compound" (green in Fig. 1) represents the final
product of the multi-step synthesis tree. Starting materials (SM, orange in
Fig. 1) are compounds for which no further precursors need to be identified,
i.e., they are the leaves of the synthesis tree. The "number of steps" denotes
the largest number of reactions from SM to the target compound (tree height).
In our approach, routes are represented as recursive dictionaries (Fig. 1)
containing the SMILES representation of the molecule and a list of other
dictionaries, containing either starting materials or trees leading to
precursors of the current node. Removing space and newline characters creates
a string representation predicted by our transformer model. We train the
DirectMultiStep (DMS) transformer, which takes the target compound,
(optionally) starting material, and number of steps as input, and predicts the
multi-step route represented as a string.

###  2.2 Data curation and preparation

A second version of the PaRoutes [39] dataset containing 450k routes is
stripped of all metadata and stored as recursive dictionaries representing
multi-step routes. The two evaluation sets, n1 and n5, are processed
similarly. A training partition is created by removing all permutations
(swapping left and right subtrees) of routes in n1 and n5 from the full
dataset, resulting in 407k routes. The training dataset is augmented by adding
2 permutations for each route. To train the DMS model which takes starting
material as an input, we find all starting materials (leaves) for each tree
and store a combination of the target compound with each starting material as
a separate entry. As a result, DMS with SM is trained on 3 349 118 inputs, and
the DMS without SM on 1 078 951 inputs. All SMILES strings are tokenized by
treating each character as a token, and the string representation of the
multi-step route is tokenized similarly, treating delimiters of the tree
(`smiles’`, `children’`, `[’`, `]’`, `{’`, and `}’`) as separate tokens. The
final vocabulary size is 52 (including start, end, and padding tokens), the
largest multi-step route has 1074 tokens, while the largest target compound
and starting material have 145 and 135 tokens, respectively.

Figure 1: The workflow of DirectMultiStep. (a) The SMILES representation of
the target compound (green), (optional) starting material (orange), and the
number of steps are tokenized, concatenated, and fed into our transformer
model. The model predicts a string representation of the multi-step synthesis
tree. Spaces are added for clarity, and indentations indicate the levels in
the synthesis route (tree). (b) Molecular structures corresponding to the
target compound (green), (optional) starting material (orange), and the
predicted synthesis tree with structures of all molecules. Figure 2:
Distribution of the relative frequencies of route lengths (in terms of number
of steps) in the training dataset before augmentation with permutations (407
148 routes, black), n1 test set (10 000 routes, blue), and n5 test set (10 000
routes, purple). Distribution is split into routes shorter (left subplot) and
longer than 6 steps (right subplot).

###  2.3 Model architecture and training

We present two variants of the DirectMultiStep (DMS) model: DMS-10M and
DMS-60M, with approximately 10 million and 60 million parameters,
respectively. Both models employ the transformer encoder-decoder architecture
[41]. DMS-10M has 6 layers and 8 attention heads in both encoder and decoder,
with inputs embedded in 256 dimensions and amplified 3x in the feedforward
blocks. DMS-60M has 8 layers and 8 attention heads, with inputs embedded in
512 dimensions and amplified 4x in the feedforward blocks. Both models use
GeLU activations [42]. The encoder of DMS with SM receives a concatenation of
the tokenized target compound, starting material, and number of steps, while
the encoder of DMS without SM receives only the tokenized target compound and
number of steps. In all models, the decoder receives output from the encoder
and the start token of the multi-step route.

During training, a 10% dropout rate is used in embedding, attention, and
feedforward layers. The ADAM optimizer [43] is used with a learning rate
scheduler that warms up from 0 to $3\times 10^{-4}$ for DMS-10M and $1\times
10^{-4}$ for DMS-60M during the first 10% of training steps and undergoes
cosine decay to $3\times 10^{-5}$ over 20 epochs. Gradients are clipped to
1.0. Finally, during training, each token in the input to the encoders has a
5% chance of being masked. Training requires approximately 200 clock hours
using 4 NVIDIA A100-40 GB GPUs on an external cluster. An additional 450 hours
were used for experimentation with model architectures and optimization of
hyperparameters.

###  2.4 Beam search for route generation

We utilize beam search with a width of 50 to generate several candidates for
the same target compound (and starting material, if provided). The cumulative
probability of each beam serves as the ranking metric for prediction accuracy
(Top-$K$). After beam search, predicted routes are checked for the validity of
all SMILES strings (which are canonicalized), the correctness of the target
compound, the presence of starting materials in the stock compound (if
applicable), and the absence of repetitions.

###  2.5 Evaluation metrics and post-processing

Evaluation is performed on all 10 000 routes in both n1 and n5, with the
deepest leaf selected as the starting material to be included in the input to
the encoder (for models that use starting materials). Top-$K$ accuracy
evaluates the percentage of top $K$ predictions precisely matching correct
multi-step routes from the test set. Predictions matching any permutation of
the correct route are counted as matches.

##  3 Results and Discussion

###  3.1 Test accuracy

Top-$K$ accuracy is shown in Table 1 for set-n1 and Table 2 for set-n5. Our
DirectMultiStep-SM-10M model achieves state-of-the-art performance, surpassing
all previous methods in all Top-$K$ accuracies, with the most significant
increase in Top-1 accuracy: 2.2x on n1 and 3.3x on n5. Remarkably, our model’s
Top-3 accuracy exceeds the Top-10 accuracy of top tree search methods using
single-step retrosynthesis (SSR), highlighting the advantages of the multi-
step first approach.

The DirectMultiStep-SM-60M model, despite having six times more parameters
than the 10M variant, exhibits lower performance on higher Top-$K$ accuracies,
suggesting that the 10M model size is optimal for the current dataset when
starting material structure is provided. However, achieving comparable
performance without starting material information requires a larger model, as
demonstrated by the DirectMultiStep-noSM-60M model’s performance. These
findings reveal a trade-off between model generality and size at the current
dataset scale: less input information necessitates a larger model to maintain
high performance. In all subsequent analysis, we use the routes predicted by
the DirectMultiStep-SM-10M model.

Table 1: Top-$K$ accuracy on route test set-n1 (10 000 routes).

Method | Top-1 | Top-2 | Top-3 | Top-4 | Top-5 | Top-10
---|---|---|---|---|---|---
MCTSa | 0.17 |  |  |  | 0.46 | 0.49
Retro*a | 0.15 |  |  |  | 0.41 | 0.45
DFPNa | 0.11 |  |  |  | 0.17 | 0.17
DirectMultiStep-SM-10M | 0.38 | 0.46 | 0.50 | 0.52 | 0.53 | 0.56
DirectMultiStep-SM-60M | 0.38 | 0.46 | 0.48 | 0.50 | 0.51 | 0.52
DirectMultiStep-noSM-60M | 0.36 | 0.42 | 0.44 | 0.45 | 0.45 | 0.46

aThese data are collected from the 2.0 version of PaRoutes in their GitHub
repository [39].

Table 2: Top-$K$ accuracy on route test set-n5 (10 000 routes).

Method | Top-1 | Top-2 | Top-3 | Top-4 | Top-5 | Top-10
---|---|---|---|---|---|---
MCTSa | 0.10 |  |  |  | 0.28 | 0.33
Retro*a | 0.10 |  |  |  | 0.30 | 0.36
DFPNa | 0.05 |  |  |  | 0.07 | 0.07
DirectMultiStep-SM-10M | 0.33 | 0.40 | 0.43 | 0.44 | 0.46 | 0.48
DirectMultiStep-SM-60M | 0.34 | 0.41 | 0.43 | 0.44 | 0.45 | 0.46
DirectMultiStep-noSM-60M | 0.32 | 0.37 | 0.39 | 0.40 | 0.41 | 0.41

aThese data are collected from the 2.0 version of PaRoutes in their GitHub
repository [39]

Importantly, our model requires only one call to obtain multiple routes for
each target compound, instead of hundreds of single-step model calls and even
more reaction template applications (as shown in PaRoutes [39]). Furthermore,
our model’s inference time (beam size 50) on a single GPU is comparable (3-16
seconds depending on the number of steps) to the time required to find the
first successful route using previous methods (7-50 seconds [39])

Fig. 3 shows the distribution of Top-1 and Top-10 accuracy over different
route lengths on both n1 and n5. Given that 90% of routes in the training
partition have 4 or fewer steps Fig. 2), one would expect accuracy to decrease
dramatically with increasing route length. However, the performance on routes
with 5-8 steps is comparable to that of shorter routes, and the performance on
9-step routes is even comparable to that of 2-step routes. The exceptional
accuracy on 9-step routes could be rationalized by the fact that their
relative abundance in n1 and n5 does not exceed the relative abundance in the
training partition (unlike for routes with length 3-7, as seen in Fig. 3).

Figure 3: Distribution of Top-1 and Top-10 accuracy of predictions with our
model on test sets n1 and n5. There is only one route with length 10 in n1,
and our model does not predict it correctly. That route is reproduced by
splitting it in half, as shown in Fig. 4

Figure 4: Separation of a 10-step route from set-n1 into two 5-step routes.
(a) Correct prediction for the first half of the 10-step route with starting
material information. (b) Correct prediction for the second half of the
10-step route with starting material information.

The requirement to provide the starting material structure as input to the
model is a limitation, although not as significant as it may seem. Corey’s
retrosynthetic framework allows finding at least one route to starting
materials by performing one functional group transformation at a time, usually
requiring protection of other functional groups. However, such a route will be
too long to be efficient in terms of yield (overall yield is the product of
yields of individual steps). The art of organic synthesis lies in finding ways
to perform several transformations simultaneously using a minimal number of
protecting groups. We believe our model provides value as an assistant to this
task: finding better ways to get to the target compound from a starting
material. The requirement to specify the number of steps as an input becomes a
useful option in this context. Moreover, if the target compound is chiral
(i.e., it has a mirror image which cannot be superimposed), even during manual
retrosynthetic planning, chemists check if chiral centers can be obtained from
chiral starting materials (such as naturally occurring amino acids) because
developing reactions that create chiral centers in the desired conformation
remains an active research area (and subject of 2001 and 2021 Nobel Prizes in
Chemistry). An example of such a route is shown in Fig. 6.

###  3.2 Separation of Long Routes into Shorter Routes

As seen in Fig. 3, our model does not find the single 10-step route in the n1
set. However, this route can be correctly reproduced with rank 1 if split into
two 5-step routes and the intermediate compound is provided as the starting
material for the first 5 reactions (Fig. 4) and a target compound for the last
5 reactions (Fig. 4). While ideally the necessity to provide an intermediate
compound should be avoided, this resembles the bidirectional search technique
commonly employed in manual retrosynthetic analysis [44].

###  3.3 Retrosynthetic planning of pharmaceutical compounds

We demonstrate the generalizability of our model by testing on three FDA-
approved drugs: Vonoprazan, Mitapivat, and Daridorexant, which were used for
evaluation by Xiong, et al. [45]. These drugs and their intermediates are
absent from the training set. Vonoprazan, a potassium competitive acid blocker
for Helicobacter pylori infections, was initially proposed with a 2-step route
(Fig. 5) [46]. However, one of its starting materials is unstable, and the
final reaction leads to significant byproducts. Subsequently, an improved
4-step route (Fig. 5) was introduced [47]. Mitapivat, a pyruvate kinase
activator for treating hemolytic anemia, has been associated with both a
5-step route (Fig. 5) and a 3-step route (Fig. 5) [48, 49]. Daridorexant, an
orexin receptor antagonist for adult insomnia, is linked to a 4-step route
(Fig. 6) [50].

Figure 5: Literature routes for Vonoprazan and Mitapivat correctly reproduced
by our model. Target compounds are in green, starting materials that are given
as inputs are colored in orange. Ranks denote the rank of this route when the
specified starting material is provided. (a) First literature route for
Vonoprazan. The model predicts correctly no matter which starting material is
given. (b) Second literature route for Vonoprazan. The model predicts the
route correctly only when an immediate precursor to Vonoprazan is given as the
target compound. (c) First literature route for Mitapivat. (d) Second
literature route for Mitapivat.

Our model correctly reproduces Vonoprazan’s first literature route (Fig. 5)
regardless of the starting material provided, with high ranks. However, the
model struggles to find the first transformation (counting from the target
compound) from the second route. In the first route, the first step is a
reductive amination, while in the second route, it’s a reduction of the
carbonyl bond in the amide group. Reductive amination is much more common for
creating C-N bonds, so the training set is heavily biased toward it. As a
result, all routes predicted by our model start as in Fig. 5. However, if we
provide the precursor to Vonoprazan in the second route as the target
compound, the model correctly reproduces the full route (Fig. 5). This
highlights an important limitation of our model: it may not perform as well on
routes involving transformations that are inadequately represented in the
training dataset. We envision that as larger multi-step datasets become
available, this problem will naturally disappear.

Figure 6: Daridorexant literature multi-step routes and predictions from the
model with starting material information. The compounds that are missing from
the prediction are in red.

Our model correctly reproduces both routes for Mitapivat (Fig. 5 and Fig. 5).
In contrast, the route for Daridorexant (Fig. 6) is only partially predicted:
our model does not split the right child of the root node (Daridorexant
itself) into the two compounds highlighted in red. All other intermediates and
starting materials are predicted correctly. The reason for this behavior is
not well understood, but it can be rationalized by looking at the distribution
of the number of leaves at root nodes: 73% of root nodes (target compounds)
have at least 1 leaf (see SI Fig. 7). Such dependence on the distributions of
the training set comprises an important limitation of the multi-step first
approach. This issue would not occur with methods based on single-step
prediction because the tree search framework would not halt at compounds
outside the stock compound set before reaching the maximum iterations.
Ensuring that multi-step first models reproduce that aspect of SSR methods
constitutes a direction for future research.

##  4 Conclusion

In this work, we introduced DirectMultiStep, a transformer-based model that
directly generates multi-step synthetic routes for a given target compound. By
predicting the entire route as a single string, our approach bypasses the need
for iterative single-step predictions and exponential search space traversal.
The model accommodates specific conditions such as the desired number of steps
and starting materials, enabling efficient route planning.

DirectMultiStep outperforms state-of-the-art methods on the PaRoutes dataset,
achieving a 2.2x improvement in Top-1 accuracy on the n1 test set and a 3.3x
improvement on the n5 test set. The model also successfully predicts routes
for FDA-approved drugs not included in the training data, showcasing its
generalization capabilities. However, the model’s performance may be impacted
by the suboptimal diversity of the current training set, particularly for less
common reaction types.

Despite these limitations, the DirectMultiStep approach presents a promising
direction towards fully automated retrosynthetic planning. As larger and more
diverse multi-step datasets become available, the model’s performance is
expected to improve further. Future research should focus on addressing the
limitations identified, such as the need to split longer routes into sub-
routes and the dependence on the distributions of the training set. Overall,
this work demonstrates the potential of transformer-based models for efficient
and scalable computer-aided synthesis planning.

##  5 Broad Impact

Chemical synthesis of target compounds of interest forms the routine for
millions of organic chemists around the world, including graduate students in
PhD programs and synthetic chemists in pharmaceutical companies. The ability
to find shorter (and thus more efficient) routes to any desired compound is an
enticing challenge with the potential of saving on labor costs and chemical
waste. As a result, multi-step prediction methods have the potential to
improve the lives of millions of people, and potentially contribute to faster
advances in the development of new drugs and materials. We hope that our work
constitutes an important advancement in that direction.

Quite unfortunately, however, there is no fundamental chemical difference
between how safe and useful compounds (such as drugs) and malign compounds
(such as chemical weapons) are synthesized. Both require application of the
same set of rules of organic synthesis, and so the ability to make drugs
faster and more efficiently inherently brings along the ability to make
poisons faster and more efficiently. One cannot impose restrictions on the
tokens (chemical elements) learned by the model, as even some combinations of
nitrogen and oxygen (which are widespread in proteins and DNA) may be
cytotoxic. Therefore, unless the humanity at large decides that the benefits
of automated retrosynthesis do not outweigh its potential harms, the most we
can do is communicate the dangers clearly and perhaps impose restrictions that
would still require a trained (which often accompanies training in ethics)
chemist to use our model. Besides, all meta-data of reaction conditions and
catalysts are removed from the predictions of our model, which would present a
deliberate barrier against malicious uses.

##  6 Code and Data Availability

The reaction route dataset is from the 2.0 version of PaRoutes in their GitHub
repository [39], available under the Apache License 2.0. Code to process
dataset, implementation of model architecture, code for training, generation,
and evaluation are available under MIT License at
https://github.com/batistagroup/DirectMultiStep.

##  7 Acknowledgments

The authors acknowledge a generous allocation of high-performance computing
time from NERSC. The development of the methodology was supported by the NSF
CCI grant (VSB, Award Number 2124511).

## References

  * [1] E. J. Corey and W. Todd Wipke.  Computer-Assisted Design of Complex Organic Syntheses: Pathways for molecular synthesis can be devised with a computer and equipment for graphical communication.  Science, 166(3902):178–192, 1969.
  * [2] A. Filipa De Almeida, Rui Moreira, and Tiago Rodrigues.  Synthetic organic chemistry driven by artificial intelligence.  Nat Rev Chem, 3(10):589–604, 2019.
  * [3] Thomas J. Struble, Juan C. Alvarez, Scott P. Brown, Milan Chytil, Justin Cisar, Renee L. DesJarlais, Ola Engkvist, Scott A. Frank, Daniel R. Greve, Daniel J. Griffin, Xinjun Hou, Jeffrey W. Johannes, Constantine Kreatsoulas, Brian Lahue, Miriam Mathea, Georg Mogk, Christos A. Nicolaou, Andrew D. Palmer, Daniel J. Price, Richard I. Robinson, Sebastian Salentin, Li Xing, Tommi Jaakkola, William. H. Green, Regina Barzilay, Connor W. Coley, and Klavs F. Jensen.  Current and Future Roles of Artificial Intelligence in Medicinal Chemistry Synthesis.  J. Med. Chem., 63(16):8667–8682, 2020.
  * [4] Zhongliang Guo, Stephen Wu, Mitsuru Ohno, and Ryo Yoshida.  Bayesian Algorithm for Retrosynthesis.  J. Chem. Inf. Model., 60(10):4474–4486, October 2020.
  * [5] Hankook Lee, Sungsoo Ahn, Seung-Woo Seo, You Young Song, Eunho Yang, Sung-Ju Hwang, and Jinwoo Shin.  Retcl: A selection-based approach for retrosynthesis via contrastive learning, 2021.
  * [6] Marwin H. S. Segler and Mark P. Waller.  Neural-Symbolic Machine Learning for Retrosynthesis and Reaction Prediction.  Chem. Eur. J., 23(25):5966–5971, May 2017.
  * [7] Connor W. Coley, Luke Rogers, William H. Green, and Klavs F. Jensen.  Computer-Assisted Retrosynthesis Based on Molecular Similarity.  ACS Cent. Sci., 3(12):1237–1245, December 2017.
  * [8] Shoichi Ishida, Kei Terayama, Ryosuke Kojima, Kiyosei Takasu, and Yasushi Okuno.  Prediction and Interpretable Visualization of Retrosynthetic Reactions Using Graph Convolutional Networks.  J. Chem. Inf. Model., 59(12):5026–5033, December 2019.
  * [9] Michael E. Fortunato, Connor W. Coley, Brian C. Barnes, and Klavs F. Jensen.  Data Augmentation and Pretraining for Template-Based Retrosynthetic Prediction in Computer-Aided Synthesis Planning.  J. Chem. Inf. Model., 60(7):3398–3407, July 2020.
  * [10] Hanjun Dai, Chengtao Li, Connor W. Coley, Bo Dai, and Le Song.  Retrosynthesis prediction with conditional graph logic network, 2020.
  * [11] Shuan Chen and Yousung Jung.  Deep Retrosynthetic Reaction Prediction using Local Reactivity and Global Attention.  JACS Au, 1(10):1612–1620, October 2021.
  * [12] Philipp Seidl, Philipp Renz, Natalia Dyubankova, Paulo Neves, Jonas Verhoeven, Jörg K. Wegner, Marwin Segler, Sepp Hochreiter, and Günter Klambauer.  Improving Few- and Zero-Shot Reaction Template Prediction Using Modern Hopfield Networks.  J. Chem. Inf. Model., 62(9):2111–2120, May 2022.
  * [13] Chaochao Yan, Qianggang Ding, Peilin Zhao, Shuangjia Zheng, Jinyu Yang, Yang Yu, and Junzhou Huang.  Retroxpert: Decompose retrosynthesis prediction like a chemist, 2020.
  * [14] Chence Shi, Minkai Xu, Hongyu Guo, Ming Zhang, and Jian Tang.  A graph to graphs framework for retrosynthesis prediction, 2021.
  * [15] Vignesh Ram Somnath, Charlotte Bunne, Connor W. Coley, Andreas Krause, and Regina Barzilay.  Learning graph models for retrosynthesis prediction, 2021.
  * [16] Xiaorui Wang, Yuquan Li, Jiezhong Qiu, Guangyong Chen, Huanxiang Liu, Benben Liao, Chang-Yu Hsieh, and Xiaojun Yao.  RetroPrime: A Diverse, plausible and Transformer-based method for Single-Step retrosynthesis predictions.  Chemical Engineering Journal, 420:129845, September 2021.
  * [17] Yu Wang, Chao Pang, Yuzhe Wang, Junru Jin, Jingjie Zhang, Xiangxiang Zeng, Ran Su, Quan Zou, and Leyi Wei.  Retrosynthesis prediction with an interpretable deep-learning framework based on molecular assembly tasks.  Nat Commun, 14(1):6155, October 2023.
  * [18] Weihe Zhong, Ziduo Yang, and Calvin Yu-Chian Chen.  Retrosynthesis prediction using an end-to-end graph generative architecture for molecular graph editing.  Nat Commun, 14(1):3009, May 2023.
  * [19] Bowen Liu, Bharath Ramsundar, Prasad Kawthekar, Jade Shi, Joseph Gomes, Quang Luu Nguyen, Stephen Ho, Jack Sloane, Paul Wender, and Vijay Pande.  Retrosynthetic Reaction Prediction Using Neural Sequence-to-Sequence Models.  ACS Cent. Sci., 3(10):1103–1113, October 2017.
  * [20] Pavel Karpov, Guillaume Godin, and Igor Tetko.  A Transformer Model for Retrosynthesis.  preprint, Chemistry, May 2019.
  * [21] Benson Chen, Tianxiao Shen, Tommi S. Jaakkola, and Regina Barzilay.  Learning to make generalizable and diverse predictions for retrosynthesis, 2019.
  * [22] Alpha A. Lee, Qingyi Yang, Vishnu Sresht, Peter Bolgar, Xinjun Hou, Jacquelyn L. Klug-McLeod, and Christopher R. Butler.  Molecular Transformer unifies reaction prediction and retrosynthesis across pharma chemical space.  Chem. Commun., 55(81):12152–12155, 2019.
  * [23] Kangjie Lin, Youjun Xu, Jianfeng Pei, and Luhua Lai.  Automatic retrosynthetic route planning using template-free models.  Chem. Sci., 11(12):3355–3364, 2020.
  * [24] Shuangjia Zheng, Jiahua Rao, Zhongyue Zhang, Jun Xu, and Yuedong Yang.  Predicting Retrosynthetic Reactions Using Self-Corrected Transformer Neural Networks.  J. Chem. Inf. Model., 60(1):47–55, January 2020.
  * [25] Igor V. Tetko, Pavel Karpov, Ruud Van Deursen, and Guillaume Godin.  State-of-the-art augmented NLP transformer models for direct and single-step retrosynthesis.  Nat Commun, 11(1):5575, November 2020.
  * [26] Seung-Woo Seo, You Young Song, June Yong Yang, Seohui Bae, Hankook Lee, Jinwoo Shin, Sung Ju Hwang, and Eunho Yang.  GTA: Graph Truncated Attention for Retrosynthesis.  AAAI, 35(1):531–539, May 2021.
  * [27] Kelong Mao, Xi Xiao, Tingyang Xu, Yu Rong, Junzhou Huang, and Peilin Zhao.  Molecular graph enhanced transformer for retrosynthesis prediction.  Neurocomputing, 457:193–202, October 2021.
  * [28] Mikołaj Sacha, Mikołaj Błaż, Piotr Byrski, Paweł Dąbrowski-Tumański, Mikołaj Chromiński, Rafał Loska, Paweł Włodarczyk-Pruszyński, and Stanisław Jastrzębski.  Molecule Edit Graph Attention Network: Modeling Chemical Reactions as Sequences of Graph Edits.  J. Chem. Inf. Model., 61(7):3273–3284, July 2021.
  * [29] Vipul Mann and Venkat Venkatasubramanian.  Retrosynthesis prediction using grammar-based neural machine translation: An information-theoretic approach.  Computers & Chemical Engineering, 155:107533, December 2021.
  * [30] Umit V. Ucak, Taek Kang, Junsu Ko, and Juyong Lee.  Substructure-based neural machine translation for retrosynthetic prediction.  J Cheminform, 13(1):4, December 2021.
  * [31] Eunji Kim, Dongseon Lee, Youngchun Kwon, Min Sik Park, and Youn-Suk Choi.  Valid, Plausible, and Diverse Retrosynthesis Using Tied Two-Way Transformers with Latent Variables.  J. Chem. Inf. Model., 61(1):123–133, January 2021.
  * [32] Ross Irwin, Spyridon Dimitriadis, Jiazhen He, and Esben Jannik Bjerrum.  Chemformer: a pre-trained transformer for computational chemistry.  Mach. Learn.: Sci. Technol., 3(1):015022, March 2022.
  * [33] Zipeng Zhong, Jie Song, Zunlei Feng, Tiantao Liu, Lingxiang Jia, Shaolun Yao, Min Wu, Tingjun Hou, and Mingli Song.  Root-aligned SMILES: a tight representation for chemical reaction prediction.  Chem. Sci., 13(31):9023–9034, 2022.
  * [34] Umit V. Ucak, Islambek Ashyrmamatov, Junsu Ko, and Juyong Lee.  Retrosynthetic reaction pathway prediction through neural machine translation of atomic environments.  Nat Commun, 13(1):1186, March 2022.
  * [35] Yu Shee, Haote Li, Pengpeng Zhang, Andrea Nikolic, Sanil Sreekumar, Frédéric Buono, Jinhua Song, Timothy Newhouse, and Victor Batista.  Site-specific template generative approach for retrosynthetic planning.  ChemRxiv, 2024.
  * [36] Marwin H. S. Segler, Mike Preuss, and Mark P. Waller.  Planning chemical syntheses with deep neural networks and symbolic AI.  Nature, 555(7698):604–610, 2018.
  * [37] Akihiro Kishimoto, Beat Buesser, Bei Chen, and Adi Botea.  Depth-first proof-number search with heuristic edge cost and application to chemical synthesis planning.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [38] Binghong Chen, Chengtao Li, Hanjun Dai, and Le Song.  Retro*: Learning retrosynthetic planning with neural guided A* search.  In The 37th International Conference on Machine Learning (ICML 2020), 2020.
  * [39] Samuel Genheden and Esben Bjerrum.  PaRoutes: towards a framework for benchmarking retrosynthesis route predictions.  Digital Discovery, 1(4):527–539, 2022.
  * [40] Amol Thakkar, Thierry Kogej, Jean-Louis Reymond, Ola Engkvist, and Esben Jannik Bjerrum.  Datasets and their influence on the development of computer assisted synthesis planning tools in the pharmaceutical domain.  Chem. Sci., 11(1):154–168, 2020.
  * [41] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [42] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus), 2023.
  * [43] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * [44] E. J. Corey and Xue-Min Cheng.  Logic of chemical synthesis.  Wiley, New York, new ed. edition, 1995.
  * [45] Jiacheng Xiong, Wei Zhang, Zunyun Fu, Jiatao Huang, Xiangtai Kong, Yitian Wang, Zhaoping Xiong, and Mingyue Zheng.  Improve retrosynthesis planning with a molecular editing language.  preprint, Chemistry, 2023.
  * [46] Masahiro Kajino, Atsushi Hasuoka, and Haruyuki Nishida.  1-heterocyclylsulfonyl, 2 aminomethyl, 5- (hetero-) aryl substituted 1-h-pyrrole derivatives as acid secretion inhibitors., WO2007026916A1, 2007.
  * [47] Qian-Ying Yu, Huang Zeng, Kai Yao, Jian-Qi Li, and Yu Liu.  Novel and practical synthesis of vonoprazan fumarate.  Synthetic Communications, 47(12):1169–1174, 2017.
  * [48] Jeffrey O. Saunders, Francesco G. Salituro, and Shunqi Yan.  Therapeutic compounds and compositions, WO2011002817A1, 2011.
  * [49] Jacob P. SIZEMORE, Liting GUO, Mahmoud Mirmehrabi, and Yeqing Su.  Crystalline forms of n-(4-(4-(cyclopropylmethyl) piperazine-1-carbonyl)phenyl)quinoline-8-sulfonamide, WO2019104134A1, 2019.
  * [50] Zhefeng Zhang, Aiqin Zhang, and Mouli Yang.  Fused-ring compound having analgesic activity, and preparation method therefor and use thereof, WO2023160004A1, 2023.

##  8 Supplementary Information

Table 3: Top-$K$ accuracy on subset of test set-n1 (random 500 routes). Method | Top-1 | Top-2 | Top-3 | Top-4 | Top-5 | Top-10
---|---|---|---|---|---|---
DirectMultiStep-SM-10M | 0.402 | 0.476 | 0.522 | 0.550 | 0.564 | 0.582
DirectMultiStep-SM-60M | 0.398 | 0.462 | 0.496 | 0.518 | 0.526 | 0.540
DirectMultiStep-noSM-10M | 0.292 | 0.362 | 0.384 | 0.388 | 0.396 | 0.402
DirectMultiStep-noSM-60M | 0.366 | 0.438 | 0.460 | 0.478 | 0.484 | 0.490
Table 4: Top-$K$ accuracy subset of test set-n5 (random 500 routes). Method | Top-1 | Top-2 | Top-3 | Top-4 | Top-5 | Top-10
---|---|---|---|---|---|---
DirectMultiStep-SM-10M | 0.374 | 0.442 | 0.474 | 0.484 | 0.498 | 0.516
DirectMultiStep-SM-60M | 0.344 | 0.428 | 0.454 | 0.456 | 0.462 | 0.480
DirectMultiStep-noSM-10M | 0.222 | 0.264 | 0.282 | 0.288 | 0.290 | 0.298
DirectMultiStep-noSM-60M | 0.310 | 0.368 | 0.388 | 0.400 | 0.406 | 0.418

Table 3 and Table 4 present the Top-$K$ accuracy results for 500 random routes
each from set-n1 and set-n5. These results illustrate that a larger model is
necessary to maintain high performance when starting material information is
not provided, as shown by the DirectMultiStep-noSM-60M model’s superior
performance compared to DirectMultiStep-noSM-10M.

Figure 7: Distribution of the relative frequencies for the number of leaves
(nodes with no children) at root node.