# Transfer learning for affordable and high quality tunneling splittings from
instanton calculations

Silvan Käser  [  Jeremy O. Richardson  [  Markus Meuwly  [<EMAIL_ADDRESS>
###### Abstract

The combination of transfer learning (TL) a low level potential energy surface
(PES) to a higher level of electronic structure theory together with ring-
polymer instanton (RPI) theory is explored and applied to malonaldehyde. The
RPI approach provides a semiclassical approximation of the tunneling splitting
and depends sensitively on the accuracy of the PES. With second order Møller-
Plesset perturbation theory (MP2) as the low-level (LL) model and energies and
forces from coupled cluster singles, doubles and perturbative triples
(CCSD(T)) as the high-level (HL) model, it is demonstrated that CCSD(T)
information from only 25 to 50 judiciously selected structures along and
around the instanton path suffice to reach HL-accuracy for the tunneling
splitting. In addition, the global quality of the HL-PES is demonstrated
through a mean average error of 0.3 kcal/mol for energies up to 40 kcal/mol
above the minimum energy structure (a factor of 2 higher than the energies
employed during TL) and $<2$ cm-1 for harmonic frequencies compared with
computationally challenging normal mode calculations at the CCSD(T) level.

University of Basel]Department of Chemistry, University of Basel,
Klingelbergstrasse 80 , CH-4056 Basel, Switzerland ETH Zurich]Laboratory of
Physical Chemistry, ETH Zurich, 8093 Zurich, Switzerland University of
Basel]Department of Chemistry, University of Basel, Klingelbergstrasse 80 ,
CH-4056 Basel, Switzerland

##  1 Introduction

Tunneling splittings are exquisitely sensitive to the accuracy of a molecular
potential energy surface (PES). The nuclear wave-functions corresponding to
the two or multiple quantum mechanical bound states involved in the split
energy levels probe an extended region on the underlying PES. Furthermore, the
tunneling splitting also informs about the barrier height and the shape of the
PES in the region connecting the two wells, see Figure 1. Due to all the
above, tunneling splittings constitute a meaningful and stringent test of the
level of theory at which the underlying PESs were calculated and the accuracy
of their representation required for simulations from which the splittings are
determined.

Figure 1: Schematic illustration of two symmetric double-well potentials with
the same energy barrier $E_{\rm B}$ but differing curvatures and tunneling
splittings $\Delta$. Malonaldehyde and its instanton path are shown.

Even if a PES is given, accurate computation of tunneling splittings for
multidimensional systems from quantum-based methods itself is a formidable
task. The ring-polymer instanton (RPI) approach provides a semiclassical
approximation of a tunneling process and can be used to calculate tunneling
splittings in molecular systems.1, 2, 3, 4 As was shown for the formic acid
dimer4, it is necessary to include all degrees of freedom of the molecule for
a quantitative comparison with experiment. This often means that the
(semiclassical) full-dimensional instanton approximation is more accurate than
a reduced-dimensional quantum calculation. Instanton theory is based on the
path-integral formulation of quantum mechanics and is formally exact only in
the limit of $\hbar\rightarrow 0$.5 In many previous studies it has been shown
to give predictions within about 20% of fully quantum-mechanical calculations
using the same PES for typical molecular systems, as long as the barrier
height is significantly higher than the zero-point energy along the tunneling
mode.1, 4, 6 Instanton calculations which, contrary to exact quantum
calculations such as wave-function propagation,7, 8 scale well with system
size, are often used in combination with analytical, full-dimensional PESs.
Path-integral molecular dynamics (PIMD) simulations also scale well with
system size but are considerably more expensive than an instanton
calculation.9, 10

In principle it is possible to implement the instanton approach using on-the-
fly ab-initio electronic structure calculations.11, 12 However, because
energies, gradients and Hessians are needed for each ring-polymer bead, this
may be impractical for medium-sized molecules if high accuracy from coupled
cluster with perturbative triples (CCSD(T)) level of theory is sought. More
recent work has been devoted to combining machine learning (ML) and instanton
theory to reduce overall computational expense. Gaussian process regression
(GPR) has been used to obtain a local fit of the PES around the dominant
tunneling pathway to calculate rate constants.13 It has been shown that the
GPR rate constants are on par with the ab initio results, however, reducing
the number of required electronic structure calculations by an order of
magnitude. Similarly, instanton rate theory has been combined with NNs to
obtain the PES more efficiently as compared to the on-the-fly approach14, 15.

As an alternative, full-dimensional PESs can now be constructed for medium-
sized molecules from which tunneling splittings can also be determined using
the instanton approach. The generation of ML PESs based on large data sets of
ab initio data is a challenging task16 and accuracy as well as the level of
theory of the PES is of crucial importance for the accurate determination of
tunneling splittings. The “gold standard” CCSD(T) approach scales as $N^{7}$
(with $N$ being the number of basis functions)17 which quickly becomes
computationally prohibitive for generating data to build a full-dimensional
PESs even for relatively small molecules. To avoid the need for calculating
large ab initio data sets at high levels of theory transfer learning (TL) 18,
19, 20 and related $\Delta$-ML21 were shown to be data and cost-effective
alternatives22, 23, 24, 25, 26, 27.

The combination of TL and instanton theory appears particularly appealing as
the instanton path (IP) can be determined on a low-level PES, which gives a
rough approximation to the true tunneling path, and can be included (and
iteratively refined if needed) into the TL data set. Additionally, the IP is
inherently local and, thus, allows concentrating on improving only a small
part of a PES. While instanton theory has been used in combination with ML
schemes13, 28, 15, the present work demonstrates the first combination of
instanton theory with TL. The capability of the combined approach is
demonstrated for the extensively studied malonaldehyde system exhibiting
intramolecular hydrogen transfer (HT).

Ring Polymer Instanton theory has been employed to calculate the tunneling
splitting of malonaldehyde on a permutationally invariant polynomial (PIP) PES
fitted to 11147 near basis-set-limit frozen-core CCSD(T) electronic
energies29. The splitting was found to be 25 cm-1 with RPI 6 as well as with a
strongly related instanton method 30. The same PES was also used to calculate
tunneling splittings using the fixed-node diffusion Monte Carlo (DMC) method
giving 21.6 cm-1 with a statistical uncertainty of 2 to 3 cm-1.29 For DT,
computed values on the PIP PES are 3.3 and 3.4 cm-1 using RPI6 and the related
instanton method30, and $3.0\pm 2-3$ cm-1 from DMC simulations.29 The
tunneling splitting from RPI calculations on a LASSO fit to CCSD(T)(F12*)
energies was found to be 19.3 cm-1.6, 31 To validate such computations, direct
comparison with experiment is also of interest. Reliable tunneling splittings
from experiment are only available for a few select systems32, 33, 34, 35, 36,
37, 38, 39. For malonaldehyde, the experimentally determined tunneling
splitting is 21.583 cm-1 and 2.915 cm-1 for HT and deuterium transfer (DT),
respectively32, 33, 40. Although these results are not in perfect agreement
with experiment, they are close enough that spectroscopic assignments can be
made and provide detailed mechanistic information about the tunneling process.

The purpose of the present work is to develop and quantitatively assess an
evidence-based procedure to determine reliable tunneling splittings by
combining transfer learned PESs with instanton calculations. It is shown that
this can dramatically reduce the cost of the overall simulation in comparison
to working with _ab initio_ potentials. The work is structured as follows.
First, the methods and the generation of the data sets are presented. This is
followed by a thorough evaluation of the accuracy of the transfer leaned PESs
in terms of tunneling splittings and harmonic frequencies. Finally, the
results are discussed and conclusions are drawn.

##  2 Methods

###  2.1 Ring-Polymer Instanton Theory

In a one-dimensional model, instanton theory is strongly related to the WKB
approximation41. Its main advantage, however, is that it can also be applied
to multidimensional systems, in which it locates the uniquely defined optimal
tunneling pathway.42 This pathway, known as the instanton, is defined as a
long imaginary-time $\tau\rightarrow\infty$ path connecting two degenerate
wells which minimizes the action, $S$. In computations, the path is located
using an efficient ring-polymer optimization based on discretizing the path
into $N$ ring-polymer beads and taking the limit $N\rightarrow\infty$
(typically on the order of 1000 is sufficient for convergence). The action is
determined by the distance between neighbouring beads as well as the
potential-energy of each bead, i.e. it uses information _along_ the IP. Full
technical details are presented in previous work.1, 2 In general, the IP is
not equivalent to the minimum-energy pathway (MEP) and will not even pass
through a saddle point. This is because the instanton finds a compromise
between length and height to optimize the tunneling path. Unlike PIMD or DMC,
no random numbers or statistical errors are involved and so the instanton
(once it is converged with $\tau\rightarrow\infty$ and
$N/\tau\rightarrow\infty$) is in principle uniquely determined by the PES.

Once the IP has been located, fluctuations around the path are computed to
second order and the information is combined into the term $\Phi$, i.e. this
is based on information _around_ the IP. For this, one requires the Hessians
(second-derivative matrix of PES) at each bead. The final prediction for the
tunneling splitting (in a double-well system) is given by

$\displaystyle\Delta=\frac{2\hbar}{\Phi}\sqrt{\frac{S}{2\pi\hbar}}\,\mathrm{e}^{-S/\hbar}.$
(1)

Because $S$ appears in the exponent it is particularly important to determine
this quantity with high accuracy.

The method has also been generalized to treat tunneling in systems with
multiple (more than two) wells 43, 3 and in cases with non-degenerate wells
for instance due to asymmetric isotopic substitution.6 The approach outlined
in this work is, in principle, also applicable to these extensions.

###  2.2 Machine Learning

All PESs used in this work are represented with a high-dimensional neural
network (NN) of the PhysNet44 architecture. PhysNet is a ‘message-passing’45
NN that employs learnable descriptors of the atomic environments to predict
individual atomic energy contributions $E_{i}$ and partial charges $q_{i}$.
The descriptors are initialized as ${\bf x}_{i}^{0}={\bf e}_{Z_{i}}$, where
${\bf e}_{Z_{i}}$ corresponds to a parameter vector defined by the nuclear
charge $Z_{i}$, i.e. atoms of the same element share the same descriptor. The
descriptor is then iteratively updated and refined to best describe the local
chemical environment of each atom $i$ by passing ‘messages’ between atoms
within a cut-off $r_{\rm cut}$ following

$\displaystyle{\bm{x}}_{i}^{l+1}={\bm{x}}_{i}^{l}+\sum_{r_{ij}<r_{cut}}\mathcal{F}({\bm{x}}_{i}^{l},{\bm{x}}_{j}^{l},r_{ij}).$
(2)

where $r_{\rm cut}$ was 10 Å. Here, ${\bm{x}}_{i}^{l}$ and ${\bm{x}}_{j}^{l}$
are the descriptors of atoms $i$ and $j$ at iteration $l$, $r_{ij}$ is their
interatomic distance and
$\mathcal{F}({\bm{x}}_{i}^{l},{\bm{x}}_{j}^{l},r_{ij})$ is the ‘message-
passing’ function (for details see Ref. 44). Because only pairwise distances
are used to encode the atoms’ chemical environment and summation is
commutative, the resulting descriptors (and thus the PES) are invariant with
respect to translation, rotation and permutation of identical atoms, which is
of particular importance when describing tunneling between degenerate wells.
The descriptors are then used to predict partial charges $q_{i}$ (which are
corrected to ensure total charge conservation) and the total energy of the
chemical system by summation of the atomic contributions and explicitly
including long-range electrostatics according to

$\displaystyle
E=\sum_{i}E_{i}+k_{e}\sum_{i=1}^{N}\sum_{j>i}^{N}\frac{q_{i}q_{j}}{r_{ij}}$
(3)

