  We extensively develop a method of implementing mean-field calculations for
deformed nuclei, using the Gaussian expansion method (GEM). This GEM algorithm
has the following advantages: (i) it can efficiently describe the
energy-dependent asymptotics of the wave functions at large $r$, (ii) it is
applicable to various effective interactions including those with finite
ranges, and (iii) the basis parameters are insensitive to nuclide, thereby many
nuclei in wide mass range can be handled by a single set of bases. Superposing
the spherical GEM bases with feasible truncation for the orbital angular
momentum of the single-particle bases, we obtain deformed single-particle
wave-functions to reasonable precision. We apply the new algorithm to the
Hartree-Fock and the Hartree-Fock-Bogolyubov calculations of Mg nuclei with the
Gogny interaction, by which neck structure of a deformed neutron halo is
suggested for $^{40}$Mg.
