  We assess the predictive accuracy of perturbation theory based estimates of
changes in covalent bonding due to linear alchemical interpolations among
molecules. We have investigated $\sigma$ bonding to hydrogen, as well as
$\sigma$ and $\pi$ bonding between main-group elements, occurring in small sets
of iso-valence-electronic molecular species with elements drawn from second to
fourth rows in the $p$-block of the periodic table. Numerical evidence suggests
that first order estimates of covalent bonding potentials can achieve chemical
accuracy if (i) the alchemical interpolation is vertical (fixed geometry), (ii)
involves molecules containing elements in the third and fourth row of the
periodic table, and (iii) a reference geometry is optimized. In this case,
changes in the bonding potential become near-linear in coupling parameter,
resulting in analytical predictions with very high accuracy ($\sim$1 kcal/mol).
Second order estimates deteriorate the prediction. If initial and final
molecules differ not only in composition but also in geometry, all estimates
become substantially worse, with second order being slightly more accurate than
first order. The independent particle approximation to the second order
perturbation performs poorly when compared to the coupled perturbed or finite
difference approach. Taylor series expansions up to fourth order of the
potential energy curve of highly symmetric systems indicate a finite radius of
convergence, as illustrated for the alchemical stretching of H$_2^+$. Numerical
results are presented for covalent bonds to hydrogen in 12 molecules with 8
valence electrons; (ii) main-group single bonds in 9 molecules with 14 valence
electrons; (iii) main-group double bonds in 9 molecules with 12 valence
electrons; (iv) main-group triple bonds in 9 molecules with 10 valence
electrons; (v) H$_2^+$ single bond with 1 electron.