Here, $k_{e}$ represents Coulomb’s constant and the second term involving
$\frac{q_{i}q_{j}}{r_{ij}}$ is damped to avoid numerical instabilities caused
by the singularity at $r_{ij}=0$ (for details refer to Ref. 44). The forces
${\bm{F}}$ and Hessians ${\bm{H}}$ can be obtained analytically using reverse
mode automatic differentiation46 as implemented in Tensorflow47.

The learnable parameters of PhysNet are fitted to reference ab initio
energies, forces and dipole moments following the strategy outlined in
Reference 44. The partial charges $q_{i}$ are fitted to the ab initio dipole
moment ($\mu=\sum_{i}^{N}q_{i}{\bm{r}_{i}}$) and explicitly enter the energy
expression (see equation 3). In the present work, the TL scheme is employed
whereby the parameters of a low-level (LL) treatment are used as a meaningful
initial guess and are fine-tuned using higher-level information. For TL, the
learning rate is reduced from $10^{-3}$ (as for learning a model from scratch)
to $10^{-4}$. The LL in the present work is the full-dimensional PES for
malonaldehyde at the MP2/aug-cc-pVTZ level of theory (henceforth, PhysNet MP2
PES) which is available from previous work and was trained on $\sim 70000$
reference structures.22 This PhysNet PES has a barrier for HT of 2.79 kcal/mol
which compares to a reference value of 2.74 kcal/mol calculated at the
MP2/aug-cc-pVTZ level of theory and the reference harmonic frequencies are
reproduced with a root-mean-square deviation (RMSD) of 3.6 cm-1. The high-
level (HL) treatment is the considerably higher and computationally much more
demanding CCSD(T)/aug-cc-pVTZ level of theory at which energies, forces and
dipole moments are calculated using Molpro48 for all data points used in TL.

###  2.3 Data Set Generation

Transfer learning requires high-level energies, forces and dipole moments for
selected geometries of the system considered and ideally cover all spatial
regions relevant for the observable(s) of interest. Without additional a
priori information it is advantageous to generate an initial pool of
structures which can be used for TL to fine-tune the LL treatment. When
selecting configurations, it is not necessary to sample both potential wells
since PhysNet handles this symmetry by construction. Here, the initial pool
contained 862 malonaldehyde configurations consisting of:

  * •

111 geometries along the MEP of the PhysNet MP2 PES.

  * •

110 geometries along the IP of the PhysNet MP2 PES.

  * •

111 geometries along the IP determined on a PES that was transfer learned by
using CCSD(T) information of the 111 MEP geometries (see above) to have an
energy barrier closer to the ab initio CCSD(T) barrier.

  * •

280 geometries obtained from normal mode sampling (NMS) around the equilibrium
geometry. For this purpose, normal mode vectors and corresponding force
constants are determined ab initio at the MP2/aug-cc-pVTZ level of theory.

  * •

240 geometries around the IP as obtained from NMS.

  * •

10 geometries along the IP of TL1 (see Section 3.2).

This data set is referred to as the “Extended Data Set” and transfer learned
models using it are called TLext. To probe the dependence of barrier heights
and tunneling splittings on details of the training, ten independent models
were trained on different splits of the data for TLext (and all the subsequent
TLs). For each of the ten resulting PESs an instanton calculation was carried
out. From this information, averages and standard deviations for the barrier
heights and tunneling splittings were determined.

After validating the performance of TLext from instanton calculations on each
of the independently trained models, smaller subsets of the Extended Data Set
were selected, employed for TL and subsequent tunneling splitting
calculations.

##  3 Results

To set the stage, the tunneling splittings for malonaldehyde were calculated
on the PhysNet MP2 PES using RPI theory. The tunneling splitting calculations
were carried out with three different values of the imaginary time, $\tau$,
corresponding to effective ‘temperatures’ $T=\hbar/k_{B}\tau\in[50,25,12.5]$ K
and with different numbers of beads $N\in[2^{5},..,2^{12}]$ to ensure
convergence. Formally the instanton result is defined in the low-temperature
limit, which is equivalent to infinitely-long imaginary times. The results are
summarized in Table S1. A tunneling splitting of 96 cm-1 is obtained compared
with 25 cm-1 from instanton calculations30, 6 on the PIP-representation29 of
the CCSD(T) reference data and 21.6 cm-1 from experiments.32, 33, 40 This
illustrates the insufficient quality of the MP2 level of theory to capture
tunneling splittings correctly. For the following, all instanton calculations
were carried out with $N=4096$ beads at an effective temperature of $T=25$ K,
which was found to be more than sufficient for convergence of $\Delta$ to two
significant figures.

Figure 2: Data sets used for TL projected onto a 2D cut through the TLext PES
spanned by the O–O distance and the reaction coordinate $q=r_{\rm
O_{A}H}-r_{\rm O_{B}H}$, for labels see Fig. 1. The Extended Data Set (862
structures, gray circles) is shown together with sets for TL0 (25 structures,
yellow crosses), TL1 (50 structures, turquoise circles), and TL2 (100
structures, salmon squares). The MEP and the instanton path (as determined on
the TLext PES) are marked with a dash-dotted and a dashed line, respectively.
The IP on the PhysNet MP2 PES is the white dashed line and clearly differs
from that (black dashed) on the CCSD(T) PES.

###  3.1 Performance of TLext

As a reference for the following exploration, the performance of TLext using
the full set of 862 energies, forces and dipole moments determined at the
CCSD(T)/aug-cc-pVTZ level is first assessed. These geometries are shown as a
projection onto the PES spanned by the O–O distance and the reaction
coordinate $q=r_{\rm O_{A}H}-r_{\rm O_{B}H}$ in Figure 2 (gray circles). As
the PES is symmetric with respect to $q=0$, the same geometry appears to the
left and to the right of the mirror plane. The Extended Data Set was split
according to $\sim 80/10/10$ % into training/validation/test sets, from which
the test sets were used only for testing. Across the 10 TL models, the
separate test sets were predicted on average with MAE($E$) $\approx 0.006$,
RMSE($E$) $\approx 0.009$ kcal/mol, MAE($F$) $\approx 0.03$ and RMSE($F$)
$\approx 0.07$ kcal/mol/Å. The average barrier height on the ten transfer-
learned PESs ($\langle{\rm TL}_{\rm ext}\rangle$) was $E_{\rm B}=3.8945\pm
0.0006$ kcal/mol which compares with an ab initio barrier of 3.8948 kcal/mol
determined at the CCSD(T)/aug-cc-pVTZ level of theory determined from present
calculations. The RPI tunneling splittings for HT and DT were $\Delta_{\rm
H}=25.3\pm 0.2$ cm-1 and $\Delta_{\rm D}=3.7\pm 0.03$ cm-1, respectively, see
Table 1 and S2. These results compare with computed splittings on the PIP
representation of CCSD(T)/aug-cc-pVTZ reference calculations using instanton
calculations that yield 25/3.4 cm-130, 6 and experimental splittings of
21.6/2.9 cm-1.32, 33, 40 Hence, the transfer-learned PES using fewer than 1000
higher level CCSD(T)/aug-cc-pVTZ energies and forces together with the same
method for determining the tunneling splitting performs on par with
calculations on the PIP representation of the $\sim 11000$ CCSD(T)/aug-cc-pVTZ
energies.29 Based on this it is of much practical interest to further reduce
the number of HL calculations required to achieve the same result. Therefore,
in a next step, different subsets of the Extended Data Set are considered and
used for TL to arrive at an ideally small number of HL points while still
retaining the accuracy in tunneling splittings from instanton calculations.

| $N_{\rm data}$ | $E_{\rm B}$ | $\Delta_{\rm H}$ | $S/\hbar$ | $\Phi$
---|---|---|---|---|---
MP2 | 70k | 2.7889 | 96.3 | 4.502 | 42.770
$\langle{\rm TL}_{0}\rangle$ | 25 | $3.8925\pm 0.0119$ | $23.4\pm 1.9$ | $5.749\pm 0.026$ | $57.405\pm 4.486$
$\langle{\rm TL}_{1}\rangle$ | 50 | $3.8974\pm 0.0161$ | $24.9\pm 1.1$ | $5.764\pm 0.012$ | $53.096\pm 2.375$
$\langle{\rm TL}_{2}\rangle$ | 100 | $3.8941\pm 0.0021$ | $25.2\pm 0.5$ | $5.748\pm 0.005$ | $53.053\pm 1.156$
$\langle{\rm TL}_{\rm ext}\rangle$ | 862 | $3.8945\pm 0.0006$ | $25.3\pm 0.2$ | $5.743\pm 0.002$ | $53.241\pm 0.454$
Table 1: Energy barriers $E_{\mathrm{B}}$ (kcal/mol), tunneling splittings
$\Delta$ (in cm-1) at $T=25$ K and with $N=4096$, action $S/\hbar$ and
fluctuation factor $\Phi$ (in a.u. of time $\hbar/E_{h}$) for malonaldehyde
determined from TLext and TL0,1,2 PESs. The ab initio barrier height for HT at
the CCSD(T)/aug-cc-pVTZ level of theory is 3.8948 kcal/mol.

###  3.2 Performance on Smaller Datasets: TL0, TL1, and TL2

From the Extended Data Set containing 862 geometries, different subsets were
extracted. The size of the data set has to be chosen small enough for
efficient computation but sufficiently large to still cover the appropriate
regions of configurational space probed by the instanton calculation.

TL0: To check whether a considerably smaller data set suffices as a starting
point, ten TLs were preformed on a data set containing only 25 geometries: a)
5 IP geometries (approximately equally spaced) determined on a PES that was
transfer learned to have a barrier closer to the ab initio CCSD(T) barrier; b)
10 geometries, each, selected from the NMS around the equilibrium geometry and
the IP. This can be done by selecting geometries based on a RMSD criterion on
the $n$ atomic positions
($\sqrt{\frac{1}{n}\sum_{i=1}^{n}||v_{i}-w_{i}||^{2}}$ where $v_{i}$ and
$w_{i}$ are two sets of Cartesian coordinates of atom $n$), which was done as
follows for both groups of geometries that were generated with NMS. Starting
from a random geometry, new geometries are added iteratively if the RMSD with
respect to the selected ones is larger than a threshold. For this reason, the
threshold is maximized to include 10 geometries. Note that no MEP geometries
are added. The data set for TL0 are the yellow crosses in Figure 2.

With this smallest subset the barrier height of the (ensemble of the) transfer
learned PES $\langle{\rm TL}_{0}\rangle$ is $E_{\rm B}=3.8925\pm 0.012$
kcal/mol which is, within errors, close to the target value of 3.8948 kcal/mol
determined at the CCSD(T)/aug-cc-pVTZ level. From 10 independent instanton
calculations the average splitting is $\Delta_{\rm H}=23.4\pm 1.9$ cm-1 which
is $\sim 2$ cm-1 below that from the simulations on the TLext PES but still
within statistical fluctuation. Comparing the action $S/\hbar$ of the IPs on
TL0 and TLext shows that they are comparable ($5.743\pm 0.002$ vs. $5.749\pm
0.026$) and even identical within the error bars. However, the uncertainty on
TL0 is larger by an order of magnitude compared with that on TLext. Thus, the
action, $S/\hbar$, of the path is clearly less well defined on TL0. For the
fluctuation factor $\Phi$ the differences are considerably larger between the
two families of PESs. Still, the values themselves are within error bounds but
again, the fluctuation around the mean for TL0 is ten times larger than that
for the TLext PESs. This conclusion also hold for DT, see Table S2. Overall,
using only 25 additional data points as done for TL0 already yields
encouraging results for the barrier height and tunneling splittings. To
explore further improvements new points were added and the process was
repeated.

Figure 3: Tunneling splittings for H- and D-transfer (left and right panels)
from all TL PESs (transparent circles). The corresponding averages (opaque
circle) and standard deviations (error bars as $\pm\sigma$) as obtained from
TL0 (blue), TL1 (red), TL2 (green) and TLext are reported, too. Comprehensive
lists of the tunneling splittings for all TL PESs are given in Tables S3 to
S6.

