  We develop the isospin-invariant Skyrme-EDF method by considering local
densities in all possible isospin channels and proton-neutron (p-n) mixing
terms as mandated by the isospin symmetry. The EDF employed has the most
general form that depends quadratically on the isoscalar and isovector
densities. We test and benchmark the resulting p-n EDF approach, and study the
general properties of the new scheme by means of the cranking in the isospin
space. We extend the existing axial DFT solver HFBTHO to the case of
isospin-invariant EDF approach with all possible p-n mixing terms. Explicit
expressions have been derived for all the densities and potentials that appear
in the isospin representation. In practical tests, we consider the Skyrme EDF
SkM* and, as a first application, concentrate on Hartree-Fock aspects of the
problem, i.e., pairing has been disregarded. Calculations have been performed
for the (A=78, T~11), (A=40, T~8), and (A=48, T~4) isobaric analog chains.
Isospin structure of self-consistent p-n mixed solutions has been investigated
with and without the Coulomb interaction, which is the sole source of isospin
symmetry breaking in our approach. The extended axial HFBTHO solver has been
benchmarked against the symmetry-unrestricted HFODD code for deformed and
spherical states. We developed and tested a general isospin-invariant
Skyrme-EDF framework. The new approach permits spin-isospin densities that may
give rise to, hitherto, unexplored modes in the excitation spectrum. The new
formalism has been tested in the Hartree-Fock limit. A systematic comparison
between HFODD and HFBTHO results show a maximum deviation of about 10 keV on
the total binding energy for deformed nuclei when the Coulomb term is included.
Without this term, the results of both solvers agree down to a ~10 eV level.
