  We report an implementation of the nuclear magnetic resonance (NMR) shielding
($\sigma$), isotope-independent indirect spin-spin coupling ($K$) and the
magnetizability ($\xi$) tensors in the frozen density embedding (FDE) scheme
using the four-component (4c) relativistic Dirac--Coulomb (DC) Hamiltonian and
the non-collinear spin density functional theory (SDFT). The formalism takes
into account the magnetic balance between the large and the small components of
molecular spinors and assures the gauge-origin independence of NMR shielding
and magnetizability results. This implementation has been applied to
hydrogen-bonded HXH$\cdots$OH$_2$ complexes (X = Se, Te, Po) and compared with
the supermolecular calculations and with the approach based on the integration
of the magnetically induced current density vector. A comparison with the
approximate Zeroth-Order Regular Approximation (ZORA) Hamiltonian indicates
non-negligible differences in $\sigma$ and $K$ in the HPoH$\cdots$OH$_2$
complex, and calls for a thourough comparison of ZORA and DC in the description
of environment effects on NMR parameters for molecular systems with heavy
elements.