TL1: For TL1 the points used for TL0 were extended and increased to 50 points,
see turquoise circles in Figure 2. The data set for TL1 contained: a) 5
geometries along the MEP of the PhysNet MP2 PES; b) 5 geometries along the IP
of a transfer learned PES from using the MEP points calculated at CCSD(T); c)
20 geometries, each, selected from the NMS around the equilibrium geometry and
the IP. The MEP and IP geometries are chosen with a uniform spacing along the
respective path and the geometries from NMS are selected following the RMSD
approach outlined for TL0.

Figure 4: Action $S/\hbar$ (left panel) and fluctuation factor $\Phi$ (right
panel) from all TL PESs (transparent circles) for HT (for DT see Figure S5).
The corresponding averages (opaque circle) and standard deviations (error bars
as $\pm\sigma$) as obtained from TL0 (blue), TL1 (red), TL2 (green) and TLext
are indicated as well.

TL1 yields an averaged barrier height $E_{\rm B}=3.8974\pm 0.0161$ which
agrees within error bars with that of TLext and the ab initio value (3.8948
kcal/mol). The splitting $\Delta_{\rm H}=24.9\pm 1.1$ cm-1 is only 0.4 cm-1
below that of TLext. The key improvement is that the fluctuation factor $\Phi$
agrees considerably better with TLext than for TL0 and the remaining 2%
discrepancy can be traced to the $0.02$ absolute difference in $S/\hbar$.
Overall, increasing the number of geometries used for TL to $N_{\rm data}=50$
in this fashion leads to a PES which reproduces the barrier height and
splittings from TLext. To further probe convergence of these results yet a
larger data set was considered.

TL2: While the accuracy of TL1 might appear satisfactory, convergence of TL1
cannot be checked without TLext. Thus, TL1 was further extended to yield TL2.
This was accomplished from a strategy related to adaptive sampling.49 Two
independent models from TL1 were used to predict the energy $E$ of the
remaining pool of geometries obtained from NMS. From these geometries, the 40
geometries with the largest deviation (here $\sim 0.02$ kcal/mol) between the
prediction of the two NNs are added to the data set. If a large deviation
between the energy predictions of the two models is found, it is likely that
no or too little reference data has been included in TL1. In addition, 5
geometries along the MEP and instanton of the TL1 were added such that TL2
contained a total of 100 geometries (see salmon squares in Figure 2).

With TL2 the barrier height $E_{\mathrm{B}}$, splitting $\Delta$, action
$S/\hbar$ and fluctuation $\Phi$ further improve over those from TL1 and are
closer to the results from TLext, see Table 1 and S2. Within $1\sigma$ all
values for HT and DT agree with those from TLext and with the reference from
the literature except for the tunneling splitting for DT. Overall, addition of
50 to 100 additional points from the HL treatment appears to suffice to arrive
at a quantitatively correct PES transfer-learned from the LL treatment
(MP2/aug-cc-pVTZ). Figure 3 illustrates the gradual convergence for TL0 to TL2
towards the results found for TLext, which is promising. While exceptional
agreement regarding the energy barrier $E_{\mathrm{B}}$ (as the “simpler”
property) is found for all the TLs the standard deviation of the splittings
$\Delta_{\rm H}$ (which is more challenging to obtain) can be reduced from 1.9
cm-1 (7.5%) for TL0 to 0.5 cm-1 (2%) for TL2. The results from TL2 are
accurate to within 0.1 cm-1 for HT and 0.04 cm-1 for DT as compared to ${\rm
TL}_{\rm ext}$. This corresponds to deviations of 0.4% and 1.1%. The IPs
themselves are reported in Figures S1 to S4 and show slight deviations between
the different TLs on the smallest data set (TL0), but for TL1, TL2 and TLext
they are hardly distinguishable.

In summary, it has been found that with between 25 and 50 HL energies and
forces for judiciously chosen structures the correct barrier height, tunneling
splitting, action and fluctuation can be obtained, see Table 1. This is the
foundation to further optimize the procedure, ideally based entirely on
information available from the LL surface by minimizing the amount of data
required from the HL treatment.

###  3.3 Towards an Optimized Procedure

For an even more efficient procedure, an approach is sought that is based on
information about the LL-PES only. Hence, an attempt is made to further reduce
the computational effort by minimizing the number of structures for which HL
calculations need to be carried out for an improved PES and tunneling
splittings. Therefore, it is explored which elements of the procedure are most
important for obtaining high accuracy to cost ratios. For moving towards a
more evidence-guided, optimized procedure, TL is carried out from LL-
information that is only contained in the MEP and IP as follows: a) only the
MEP (TLMEP) using 111 geometries along the MEP of the PhysNet MP2 PES; b) only
the IP (TLIP) using 110 geometries along the IP of the PhysNet MP2 PES; and c)
a combination of the MEP and the IP (TLMI) including 111 MEP and 110 instanton
geometries as obtained from the PhysNet MP2 PES. The HL information for TL
consisted again of energies, forces and dipole moment determined at the
CCSD(T)/aug-cc-pVTZ level of theory. A total of 5 TLs were performed, each on
different splits of the data. The two best NNs (judged from the performance on
the validation set) were used for further analysis.

|  $E_{\mathrm{B}}$ [kcal/mol] |  $\Delta_{\rm H}$ [cm-1] | $S/\hbar$ | $\Phi$
---|---|---|---|---
TL${}_{\rm MEP}^{a}$ | 3.94 | 15.4 | 5.887 | 76.482
TL${}_{\rm MEP}^{b}$ | 3.95 | 14.1 | 5.912 | 81.921
TL${}_{\rm IP}^{a}$ | 3.49 | 16.7 | 5.792 | 76.975
TL${}_{\rm IP}^{b}$ | 3.50 | 16.9 | 5.794 | 75.858
TL${}_{\rm MI}^{a}$ | 3.90 | 18.2 | 5.768 | 72.262
TL${}_{\rm MI}^{b}$ | 3.90 | 16.9 | 5.774 | 77.374
$\langle{\rm TL}_{\rm ext}\rangle$ | 3.8945 | 25.3 | 5.743 | 53.241
Table 2: Energy barriers $E_{\mathrm{B}}$, tunneling splittings $\Delta_{\rm
H}$ (at $T=25$ K and $N=4096$), action $S/\hbar$ and fluctuation factor $\Phi$
(in a.u. of time $\hbar/E_{h}$) for HT determined from TL PESs using MEP
points only (TLMEP), instanton points only (TLIP), and a combination of MEP
and instanton points (TLMI). TL${}_{\rm x}^{a}$ and TL${}_{\rm x}^{b}$
correspond to two NNs that are trained on different splits of the data. Note
that the energy barrier for TLIP is inaccurate, as expected, because the IP
misses the transition state of the PES. The ab initio barrier at the
CCSD(T)/aug-cc-pVTZ level of theory is 3.8948 kcal/mol.

First, it is noted that for TLMEP and TLMI the barrier $E_{\mathrm{B}}$ for HT
agrees well with the target value of 3.8948 kcal/mol (ab initio CCSD(T)/aug-
cc-pVTZ level value) whereas this is not the case for TLIP, as expected,
because the instanton path does not pass through the transition state of the
MEP (see Table 2). Also, despite starting from MP2 information only, TL to the
HL model yields considerably improved tunneling splittings $\Delta_{\rm H}$
for all three models, ranging from 15 to 18 cm-1, compared to those from the
PhysNet MP2 PES (96 cm-1).

Considering the action $S/\hbar$ from TLMEP, TLIP and TLMI it is seen that it
progressively approaches that from TLext. For TLMEP the action overshoots the
target value from TLext by $\sim 0.15$ which leads to an error of $\sim 15$ %
in the splitting because $S/\hbar$ appears in the exponential factor in Eq. 1.
Conversely, with a difference of 0.04 compared with TLext, the error for TLMI
due to $S/\hbar$ is only $\sim 4$ %. The influence of $\Phi$ on the difference
between TLext and the three models considered in Table 2 is minor because for
all of them $\Phi$ is uniformly too large by $\sim 40$ % compared with that
from TLext.

Table 2 suggests that a combination of information from the MEP and the IP
used for transfer learning, i.e. TLMI, yields $E_{\rm B}$ and $\Delta_{\rm H}$
closest to the results from TLext. However, the tunneling splittings still
differ by more than 10 % from the target value obtained on the TLext PESs.
Considering the actions $S/\hbar$ and fluctuations $\Phi$ for all the transfer
learned PESs in Table 2 it is found that in particular the values of $\Phi$,
which are sensitive to fluctuations around the IP, differ considerably from
that on TLext. Hence as a last improvement points _along_ the MEP and IP are
combined with structures _around_ the IP.

| $N_{\rm data}$ | $E_{\rm B}$ | $\Delta_{\rm H}$ | $S/\hbar$ | $\Phi$
---|---|---|---|---|---
MP2 | 70k | 2.7889 | 96.3 | 4.502 | 42.770
$\langle{\rm TL}_{0}\rangle$ | 25 | $3.8925\pm 0.0119$ | $23.4\pm 1.9$ | $5.749\pm 0.026$ | $57.405\pm 4.486$
$\langle{\rm TL}_{1}\rangle$ | 50 | $3.8974\pm 0.0161$ | $24.9\pm 1.1$ |  $5.764\pm 0.01$2 | $53.096\pm 2.375$
$\langle{\rm TL}_{2}\rangle$ | 100 | $3.8941\pm 0.0021$ | $25.2\pm 0.5$ | $5.748\pm 0.005$ | $53.053\pm 1.156$
$\langle{\rm TL}_{\rm ext}\rangle$ | 862 | $3.8945\pm 0.0006$ | $25.3\pm 0.2$ | $5.743\pm 0.002$ | $53.241\pm 0.454$
$\langle{\rm TL}_{\rm EB1}\rangle$ | 25 | $3.9025\pm 0.0232$ | $23.7\pm 1.1$ | $5.740\pm 0.013$ | $57.099\pm 3.249$
$\langle{\rm TL}_{\rm EB2}\rangle$ | 25 | $3.9041\pm 0.0189$ | $23.5\pm 2.1$ | $5.733\pm 0.025$ | $58.241\pm 4.850$
Table 3: Energy barriers $E_{\mathrm{B}}$ (kcal/mol), tunneling splittings
$\Delta_{\rm H}$ (in cm-1) at $T=25$ K and with $N=4096$, action $S/\hbar$ and
fluctuation factor $\Phi$ (in a.u. of time $\hbar/E_{h}$) for malonaldehyde
determined from TLEB. The ab initio barrier at the CCSD(T)/aug-cc-pVTZ level
of theory is 3.8948 kcal/mol.

For two final, evidence-based TLs (TLEB1 and TLEB2), a set of 25 data points
was generated as follows. A total of 5 points was selected along the MEP and
IP (one close to the minimum and 2 points along the PhysNet MP2 PES MEP and IP
each, see black points in Figure 5). These were supplemented by 20 geometries
from NMS around the equilibrium structure and the IP that are selected by
means of an RMSD criterion. For TLEB1 (orange circles in Figure 5), the
geometries with largest RMSD are selected (single geometries which occupy the
same $r_{\rm OO}$ and $q$ coordinates are eliminated). For TLEB2 (green
crosses in Figure 5), besides the RMSD criterion, the geometries were selected
to cover the important configurational space more regularly (as judged by
Figure 5).

