During major disruptions, an induced loop voltage accelerates runaway electrons (REs) towards high energy, being in the order of 1–100 MeV in present tokamaks and ITER. The stochastization mechanisms of such high-energy RE drift orbits are investigated by three-dimensional (3D) orbit following in tokamak plasmas. Drift resonance is shown to play an important role in determining the onset of stochastic drift orbits for different electron energies, particularly in cases with low-order perturbations that have radially global eigenfunctions of the scale of the plasma minor radius. The drift resonance due to the coupling between the cross-field drift motion with radially global modes yields a secondary island structure in the RE drift orbit, where the width of the secondary drift islands shows a square-root dependence on the relativistic gamma factor γ. Only for highly relativistic REs (γ ≫ 1), the widths of secondary drift islands are comparable with those of magnetic islands due to the primary resonance, thus the stochastic threshold becoming sensitive to the RE energy. Because of poloidal asymmetry due to toroidicity, the threshold becomes sensitive not only to the relative amplitude but also to the phase difference between the modes. In this paper, some examples of 3D orbit-following calculations are presented for analytic models of magnetic perturbations with multiple toroidal mode numbers, for both possibilities that the drift resonance enhances and suppresses the stochastization being illustrated.