In the present work, we studied HD+ (dpe) molecular ion in the framework of variational method without assuming the Born-Oppenheimer approximation. The non-adiabatic wave function was expanded in terms of explicitly correlated Gaussian basis functions. An algorithm for calculating the nucleus-nucleus correlation function (i.e., the probability density of one nucleus in the reference frame where the other one is at the origin) was derived, implemented, and used to depict all bound states of HD+ with zero rotational energy. (C) 2005 Elsevier B.V. All rights reserved.