For both sets of points the corresponding CCSD(T)/aug-cc-pVTZ energies, forces
and dipole moments were used for TL, resulting in two sets of transfer learned
PESs: TLEB1 and TLEB2. For both of them the barrier height (3.90 kcal/mol) is
within error bars of TLext for which it was 3.89 kcal/mol. The action
$S/\hbar$ for both EB-models agree with TLext within error bounds although the
fluctuations around the mean is larger by almost an order of magnitude, see
Table 3. For the fluctuation $\Phi$ the differences compared with TLext are
$\sim 5$ %, commensurate with TL0 and evidently improved over those using MEP,
IP or MI, see Table 2 which do not train on geometries around the path. The
tunneling splittings are $\Delta_{\rm H,EB1}=23.7\pm 1.1$ cm-1 and
$\Delta_{\rm H,EB2}=23.5\pm 2.1$ cm-1, both of which are close to/within error
bounds of the reference value $(25.3\pm 0.2)$, see Table 3. These results are
comparable and slightly better to those on TL0 which also was based on only 25
points for TL. However, the training data for TLEB1,EB2 are selected based
entirely on the PhysNet MP2 PES whereas TL0 made use of HL information in that
it employed geometries along the IP of a PES with corrected barrier. Hence,
from a computational perspective, the EB models are considerably more cost-
effective.

Figure 5: The data sets used for TLEB1 (orange circles) and TLEB2 (green
crosses) are shown on a 2D projection of the PES spanned by the O–O distance
and the reaction coordinate $q=r_{\rm O_{A}H}-r_{\rm O_{B}H}$. Both data set
contain 25 geometries. The MEP and the IP (as determined on the TLext PES) are
marked with a dash-dotted and a dashed line, respectively. The IP on the MP2
PES is illustrated as white dashed line.

Overall, it is found that TL to the HL with 25 additional points yields a
barrier height that agrees with the full HL treatment and the tunneling
splitting differs by only $\sim 1$ cm-1. Any further improvement requires
additional points. Based on the results in Table 3 it is expected that using a
TL model trained on fewer than 50 judiciously selected HL data points yields
results within 1% of the HL reference TLext. This needs also be contrasted
with an expected accuracy of instanton calculations for tunneling splittings
of $\sim 20$ %.

##  4 Discussion and Conclusions

The present work aimed at developing a computationally efficient and accurate
road-map for how to improve a given LL model - which was assumed to be
“comprehensive” (here $7\times 10^{4}$ MP2 energies and forces were used) for
the observable of interest - to a HL model by providing a small amount of
additional information at the higher level of theory considering a particular
observable. Here, the observable was the tunneling splitting for HT/DT in
malonaldehyde for which the LL model (MP2/aug-cc-pVTZ) yielded $\Delta=96$
cm-1, compared with a literature value of $\Delta_{\rm H}\sim 25$ cm-1 from a
PIP-represented PES of CCSD(T)/aug-cc-pVTZ reference energies using a range of
methods for computing $\Delta_{\rm H/D}$, see Table S7. Most HL models
generated were based on TL using 10s to 100s of HL points and yield
$\Delta\sim 25$ cm-1 which is a substantial improvement over the LL model and
consistent with computations in the literature at the same level of theory but
employing computationally much more demanding approaches. The remaining
differences between computations and experiments are due to a) shortcomings of
the CCSD(T) level of theory compared with a full CI treatment, b) the
incompleteness of the basis set, and c) inherent semiclassical approximations
of instanton theory (e.g. neglect of coupling to overall molecule rotation and
anharmonicity perpendicular to the instanton path).

A typically used shortcut is to optimize the instanton using a LL ab initio
method, e.g. DFT or MP2, and then compute the CCSD(T) properties along the
path to correct the action $S/\hbar$. Such a hybrid approach was assessed
using the PhysNet MP2 PES to optimize the instanton and then calculate the
action $S/\hbar$ on the TLext PES.11, 50, 12, 51 The results are summarized in
Table 4 and illustrate that although the hybrid approach is, in this case,
able to infer the correct value for $S/\hbar$, the TL approach additionally
improves $\Phi$. Using Equation 1 with the action as determined by the hybrid
approach ($S/\hbar=5.7401$) and the fluctuation factor determined on the
PhysNet MP2 PES ($\Phi=42.7705$) yields $\Delta_{\rm H}\sim 31.5$ cm-1 which
overestimates the value of 25.1 cm-1 from TLext. The TL approach is thus able
to provide a more accurate prediction of $\Delta$ for a similar computational
cost.

| H | D
---|---|---
PES | $\Delta_{\rm H}$ | $S/\hbar$ | $\Phi$ | $\Delta_{\rm D}$ | $S/\hbar$ | $\Phi$
MP2 | 96.3 | 4.502 | 42.771 | 18.8 | 5.705 | 73.964
Hybrid | 31.5 | 5.740 | 42.771 | 4.4 | 7.276 | 73.964
$\langle{\rm TL}_{\rm EB2}\rangle$ | 23.5 | 5.733 | 58.241 | 3.5 | 7.266 | 94.777
TLext | 25.1 | 5.744 | 53.586 | 3.6 | 7.274 | 89.957
$\langle{\rm TL}_{\rm ext}\rangle$ | 25.3 | 5.743 | 53.241 | 3.7 | 7.273 | 89.350
Lit.6 | 25 | 6.129 | 37.794 | 3.3 | 7.790 | 61.392
Table 4: Actions $S/\hbar$ and fluctuation factor $\Phi$ as obtained from the MP2 PhysNet PES (MP2), on a representative TLext PES (TLext) and a hybrid approach optimizing the IP on the PhysNet MP2 PES and using TLext to obtain CCSD(T) properties along resulting IP6 Mode | MP2 | $\langle{\rm TL}_{0}\rangle$ | $\langle{\rm TL}_{1}\rangle$ | $\langle{\rm TL}_{2}\rangle$ | $\langle{\rm TL}_{\rm EB1}\rangle$ | $\langle{\rm TL}_{\rm ext}\rangle$ | CCSD(T)
---|---|---|---|---|---|---|---
1 | 277.49 | 266.68 | 265.48 | 265.28 | 266.93 | 265.17 | 264.71
2 | 286.59 | 285.58 | 283.96 | 280.88 | 286.01 | 283.77 | 281.85
3 | 394.33 | 389.41 | 387.06 | 389.09 | 392.05 | 389.70 | 389.12
4 | 514.07 | 501.82 | 502.93 | 503.23 | 505.62 | 505.50 | 505.07
5 | 789.38 | 772.93 | 774.86 | 775.82 | 773.14 | 775.33 | 775.06
6 | 888.62 | 885.17 | 886.17 | 886.92 | 885.18 | 886.39 | 886.07
7 | 937.61 | 906.84 | 907.54 | 908.64 | 908.10 | 908.03 | 907.18
8 | 1012.29 | 988.09 | 989.05 | 992.64 | 991.79 | 990.78 | 989.73
9 | 1023.82 | 1005.95 | 1004.01 | 1002.71 | 1006.17 | 1002.93 | 1002.60
10 | 1048.75 | 1039.62 | 1039.69 | 1038.19 | 1039.50 | 1037.56 | 1037.85
11 | 1109.69 | 1104.46 | 1104.81 | 1102.64 | 1102.08 | 1102.20 | 1101.03
12 | 1288.28 | 1274.86 | 1271.94 | 1272.65 | 1272.27 | 1273.41 | 1272.73
13 | 1403.09 | 1400.18 | 1399.97 | 1401.56 | 1402.15 | 1401.73 | 1400.73
14 | 1407.97 | 1404.80 | 1408.40 | 1407.17 | 1406.59 | 1408.81 | 1406.94
15 | 1482.06 | 1458.89 | 1463.84 | 1463.95 | 1462.21 | 1467.70 | 1469.33
16 | 1641.52 | 1624.37 | 1627.95 | 1630.14 | 1631.91 | 1633.35 | 1632.52
17 | 1692.91 | 1681.91 | 1687.06 | 1691.27 | 1686.72 | 1693.59 | 1693.63
18 | 3039.02 | 3003.37 | 2999.68 | 2999.67 | 3005.82 | 3000.03 | 3001.26
19 | 3107.12 | 3176.46 | 3178.76 | 3180.11 | 3176.23 | 3176.74 | 3176.35
20 | 3217.85 | 3229.25 | 3228.54 | 3228.74 | 3232.47 | 3226.75 | 3227.12
21 | 3267.30 | 3259.40 | 3260.39 | 3262.98 | 3254.33 | 3263.25 | 3260.27
MAE | 14.62 | 3.01 | 1.93 | 1.54 | 2.61 | 0.89 |
Table 5: Averaged harmonic frequencies calculated from PhysNet potentials and
using ab initio techniques are given in cm-1. As judged from the MAE($\omega$)
the PhysNet model trained on the Extended Data Set containing 862 geometries
is the most accurate, followed by the TL2 (100 data points), TL1 (50 data
points) TLEB1 (25 data points) and TL0 (25 data points). The ab initio
harmonic frequencies obtained at the MP2 level are shown for comparison.

For calculating the tunneling splittings based on the instanton approach it
was found that an evidence-based approach starting from MEP and IP on the LL
PES, augmented with geometries drawn from a pool of structures selected such
that their RMSD is maximal with respect to an existing set of structures
requires of the order of 50 points at the HL for TL. Therefore, only local and
not global knowledge of the PES is required as would, e.g. be necessary for
fully quantum-mechanical methods such as wavepackets. The bottleneck to a
“direct” ab initio-based instanton approach is typically the calculation of
$N$ Hessians, as these are rather expensive to compute.11, 12 Earlier work on
instanton rate theory combined with machine-learning techniques for the H +
CH4 and H + C2H6 reactions required $\sim 50$ energies and forces and 8
Hessians in the training set to converge the rate constant to within $1$ % of
the ab initio result at 200 K.13 Using TL, calculating any high-level _ab
initio_ Hessians at all has been avoided. As is demonstrated here, this can
significantly lower the computational expense with no loss of accuracy.

With regards to the accuracy of the TL PESs it is of interest to compare their
performance on out-of-sample structures. For this a test set was generated
from MD simulations at 700 K on one of the TLext PESs from which 100
geometries were randomly extracted. In addition, 10 equally spaced off-grid
geometries along the IP on the same PES were selected. The CCSD(T)/aug-cc-pVTZ
energies of these 110 geometries cover a range from $\sim 5$ to $40$ kcal/mol
above the global minimum. The energies for these structures were computed
based on TLext (most rigorous TL using 862 HL structures) and TLEB1 (following
the recommended procedure; TL with 25 HL energies and forces), respectively,
and the [MAE100($E$),MAE10($E$)] for the two out-of-sample sets are
[0.21,0.004] kcal/mol and [0.34,0.005] kcal/mol. Notably, the energies of the
geometries used for TLext and TLEB1 only cover a range 20 kcal/mol above the
global minimum whereas the out-of-sample energies reach twice as high, up to
40 kcal/mol above the minimum. Hence, the out-of-sample structures contain
true predictions on the HL-PES. As a comparison, for the PIP PES, which used
energies only and no forces, the reported fitting errors (i.e. in-sample) are
32 cm-1 (0.09 kcal/mol) for energies below 2000 cm-1 (5.7 kcal/mol) above the
global minimum and 211 cm-1 (0.6 kcal/mol) for energies up to 20000 cm-1 (51.2
kcal/mol).29

TL as used in the present work - namely as a local refinement of a LL-PES -
can also be regarded as a variant of the more global “morphing” approach for
PESs.52 It is therefore of interest to consider in what way observables other
than the tunneling splitting change upon TL from LL to HL. For this, harmonic
frequencies were determined for a number of transfer learned PESs. The
harmonic frequencies averaged over the 10 individually trained NNs for
different TLs are reported in Table 5, where they are compared with
frequencies determined from CCSD(T)/aug-cc-pVTZ calculations at the
corresponding equilibrium structure of malonaldehyde. As judged from the
MAE($\omega$) the PhysNet model for TLext is most accurate (MAE $<1$ cm-1),
followed by TL2 (MAE $<2$ cm-1), TL1 (MAE $\sim 2$ cm-1)and TL0 (MAE $\sim 3$
cm-1), as expected, and show a considerable improvement over the MP2
frequencies. For TLEB1 the MAE is $<3$ cm-1. Thus, TL to the HL model also
improves the shape of the PES in degrees orthogonal to the two reaction
coordinates considered for the tunneling splitting.

