Explicitly correlated all-particle Gaussian functions with shifted centers (ECGs) are implemented within the earlier proposed effective variational non-Born-Oppenheimer method for calculating bound states of molecular systems in magnetic field (Adamowicz et al., 2015). The Hamiltonian used in the calculations is obtained by subtracting the operator representing the kinetic energy of the center-of-mass motion from the total laboratory-frame Hamiltonian. Test ECG calculations are performed for the HD molecule. (C) 2017 Elsevier B.V. All rights reserved.