  Chemical bonding is a ubiquitous concept in chemistry and it provides a
common basis for experimental and theoretical chemists to explain and predict
the structure, stability and reactivity of chemical species. Among others, the
Energy Decomposition Analysis (EDA, also known as the Extended Transition State
method) in combination with Natural Orbitals for Chemical Valence (EDA-NOCV) is
a very powerful tool for the analysis of the chemical bonds based on a charge
and energy decomposition scheme within a common theoretical framework. While
the approach has been applied in a variety of chemical contexts, the current
implementations of the EDA-NOCV scheme include relativistic effects only at
scalar level, so simply neglecting the spin-orbit coupling effects and de facto
limiting its applicability. In this work, we extend the EDA-NOCV method to the
relativistic four-component Dirac-Kohn-Sham theory that variationally accounts
for spin-orbit coupling. Its correctness and numerical stability have been
demonstrated in the case of simple molecular systems, where the relativistic
effects play a negligible role, by comparison with the implementation available
in the ADF modelling suite (using the non-relativistic Hamiltonian and the
scalar ZORA approximation). As an illustrative example we analyse the
metal-ethylene coordination bond in the group 6-element series
(CO)$_5$TM-C$_2$H$_4$, with TM =Cr, Mo, W, Sg, where relativistic effects are
likely to play an increasingly important role as one moves down the group. The
method provides a clear measure (also in combination with the CD analysis) of
the donation and back-donation components in coordination bonds, even when
relativistic effects, including spin-orbit coupling, are crucial for
understanding the chemical bond involving heavy and superheavy atoms.