Figure 6: Flowchart of the recommended TL + RPI procedure.

In terms of a recommended procedure it is noted that the strategy outlined in
going from TL1 to TL2 (adaptive sampling/active learning) can also be pursued
recursively from a “pool” of geometries generated from sampling the PhysNet
MP2 PES. This procedure can be repeated until convergence of the barrier
height and the tunneling splittings. The approach proposed for future
application is (see Figure 6): i) create a LL-PES from a fine grid (here
$7\times 10^{4}$ points) and train a ML model (here PhysNet) ii) generate a
pool of structures based on the LL-PES (including MEP, Instanton, NMS) iii)
choose $N\sim 25$ geometries following EB1/2 and determine energies and forces
from HL ab initio calculations iv) perform TL and instanton calculations on
the HL-PES v) refine the data set using adaptive sampling and structures along
the new IP. Then repeat TL and instanton calculations vi) repeat iv) and v)
until convergence. The present work demonstrates that following such a road-
map requires $\sim 50$ HL energy and force evaluations to determine an
accurate tunneling splitting of malonaldehyde which is manifestly more
efficient than previously explored approaches.

In summary, given that LL models can be constructed efficiently even for
moderately sized molecules such as malonaldehyde or larger,53 the present work
confirms that with specific, evidence-based information grounded in physical
understanding of the process in question, several 10 points from a HL
treatment are sufficient to generate high-quality PESs for a target observable
which was the tunneling splitting in malonaldehyde in the present work. The
MAE($E$) for the TL-PESs trained on energies spanning 20 kcal/mol above the
global minimum is $\sim 0.3$ kcal/mol on off-grid structures spanning 40
kcal/mol and the resulting harmonic frequencies agree to within 1 to 3 cm-1
with rigorous and very time consuming normal mode analysis at the CCSD(T)/aug-
cc-pVTZ level of theory. The recommended approach deduced from the present
work is based on information about the MEP, the IP, and fluctuations around
the IP determined on the LL-PES for which HL calculations are required for TL
to determine the HL-PES. It is expected that - with suitable adaptations due
to the particular observable considered - the present approach can also be
applied to other observables that are computationally expensive to determine
for a given PES, e.g. the quantum bound states of molecules or scattering
cross sections for gas phase reactions from wavepacket propagation.

## Acknowledgment

This work has been financially supported by the Swiss National Science
Foundation NCCR-MUST (MM and JOR), the AFOSR, and the University of Basel.

## References

  * Richardson and Althorpe 2011 Richardson, J. O.; Althorpe, S. C. Ring-polymer instanton method for calculating tunneling splittings. _J. Chem. Phys._ 2011, _134_ , 054109
  * Richardson 2018 Richardson, J. O. Ring-polymer instanton theory. _Int. Rev. Phys. Chem._ 2018, _37_ , 171–216
  * Richardson et al. 2016 Richardson, J. O.; Pérez, C.; Lobsiger, S.; Reid, A. A.; Temelso, B.; Shields, G. C.; Kisiel, Z.; Wales, D. J.; Pate, B. H.; Althorpe, S. C. Concerted Hydrogen-Bond Breaking by Quantum Tunneling in the Water Hexamer Prism. _Science_ 2016, _351_ , 1310–1313
  * Richardson 2017 Richardson, J. O. Full-and reduced-dimensionality instanton calculations of the tunnelling splitting in the formic acid dimer. _Phys. Chem. Chem. Phys._ 2017, _19_ , 966–970
  * Coleman 1977 Coleman, S. The uses of instantons. 1977; Also in S. Coleman, _Aspects of Symmetry_ , chapter 7, pp. 265–350 (Cambridge University Press, 1985)
  * Jahr et al. 2020 Jahr, E.; Laude, G.; Richardson, J. O. Instanton theory of tunneling in molecules with asymmetric isotopic substitutions. _J. Chem. Phys._ 2020, _153_ , 094101
  * Hammer and Manthe 2011 Hammer, T.; Manthe, U. Intramolecular proton transfer in malonaldehyde: Accurate multilayer multi-configurational time-dependent Hartree calculations. _J. Chem. Phys._ 2011, _134_ , 224305
  * Schröder et al. 2011 Schröder, M.; Gatti, F.; Meyer, H.-D. Theoretical studies of the tunneling splitting of malonaldehyde using the multiconfiguration time-dependent Hartree approach. _J. Chem. Phys._ 2011, _134_ , 234307
  * Mátyus et al. 2016 Mátyus, E.; Wales, D. J.; Althorpe, S. C. Quantum tunneling splittings from path-integral molecular dynamics. _J. Chem. Phys._ 2016, _144_ , 114108
  * Vaillant et al. 2018 Vaillant, C. L.; Wales, D. J.; Althorpe, S. C. Tunneling splittings from path-integral molecular dynamics using a Langevin thermostat. _J. Chem. Phys._ 2018, _148_ , 234102
  * Mil’nikov et al. 2003 Mil’nikov, G. V.; Yagi, K.; Taketsugu, T.; Nakamura, H.; Hirao, K. Tunneling splitting in polyatomic molecules: Application to malonaldehyde. _J. Chem. Phys._ 2003, _119_ , 10–13
  * Sahu et al. 2021 Sahu, N.; Richardson, J. O.; Berger, R. Instanton calculations of tunneling splittings in chiral molecules. _J. Comput. Chem._ 2021, _42_ , 210–221
  * Laude et al. 2018 Laude, G.; Calderini, D.; Tew, D. P.; Richardson, J. O. Ab initio instanton rate theory made efficient using Gaussian process regression. _Faraday Discuss._ 2018, _212_ , 237–258
  * Cooper et al. 2018 Cooper, A. M.; Hallmen, P. P.; Kästner, J. Potential energy surface interpolation with neural networks for instanton rate calculations. _J. Chem. Phys._ 2018, _148_ , 094106
  * McConnell and Kästner 2019 McConnell, S. R.; Kästner, J. Instanton rate constant calculations using interpolated potential energy surfaces in nonredundant, rotationally and translationally invariant coordinates. _J. Comput. Chem._ 2019, _40_ , 866–874
  * Unke et al. 2020 Unke, O. T.; Koner, D.; Patra, S.; Käser, S.; Meuwly, M. High-dimensional potential energy surfaces for molecular simulations: from empiricism to machine learning. _Mach. learn.: sci. technol._ 2020, _1_ , 013001
  * Friesner 2005 Friesner, R. A. Ab initio quantum chemistry: Methodology and applications. _Proc. Natl. Acad. Sci. USA_ 2005, _102_ , 6648–6653
  * Pan and Yang 2009 Pan, S. J.; Yang, Q. A survey on transfer learning. _IEEE Trans. Knowl. Data Eng._ 2009, _22_ , 1345–1359
  * Taylor and Stone 2009 Taylor, M. E.; Stone, P. Transfer learning for reinforcement learning domains: A survey. _J. Mach. Learn. Res._ 2009, _10_ , 1633–1685
  * Smith et al. 2019 Smith, J. S.; Nebgen, B. T.; Zubatyuk, R.; Lubbers, N.; Devereux, C.; Barros, K.; Tretiak, S.; Isayev, O.; Roitberg, A. E. Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning. _Nat. Commun._ 2019, _10_ , 1–8
  * Ramakrishnan et al. 2015 Ramakrishnan, R.; Dral, P.; Rupp, M.; von Lilienfeld, O. A. Big Data meets quantum chemistry approximations: The $\Delta$-machine learning approach. _J. Chem. Theory Comput._ 2015, _11_ , 2087–2096
  * Käser et al. 2020 Käser, S.; Unke, O. T.; Meuwly, M. Reactive dynamics and spectroscopy of hydrogen transfer from neural network-based reactive potential energy surfaces. _New J. Phys._ 2020, _22_ , 055002
  * Käser et al. 2020 Käser, S.; Koner, D.; Christensen, A. S.; von Lilienfeld, O. A.; Meuwly, M. Machine Learning Models of Vibrating H2CO: Comparing Reproducing Kernels, FCHL, and PhysNet. _J. Phys. Chem. A_ 2020, _124_ , 8853–8865
  * Nandi et al. 2021 Nandi, A.; Qu, C.; Houston, P. L.; Conte, R.; Bowman, J. M. $\Delta$-machine learning for potential energy surfaces: A PIP approach to bring a DFT-based PES to CCSD(T) level of theory. _J. Chem. Phys._ 2021, _154_ , 051102
  * Käser et al. 2021 Käser, S.; Boittier, E. D.; Upadhyay, M.; Meuwly, M. Transfer Learning to CCSD(T): Accurate Anharmonic Frequencies from Machine Learning Models. _J. Chem. Theory Comput._ 2021, _17_ , 3687–3699
  * Qu et al. 2021 Qu, C.; Houston, P. L.; Conte, R.; Nandi, A.; Bowman, J. M. Breaking the Coupled Cluster Barrier for Machine-Learned Potentials of Large Molecules: The Case of 15-Atom Acetylacetone. _J. Phys. Chem. Lett._ 2021, _12_ , 4902–4909
  * Käser and Meuwly 2022 Käser, S.; Meuwly, M. Transfer learned potential energy surfaces: accurate anharmonic vibrational dynamics and dissociation energies for the formic acid monomer and dimer. _Phys. Chem. Chem. Phys._ 2022, _24_ , 5269–5281
  * Laude et al. 2020 Laude, G.; Calderini, D.; Welsch, R.; Richardson, J. O. Calculations of quantum tunnelling rates for Muonium reactions with methane, ethane and propane. _Phys. Chem. Chem. Phys._ 2020, _22_ , 16843–16854
  * Wang et al. 2008 Wang, Y.; Braams, B. J.; Bowman, J. M.; Carter, S.; Tew, D. P. Full-dimensional quantum calculations of ground-state tunneling splitting of malonaldehyde using an accurate ab initio potential energy surface. _J. Chem. Phys._ 2008, _128_ , 224314
  * Cvitas and Althorpe 2016 Cvitas, M. T.; Althorpe, S. C. Locating instantons in calculations of tunneling splittings: The test case of malonaldehyde. _J. Chem. Theory Comput._ 2016, _12_ , 787–803
  * Mizukami et al. 2014 Mizukami, W.; Habershon, S.; Tew, D. P. A compact and accurate semi-global potential energy surface for malonaldehyde from constrained least squares regression. _J. Chem. Phys._ 2014, _141_ , 144310
  * Firth et al. 1991 Firth, D.; Beyer, K.; Dvorak, M.; Reeve, S.; Grushow, A.; Leopold, K. Tunable far-infrared spectroscopy of malonaldehyde. _J. Chem. Phys._ 1991, _94_ , 1812–1819
  * Baba et al. 1999 Baba, T.; Tanaka, T.; Morino, I.; Yamada, K. M.; Tanaka, K. Detection of the tunneling-rotation transitions of malonaldehyde in the submillimeter-wave region. _J. Chem. Phys._ 1999, _110_ , 4131–4133
  * Redington et al. 2008 Redington, R. L.; Redington, T. E.; Sams, R. L. Infrared Absorption Spectra in the Hydroxyl Stretching Regions of Gaseous Tropolone OHO Isotopomers. _Z. Phys. Chem._ 2008, _222_ , 1197–1211
  * Murdock et al. 2010 Murdock, D.; Burns, L. A.; Vaccaro, P. H. Vibrational specificity of proton-transfer dynamics in ground-state tropolone. _Phys. Chem. Chem. Phys._ 2010, _12_ , 8285–8299
  * Ortlieb and Havenith 2007 Ortlieb, M.; Havenith, M. Proton transfer in (HCOOH)2: an IR high-resolution spectroscopic study of the antisymmetric C–O stretch. _J. Phys. Chem. A_ 2007, _111_ , 7355–7363
  * Zhang et al. 2017 Zhang, Y.; Li, W.; Luo, W.; Zhu, Y.; Duan, C. High resolution jet-cooled infrared absorption spectra of (HCOOH)2,(HCOOD)2, and HCOOH-HCOOD complexes in 7.2 $\mu$m region. _J. Chem. Phys._ 2017, _146_ , 244306
  * Li et al. 2019 Li, W.; Evangelisti, L.; Gou, Q.; Caminati, W.; Meyer, R. The barrier to proton transfer in the dimer of formic acid: a pure rotational study. _Angew. Chem. Int. Ed. Engl._ 2019, _58_ , 859–865
  * Insausti et al. 2022 Insausti, A.; Ma, J.; Yang, Q.; Xie, F.; Xu, Y. Rotational Spectroscopy of 2-Furoic Acid and Its Dimer: Conformational Distribution and Double Proton Tunneling. _ChemPhysChem_ 2022, _23_ , e202200176
  * Baughcum et al. 1984 Baughcum, S. L.; Smith, Z.; Wilson, E. B.; Duerst, R. W. Microwave spectroscopic study of malonaldehyde. 3. Vibration-rotation interaction and one-dimensional model for proton tunneling. _J. Am. Chem. Soc._ 1984, _106_ , 2260–2265
  * Mil’nikov and Nakamura 2001 Mil’nikov, G. V.; Nakamura, H. Practical implementation of the instanton theory for the ground-state tunneling splitting. _J. Chem. Phys._ 2001, _115_ , 6881–6897
  * Richardson 2018 Richardson, J. O. Perspective: Ring-polymer instanton theory. _J. Chem. Phys._ 2018, _148_ , 200901
  * Richardson et al. 2011 Richardson, J. O.; Althorpe, S. C.; Wales, D. J. Instanton calculations of tunneling splittings for water dimer and trimer. _J. Chem. Phys._ 2011, _135_ , 124109
  * Unke and Meuwly 2019 Unke, O. T.; Meuwly, M. PhysNet: A neural network for predicting energies, forces, dipole moments, and partial charges. _J. Chem. Theory Comput._ 2019, _15_ , 3678–3693
  * Gilmer et al. 2017 Gilmer, J.; Schoenholz, S. S.; Riley, P. F.; Vinyals, O.; Dahl, G. E. Neural message passing for quantum chemistry. Proc. of the 34th Int. Conf. on Machine Learning-Volume 70. 2017; pp 1263–1272
  * Baydin et al. 2017 Baydin, A. G.; Pearlmutter, B. A.; Radul, A. A.; Siskind, J. M. Automatic differentiation in machine learning: a survey. _J. Mach. Learn. Res._ 2017, _18_ , 5595–5637
  * Abadi et al. 2016 Abadi, M.; Barham, P.; Chen, J.; Chen, Z.; Davis, A.; Dean, J.; Devin, M.; Ghemawat, S.; Irving, G.; Isard, M. et al. Tensorflow: A system for large-scale machine learning. 12th USENIX symposium on operating systems Design and Implementation (OSDI 16). 2016; pp 265–283
  * Werner et al. 2019 Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Györffy, W.; Kats, D.; Korona, T.; Lindh, R. et al. MOLPRO, version 2019, a package of ab initio programs. 2019
  * Behler 2015 Behler, J. Constructing high-dimensional neural network potentials: A tutorial review. _Int. J. Quantum. Chem._ 2015, _115_ , 1032–1050
  * Meisner and Kästner 2018 Meisner, J.; Kästner, J. Dual-Level Approach to Instanton Theory. _J. Chem. Theory Comput._ 2018, _14_ , 1865–1872
  * Heller and Richardson 2022 Heller, E. R.; Richardson, J. O. Heavy-Atom Quantum Tunnelling in Spin Crossovers of Nitrenes. _Angew. Chem. Int. Ed._ 2022,
  * Meuwly and Hutson 1999 Meuwly, M.; Hutson, J. Morphing Ab Initio Potentials: A Systematic Study of Ne-HF. _J. Chem. Phys._ 1999, _110_ , 8338–8347
  * Koner and Meuwly 2020 Koner, D.; Meuwly, M. Permutationally invariant, reproducing kernel-based potential energy surfaces for polyatomic molecules: From formaldehyde to acetone. _J. Chem. Theory Comput._ 2020, _16_ , 5474–5484


# An efficient branch-and-cut algorithm for the parallel drone scheduling
traveling salesman problem

Minh Anh Nguyen<EMAIL_ADDRESS>Hai Long Luong  Minh Hoàng
Hà<EMAIL_ADDRESS>Ha-Bang Ban  ORLab, Faculty of Computer
Science, Phenikaa University, Hanoi, Vietnam  School of Information and
Communication Technology, Hanoi University of Science and Technology, Vietnam

###### Abstract

We propose an efficient branch-and-cut algorithm to exactly solve the parallel
drone scheduling traveling salesman problem. Our algorithm can find optimal
solutions for all but two existing instances with up to 229 customers in a
reasonable running time. To make the problem more challenging for future
methods, we introduce two new sets of 120 larger instances with the number of
customers varying from 318 to 783 and test our algorithm and investigate the
performance of state-of-the-art metaheuristics on these instances.

###### keywords:

Parallel drone scheduling; traveling salesman problem; branch and cut;
benchmark instances

verbose,tmargin=2.5cm,bmargin=3cm,lmargin=2.5cm,rmargin=2.5cm,headsep=1.5cm

##  1 Introduction

The goal of this paper is to study the Parallel Drone Scheduling Traveling
Salesman Problem (PDSTSP) defined as follows. We are given a complete
undirected graph $G=(V,E)$ with the set of vertices $V=\\{0,1,...,n,n+1\\}$
and the set of edges $E=\\{(i,j):i,j\in V,i<j\\}$. Vertex 0 represents the
depot while vertex $n+1$ is a copy of the depot vertex and is defined for the
modeling purpose. The remaining vertices represent the set of customers
$V_{c}$. At the depot, there is a truck and a set $D$ of $m$ drones available
for delivery tasks. Denote a subset of drone-eligible customers
$V_{d}\subseteq V_{c}$, whose distance from the depot and parcel weight can be
served by drones. The set of customers which must be served by the truck is
denoted by $V_{t}=V_{c}\setminus V_{d}$. The travel time from node $i$ to node
$j$ of the truck is denoted as $t_{ij}$. The time required by a drone to
service customer $i$ is denoted as $t^{\prime}_{i}$, which is calculated by
doubling the time of flying from the depot to customer $i$. The PDSTSP
consists in constructing a set of routes for the truck and the fleet of drones
such that: (i) every customer vertex is visited exactly once by either truck
or drone; (ii) the truck route starts and ends at the depot; (iii) the drone
performs back-and-forth trips from the depot to serve several customers (one
customer for each single trip); (iv) the completion time required to serve all
customers and return all vehicles to the depot, is minimized.

The PDSTSP has received a significant attention from the literature since its
first introduction in Murray & Chu [5]. The authors propose a Mixed Integer
Linear Programming (MILP) formulation to model the problem and design a simple
heuristic. The latter method divides the customers into two subsets, one for
the truck and one for the drones. Two subproblems: the TSP and the Parallel
Machine Scheduling (PMS) problem, defined on two customer sets, are separately
solved to construct the truck tour and drone tours, respectively. In addition,
an improvement step is performed to reassign customers either to the drone set
or to the truck to better balance the truck and drones’ activities. The
solution approaches are tested on small instances with up to 20 customers.

Many attempts have been done to efficiently solve larger PDSTSP instances.
Mbiadou Saleu, Deroussi, Feillet, Grangeon & Quilliot [3] propose the first
metaheuristic, called iterative two-step heuristic method. The first step
constructs a customer sequence (also called a giant or TSP tour) that visits
all the customers. In the second step, a dynamic programming is used to split
the giant tour into a truck tour and a set of trips for the drones, creating a
complete PDSTSP solution. The method is tested on newly generated instances
with up to 229 customers. Dell’Amico, Montemanni & Novellani [1] later propose
a matheuristic in which the authors use the same idea as found in [3].
However, instead of dynamic programming, a MILP approach is used to handle the
second step. The obtained results show that the algorithm is slower but
provides better solutions than the heuristic of [3].

A nature-inspired algorithm called Hybrid Ant Colony Optimization (HACO) is
proposed in Dinh, Do & Hà [2]. The algorithm reuses the idea of [3] that
represents a PDSTSP solution as a permutation of all customers. However, the
authors propose a new dynamic programming with a controllable amount of states
while exploring wider solution space to split the giant tour into a complete
PDSTSP solution. Several problem-tailored components such as local search and
solution construction also play important roles in improving the performance
of the algorithm. When being tested on benchmark instances from the
literature, the HACO algorithm outperforms the two existing metaheuristics in
terms of both running time and solution quality. The most recent metaheuristic
for the PDSTSP is a Ruin-and-Recreate algorithm called Slack Induction by
String and Sweep Removals (SISSR) introduced in [6]. Exploiting PDSTSP
solution characteristics in an effective manner, the metaheuristic manages to
introduce “sufficient” rooms to a solution via new removal operators during
the ruin phase. It is expected to enhance the possibilities for improving
solutions later in the recreate phase. Also in [6], the authors introduce a
new variant of the PDSTSP called PDSVRP in which they allow multiple trucks
and consider the objective of minimizing the total transportation costs.

Raj, Lee, Lee, Walteros & Murray [9] present a variant of the PDSVRP problem
in which they consider the objective function that minimizes the completion
time all the vehicles come back to the depot. The authors propose three MILP
formulations and develop a branch-and-price (BAP) algorithm and a BAP-based
heuristic to solve the problem. The authors report that their exact methods
can solve to optimality all the small PDSTSP instances of [5] and 3 instances
with 52 customers of [3]. Most recently, Mbiadou Saleu, Deroussi, Feillet,
Grangeon & Quilliot [4] generalize the PDSTSP problem in the same way as the
PDSVRP. The problem is called the Parallel Drone Scheduling Multiple Traveling
Salesman Problem. A hybrid Iterated Local Search is proposed. A MILP
formulation and a simple branch-and-cut approach are also designed to solve
small instances.

As mentioned above, most approaches proposed to solve the PDSTSP belong to the
class of heuristic methods. To the best of our knowledge, the results of exact
algorithms for the problem are quite limited. Only optimal results obtained on
small instances with no more than 20 customers and three instances with 52
customers are provided in [5, 1, 9]. Optimal solutions for 87 over 90
benchmark instances with up to 229 customers created by [3] have never been
reported. The main motivation of this paper is to fill this methodology gap by
proposing a simple, yet efficient exact branch-and-cut (B&C) algorithm for the
PDSTSP. The experimental results on the currently largest problem instances of
[3] show that the developed solution approach can solve all but two instances
in a reasonable running time. To provide more challenging instances to test
solution methods in the future, we introduce a new set of instances with up to
600 customers. We report the results of the state-of-the-art algorithms
obtained on these instances to further investigate their performance.

The remaining of the paper is structured as follows: Section 2 describe our
MILP formulation along with numerous valid inequalities, upon which the
branch-and-cut algorithm is proposed in Section 3. Section 4 reports
experimental results and Section 5 concludes our research.

##  2 Mixed integer linear formulation and valid inequalities

The PDSTSP has been formulated as MILP models in [5, 3, 1, 9]. A common
feature of these models is that they use variables defined on directed arcs to
represent the flow of truck tour. In this study, our model is based on the
truck-flow variables defined on undirected edges to reduce a half of flow
variables. As such, our formulation use four types of variables described as
follows:

  * 1.

$y_{i}$: binary variable equal to 1 if customer $i\in V_{c}$ is served by the
truck, and 0 otherwise;

  * 2.

$x_{ij}$: binary variable equal to $1$ if edge $(i,j)$ is traveled by the
truck ($i,j\in V:i<j$).

  * 3.

$z_{ik}$: binary variable equal to $1$ if customer $i\in V_{d}$ is served by
drone $k\in D$.

  * 4.

$t_{max}$: positive real variable representing the completion time at which
all vehicles return to the depot.

The MILP model for the PDSTSP is formulated as:

minimize $\displaystyle t_{max}$ (1) $\displaystyle\textrm{subject to}\qquad
t_{max}$ $\displaystyle\geq\sum_{i<j}t_{ij}x_{ij}$ (2) $\displaystyle t_{max}$
$\displaystyle\geq\sum_{i\in V_{d}}t^{\prime}_{i}z_{ik}\quad\forall k\in D$
(3) $\displaystyle y_{i}+\sum_{k\in M}z_{ik}$ $\displaystyle=1\quad\forall
i\in V_{d}$ (4) $\displaystyle y_{i}$ $\displaystyle=1\quad\forall i\in
V_{t}\cup\\{0,n+1\\}$ (5) $\displaystyle\sum_{i>0}x_{0i}$ $\displaystyle=1$
(6) $\displaystyle\sum_{i<n+1}x_{i(n+1)}$ $\displaystyle=1$ (7)
$\displaystyle\sum_{i<j}x_{ij}+\sum_{j<k}x_{jk}$
$\displaystyle=2y_{j}\quad\forall j\in V_{c}$ (8)
$\displaystyle\sum_{\begin{subarray}{c}i\in V,j\in V\setminus S\\\ \text{or
}j\in S,i\in V\setminus S\end{subarray}}x_{ij}$ $\displaystyle\geq
2y_{k}\quad\forall k\in S,S\subseteq V_{c},|S|\geq 3$ (9) $\displaystyle
x_{ij}$ $\displaystyle\in\\{0,1\\}\quad\forall i,j\in V:i<j$ (10)
$\displaystyle y_{i}$ $\displaystyle\in\\{0,1\\}\quad\forall i\in V_{c}$ (11)
$\displaystyle z_{ik}$ $\displaystyle\in\\{0,1\\}\quad\forall i\in V_{d},k\in
D$ (12) $\displaystyle t_{max}$ $\displaystyle\geq 0$ (13)

The objective function (1) minimizes the time all vehicles coming back to the
depot. Constraints (2) and (3) define variable $t_{max}$. More precisely, they
provide lower bounds for the objective value in terms of each involved
vehicle’s completion time. Constraints (4) require that each drone-eligible
customer must be served by either the truck or a drone, while constraints (5)
ensure the truck must visit all the customers which cannot be served by the
drones. Constraints (6) and (7) make sure that the truck must travel on only
one edge linking to the depot and its copy, respectively. Constraints (8) are
degree constraints for the truck tour. They imply that if vertex $j$ is
visited by the truck, there must be two edges of the truck tour linking to
$j$. The connectivity of the truck tour is guaranteed by constraints (9),
which are also called the Subtour Elimination Constraints (SECs). They force
the presence of at least two edges between any set $S$ and $V\setminus S$, for
every subset $S$ of $V_{c}$ such that $S$ contains at least one vertex visited
by the truck. Finally, the remaining constraints define the domain of
variables. Note that in this formulation, we allow the case where
$x_{0(n+1)}=0$ since this can happen if the drones serve all the customers and
the truck does not serve any of them.

The formulation above can be enriched by the presence of valid inequalities.
The following constraints follow from the definition of the $x_{ij}$ and
$y_{k}$ variables:

$\displaystyle x_{ij}$ $\displaystyle\leq y_{i}\quad\forall i,j\in V:i<j$ (14)
$\displaystyle x_{ij}$ $\displaystyle\leq y_{j}\quad\forall i,j\in V:i<j$ (15)

The symmetry between the drones can be broken by adding the following
constraints:

$\displaystyle\sum_{i\in V_{d}}d_{i}z_{ik-1}$ $\displaystyle\geq\sum_{i\in
V_{d}}d_{i}z_{ik}\quad\forall k\in M\setminus\\{1\\}$ (16)

Then, the constraints (3) can be replaced by:

$\displaystyle t_{max}$ $\displaystyle\geq\sum_{i\in V_{d}}d_{i}z_{i1}$ (17)

It is possible to derive for the formulation several different forms of the
Subtour Elimination Constraints (SECs). Constraints (9) are equivalent to the
following constraints:

######  Proposition 1.

The constraints

$\displaystyle\sum_{i,j\in S:i<j}x_{ij}\leq\sum_{k\in S\setminus l}y_{k}$
$\displaystyle\forall l\in S,S\subseteq V_{c},|S|\geq 3$ (18)

are valid SECs for the PDSTSP.

Proof. Let $S^{*}$ be the set of vertices visited by the truck in a feasible
PDSTSP solution. The left-hand side of (18) is less than or equal to $|S\cap
S^{*}|-1$, while the right-hand side is equal to $|S\cap S*|-y_{l}$. Because
$y_{l}\leq 1\;\forall l\in S$, constraints (18) are always valid. $\square$

Summing up constraints (18) for every vertex $l\in S$, we have:

######  Proposition 2.

The constraints

$\displaystyle|S|\sum_{i,j\in S:i<j}x_{ij}\leq(|S|-1)\sum_{k\in S}y_{k}$
$\displaystyle\forall S\subseteq V_{c},|S|\geq 3$ (19)

are valid SECs for the PDSTSP.

In a similar manner, the following constraints can be derived directly from
constraints (9):

######  Proposition 3.

The constraints

$\displaystyle|S|\sum_{\begin{subarray}{c}i\in S,j\in V\setminus S\\\ \text{or
}j\in S,i\in V\setminus S\end{subarray}}x_{ij}\geq 2\sum_{k\in S}y_{k}$
$\displaystyle\forall S\subseteq V_{c},|S|\geq 3$ (20)

are valid SECs for the PDSTSP.

Connectivity constraints (9) and (18) can be strengthened if there exists at
least one customer that must be serviced by the truck.

######  Proposition 4.

The constraints

$\displaystyle\sum_{\begin{subarray}{c}i\in S,j\in V\setminus S\\\ \text{or
}j\in S,i\in V\setminus S\end{subarray}}x_{ij}\geq 2$ $\displaystyle\forall
S\subseteq V_{c},S\cap V_{t}\neq\varnothing,|S|\geq 3$ (21)

are valid SECs for the PDSTSP.

Proof. Immediate since in this case the truck tour in the solution must
contain at least one vertex in $S$ and one vertex in $V\setminus S$. $\square$

The equivalent form of constraints (21) can be written as follows:

######  Proposition 5.

The constraints

$\displaystyle\sum_{i,j\in S:i<j}x_{ij}\leq\sum_{k\in S}y_{k}-1$
$\displaystyle\forall S\subseteq V_{c},S\cap V_{t}\neq\varnothing,|S|\geq 3$
(22)

are valid SECs for the PDSTSP.

Proof. Immediate from constraints (18) since in this case there must be a
vertex $l\in S$ such that $y_{l}=1$. $\square$

Several other valid inequalities for the symmetric TSP can be applied directly
to the PDSTSP (see [8]) or can be strengthened. In this study, we use valid
inequalities based on the generation of the 2-matching constraints.

######  Proposition 6.

The constraints

$\displaystyle\sum_{i,j\in H:i<j}x_{ij}+\sum_{(i,j)\in
E^{\prime}}x_{ij}\leq\sum_{j\in H}y_{j}+\frac{1}{2}(|H^{\prime}|-1)$ (23)

for all $H\subseteq V_{c}$, $H^{\prime}\subseteq H$ and $E^{\prime}\subset E$
satisfying: (1) $|H^{\prime}|\geq 3$ and odd; (2) only one endpoint of each
edge in $E^{\prime}$ is in $H^{\prime}$ and each vertex in $H^{\prime}$ must
be an endpoint of a single edge in $E^{\prime}$; (3) there is no edge in
$E^{\prime}$ such that whose both endpoints are in $H$; and (4) two different
edges in $E^{\prime}$ do not share a same endpoint.

are valid inequalities for the PDSTSP.

Proof. Summing up constraints (8) for all $j\in H$, we obtain:

$\displaystyle 2\sum_{i,j\in H:i<j}x_{ij}+\sum_{\begin{subarray}{c}i\in H,j\in
V\setminus H\\\ \text{or }j\in H,i\in V\setminus
H\end{subarray}}x_{ij}=2\sum_{j\in H}y_{j}$ (24)

We add $\sum_{(i,j)\in E^{\prime}}x_{ij}$ to the left-hand side of (24) and
$|H^{\prime}|$ to the right-hand side. Since $x_{ij}\leq 1$ due to (10) and
two different edges in $E^{\prime}$ do not share a same endpoint, equation
(24) becomes the following inequality:

$\displaystyle 2\sum_{i,j\in H:i<j}x_{ij}+2\sum_{(i,j)\in
E^{\prime}}x_{ij}+\sum_{\begin{subarray}{c}i\in H\setminus H^{\prime},j\in
V\setminus H\\\ \text{or }j\in H\setminus H^{\prime},i\in V\setminus
H\end{subarray}}x_{ij}+\sum_{\begin{subarray}{c}i\in H^{\prime},j\in
V\setminus H:(i,j)\notin E^{\prime}\\\ \text{or }j\in H^{\prime},i\in
V\setminus H:(j,i)\notin E^{\prime}\end{subarray}}x_{ij}\leq 2\sum_{j\in
H}y_{j}+|H^{\prime}|$ (25)

Removing the third and fourth terms of the left-hand side from (25) and
dividing by $2$, we obtain:

$\displaystyle\sum_{i,j\in H:i<j}x_{ij}+\sum_{(i,j)\in
E^{\prime}}x_{ij}\leq\sum_{j\in H}y_{j}+\frac{1}{2}|H^{\prime}|$ (26)

As the left-hand side of (26) is integer and $|H^{\prime}|\geq 3$ and
$|H^{\prime}|$ is odd, we can round down the right-hand side, which yields
(23). $\square$

##  3 Branch-and-cut algorithm

The PDSTSP is solved to optimality using a standard branch-and-cut algorithm.
We solve a linear programming (LP) including the objective function (1) and
the constraints (2), (14), (15), (16), (17), and (4)-(8). We then search for
violated constraints of type (9) and (18)-(23), and the detected constraints
are added to the current LP, which is then reoptimized. This process is
repeated until all the constraints are satisfied. If there are fractional
variables, we branch. If all the variables are integer, we explore another
node.

The separation of the connectivity constraints (9) and (18)-(22) is performed
as follows. From a LP solution $(x^{*},y^{*})$ obtained on a node of the
search tree, we build a graph $G^{*}=(V^{*},E^{*})$ where set $V^{*}$ includes
vertex $i$ if the value of variable $y^{*}_{i}$ is positive and set $E^{*}$
contains edge ($i,j$) with the cost equal to $x^{*}_{ij}$. To find set
$S\subset V_{c}$ that violates the connectivity constraints, we first use the
Breadth First Search (BFS) to find all the strongly connected components of
the graph induced by the current LP solution. If the number of components is
at least two, each component is recorded as a set $S$ that we need to find. On
the other hand, if only one component exists and the LP solution is still
fractional, we solve the minimum cut problem using the procedure described in
[8]. The obtained cut separates the vertices of the induced graph into two
subsets, we record one which does not contain the depot and violates one of
SECs.

For each recorded set $S$, if it includes any customers in set $V_{t}$, we
generate SECs (21) and (22) from $S$. Otherwise, SECs (9), (18), (19), (20)
are generated. To limit the number of generated constraints of type (9) and
(18), we only apply them on integer solutions. Finally, we notice that
constraints (9) and (18), (19) and (20), (21) and (22) are pairwise
equivalent. Therefore, it is sufficient to generate only one of them whenever
needed. We select the constraint with the smaller number of nonnegative
coefficients on the left-hand side in the hope of creating the lighter model
that is easier to solve. This idea is also used to propose the Integer Linear
Program based method tackling the TSP in [7].

When the BFS mentioned above cannot find any connected components and the
solution is still fractional, a search for violated 2-matching constraints
(23) is started. One can use the procedure proposed in [8] running in
$\mathcal{O}(n^{4})$ to separate these constraints, which is quite expensive.
In this study, we propose a faster procedure described as follows.

  1. 1.

Let $G^{*}=(V^{*},E^{*})$ be the weighted graph induced by the fractional
solution that needs to be separated. First, four empty sets $\overline{E}$,
$E^{\prime}$, $H$, and $H^{\prime}$ are created. Every vertex $i\in V^{*}$ is
set unmarked.

  2. 2.

We iterate over every edge $(i,j)\in E^{*}$ in the decreasing order of value
$x^{*}_{ij}$. If both vertices $i$ and $j$ are not the endpoints of an edge in
$\overline{E}$, we put edge $(i,j)$ into $\overline{E}$. The aim of this step
is to greedily find set $E^{\prime}$ to possibly maximize the left-hand side
of (23).

  3. 3.

We search for a unmarked vertex $i\in V^{*}$. If we cannot find any unmarked
vertex, the process is terminated. Otherwise, we mark $i$ and go to Step 4.

  4. 4.

Vertex $i$ is put into set $H$. A BFS is triggered to find every unmarked
vertex connecting to $i$ through a path which does not contain any edges of
$\overline{E}$. All the vertices on the path are marked and put into set $H$.

  5. 5.

Set $H^{\prime}$ is defined as $H^{\prime}=H\cap V(\overline{E})$, where
$V(\overline{E})$ is the set of vertices which are the endpoints of all edges
in $\overline{E}$. Then, set $E^{\prime}$ is constructed such that
$E^{\prime}\subset\overline{E}$, $|E^{\prime}|=|H^{\prime}|$, and only one
endpoint of every edge in $E^{\prime}$ is a vertex in $H^{\prime}$. If the
size of $H^{\prime}$ is odd and larger than 2, we check the violation of the
constraint (23) associated with sets $H$, $H^{\prime}$, and $E^{\prime}$. We
then clear sets $H$, $H^{\prime}$ and $E^{\prime}$, and go to Step 3.

It is trivial that the complexity of steps 1, 3, and 5 is
$\mathcal{O}(n^{2})$. In step 2, ordering the edges in $E^{*}$ takes
$\mathcal{O}(n^{2}\log(n^{2}))$ and we must iterate over all edges, the
overall complexity is thus $\mathcal{O}(n^{2}\log(n))$. The BFS procedure in
step 4 runs in $\mathcal{O}(n^{2})$ because a vertex in $V^{*}$ is never
considered more than once. Hence, the complexity of the separation procedure
is $O(n^{2}\log{n})$, which is lower than $O(n^{4})$ of the procedure proposed
in [8].

##  4 Computational experiments

We test the B&C algorithm on the existing benchmark instances introduced in
[5] and [3] with up to 229 customers. Two harder instance sets are also
created to further investigate the performance of our B&C algorithm as well as
the state-of-the-art metaheuristics. For the first set, we follow the same
procedure of [3], but based on larger TSP instances from the TSPLIB library
with up to 783 customers: lin318, pr439, rat575, and rat783. From each
selected TSP instance, we derive 15 PDSTSP instances by modifying the
following parameters: the position of the depot ($dp$), the percentage of
drone-eligible customers ($el$), the drone speed ($sp$), and the number of
drones (#). The second set is generated by the instance generator proposed in
[6]. Each instance is characterized by: the number of customers ($n$), which
is either 400 or 600; the distribution of customer positions ($cd$), which is
either random (r), clustered (c), random clustered (rc); the position of the
depot ($dp$), which is either in the center of the customers (c) or in one
corner of the customers’ region (e); the number of drones ($m$), that is
between 1 and 5. For each value of $n$, 30 instances of this class are
generated.

Our branch-and-cut algorithm is implemented in C++ compiled with GCC 9.3.0 and
built around CPLEX 12.10. All the CPLEX parameters are set to their default
values, except the relative MIP gap tolerance MIPGap set to $10^{-6}$ to reach
the precision of results reported by the existing metaheuristics. All tests
are performed on a desktop operating Kubuntu 20.04, and equipped with an AMD
Ryzen 3700X processor and 32GB RAM. We employ single thread for
reproducibility purposes. The CPLEX is configured to terminate after three
hours of solving an instance. We compare the results of the B&C with two
published state-of-the-art methods, HACO in [2] and SISSR in [6]. Each
metaheuristic is run 10 times on each instance and the best solution is
recorded. We set each run of the metaheuristics to stop after 300 seconds on
instances with less than 300 customers, and 600 seconds on the larger
instances. The best solutions found by the heuristics are also used to provide
good upper bounds for the exact method. The detailed benchmark instances and
results of our experiments can be found at http://orlab.com.vn/home/download.

Table 1: Experimental results of the B&C algorithm Instance | # Opt | # Imp | gap | time
---|---|---|---|---
Name | $n$ | # Ins | Source |  |  | (%) | (s)
| 10 | 360 | [5] | 360 | 0 | 0.00 | 0.01
| 20 | 360 | [5] | 360 | 0 | 0.00 | 0.06
att48 | 48 | 15 | [3] | 15 | 0 | 0.00 | 0.89
berlin52 | 52 | 15 | [3] | 15 | 0 | 0.00 | 0.71
eil101 | 101 | 15 | [3] | 15 | 0 | 0.00 | 12.01
gr120 | 120 | 15 | [3] | 14 | 0 | 0.00 | 740.17
pr152 | 152 | 15 | [3] | 15 | 0 | 0.00 | 72.95
gr229 | 229 | 15 | [3] | 14 | 0 | 0.00 | 1680.39
lin318 | 318 | 15 | this paper | 10 | 4 | 0.01 | 4998.34
pr439 | 439 | 15 | this paper | 14 | 7 | 0.00 | 3241.04
rat575 | 575 | 15 | this paper | 7 | 12 | 0.20 | 7885.87
rat783 | 783 | 15 | this paper | 1 | 5 | 0.42 | 10704.08
| 400 | 30 | this paper | 16 | 20 | 0.05 | 5940.12
| 600 | 30 | this paper | 3 | 16 | 0.23 | 10056.52

The obtained results are summarized in Table 1 where Column ‘Name’ shows the
name of the instance class (if any), Column ‘$n$’ reports the number of
customers, Column ‘# Ins’ represents the number of instances in the class, and
Column ‘Source’ shows the reference where the instances are generated. The
performance of the algorithm can be observed via the number of successfully
solved instances (Column ‘# Opt’), the number of instances in which the B&C
algorithm improves the results of the state-of-the-art metaheuristics (Column
‘# Imp’), the averaged gap values returned by CPLEX (Column ‘gap’), and the
running time on average of the algorithm in seconds (Column ‘time’).

As can be seen in the table result, the B&C algorithm can solve to optimality
all the small instances of [5] in very short running times (less than 1 second
in all cases). It is worth recalling that, the best branch-and-cut algorithm
in the literature, which is proposed in [9], can also successfully solve all
of these instances. However, they took 0.86 seconds on average to find the
optimal solution on the 20-customers instances compared to our method which
took 0.06 seconds on average. To solve the hardest instance, they took 47.95
seconds while we took 0.68 seconds, which is much faster. Although their
algorithm is run on an Intel i7-6700 processor, which has 1.2 times lower
single-thread performance compared to our computer (Passmark Software 2021),
our algorithm is clearly faster.

For the largest existing instances introduced in [3], our algorithm provides
optimal solutions for 88 over 90 instances. For the two instances that we
cannot solve to optimality, the gaps are quite small (less than 0.01%). In
[9], the authors also test their exact methods on a subset of 18 instances
selected from 90 instances of [3], and only three optimal solutions for the
berlin52-based instances are found. This again confirms the good performance
of our B&C.

For the first class of the new instances, our algorithm is capable to
optimally solve all pr439-based instances but one, two third lin318-based
instances, and nearly a half of rat575-based instances. On the largest
instances with 783 customers, our algorithm can successfully solve only one
instance. The second class of the new instances tends to be more difficult for
the algorithm since we can solve to optimality only 16 over 30 400-customer
instances and 3 over 30 600-customer instances. All the valid cuts of the
algorithm are activated during the search and constraints (9), (21), (22), and
(23) are created more extensively. Regarding the performance of the two
metaheuristics, we observe that SISSR dominates HACO in terms of solution
quality. The performance of both heuristic methods tends to decrease markedly
when the instance size increases since the exact method can improve 64
solutions of the metaheuristics, all improved cases appear on the new
instances. Therefore, there is a huge room for designing new metaheuristics in
future researches to efficiently handle large PDSTSP instances.

##  5 Conclusions

In this research, we propose a simple, yet efficient branch-and-cut algorithm
to exactly solve the PDSTSP. The problem is first formulated as a mixed
integer linear program with truck-flow variables defined on undirected edges,
not on directed arcs as in existing models. The formulation is then
strengthened by valid inequalities and the branch-and-cut algorithm is
developed. The results obtained on the existing benchmark instances show a
good performance of the algorithm when it can successfully solve 88 over 90
instances. We also introduce two new sets of larger instances with up to 783
customers. When being tested on these instances, our exact algorithm can
successfully solve 51 over 120 instances. The results also allow to
investigate the performance of the state-of-the-art metaheuristics. The SISSIR
algorithm is the best metaheuristic and performs well on the instances of the
literature since it can find all optimal solutions provided by the B&C.
However, to efficiently solve larger instances including ones introduced in
this paper, it is necessary to design a new metaheuristic in future
researches.

## References

  * Dell’Amico et al. [2020] Dell’Amico, M., Montemanni, R., & Novellani, S. (2020).  Matheuristic algorithms for the parallel drone scheduling traveling salesman problem.  Annals of Operations Research, 289, 211–226.
  * Dinh et al. [2021] Dinh, Q. T., Do, D. D., & Hà, M. H. (2021).  Ants can solve the parallel drone scheduling traveling salesman problem.  In Proceedings of the Genetic and Evolutionary Computation Conference (p. 14–21).  New York, NY, USA: Association for Computing Machinery.
  * Mbiadou Saleu et al. [2018] Mbiadou Saleu, R. G., Deroussi, L., Feillet, D., Grangeon, N., & Quilliot, A. (2018).  An iterative two-step heuristic for the parallel drone scheduling traveling salesman problem.  Networks, 72, 459–474.
  * Mbiadou Saleu et al. [2021] Mbiadou Saleu, R. G., Deroussi, L., Feillet, D., Grangeon, N., & Quilliot, A. (2021).  The parallel drone scheduling problem with multiple drones and vehicles.  European Journal of Operational Research, forthcoming.
  * Murray & Chu [2015] Murray, C. C., & Chu, A. G. (2015).  The flying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery.  Transportation Research Part C: Emerging Technologies, 54, 86–109.
  * Nguyen et al. [2021] Nguyen, M. A., Pham, T.-H. G., Hà, M. H., & Pham, M.-T. (2021).  The min-cost parallel drone scheduling vehicle routing problem.  European Journal of Operational Research, In Press.
  * Oswin et al. [2017] Oswin, A., Fischer, A., Fischer, F., Meier, J. F., Pferschy, U., Pilz, A., & Staněk, R. (2017).  Minimization and maximization versions of the quadratic travelling salesman problem.  Optimization, 66, 521–546.
  * Padberg & Rinaldi [1990] Padberg, M., & Rinaldi, G. (1990).  Facet identification for the symmetric traveling salesman polytope.  Math. Program., 47, 219–257.
  * Raj et al. [2021] Raj, R., Lee, D., Lee, S., Walteros, J., & Murray, C. (2021).  A branch-and-price approach for the parallel drone scheduling vehicle routing problem.  URL: https://ssrn.com/abstract=3879710.