The theoretical treatment in Paper I [D. W. Miller and S. A. Adelman, J. Chem. Phys. 117, 2672, (2002), preceding paper] of the vibrational energy relaxation (VER) of low-frequency, large mass dihalogen solutes is extended to the VER of the high-frequency, small mass molecular hydrogen solutes H-2 and D-2 in a Lennard-Jones argon-like solvent. As in Paper I, values of the relaxation times T-1 predicted by the theory are tested against molecular dynamics (MD) results and are found to be of semiquantitative accuracy. To start, it is noted that standard Lennard-Jones site-site potentials derived from macroscopic data can be very inaccurate in the steep repulsive slope region crucial for T-1. Thus, the H-Ar Lennard-Jones diameter sigma(UV) is not taken from literature values but rather is chosen as sigma(UV)=1.39 Angstrom, the value needed to make the theory reproduce the experimental H-2/Ar gas phase VER rate constant. Next, by MD simulation it is shown that the vibrational coordinate fluctuating force autocorrelation function <(F) over tilde (t)(F) over tilde >(0) of Paper I decays roughly an order of magnitude more rapidly for the molecular hydrogen solutions than for the dihalogen solutions. This result implies a relatively slow decay for the molecular hydrogen friction kernels beta(omega)=(k(B)T)(-1)integral(0)(infinity)<(F) over tilde (t)(F) over tilde >(0) cos omega tdt, yielding for the H-2/Ar and D-2/Ar systems at T=150 K physical millisecond values for T-1=beta(-1)(omega(l)) despite the high liquid phase vibrational frequencies omega(l) of H-2 and D-2. The rapid decay of <(F) over tilde (t)(F) over tilde >(0) is due to both the steepness of the repulsive slope of the H-Ar potential and the small masses of H and D. Thus, the small value chosen for sigma(UV) is needed to avoid unphysically long T-1's. Next, an analytical treatment of the H-2/D-2 isotope effect on T-1, based on the theory, is found to predict that the H-2/Ar and D-2/Ar T-1's are close in value due to the compensating effects of lower omega(l) but slower decay of <(F) over tilde (t)(F) over tilde >(0) for D-2/Ar, a result in qualitative agreement with experiments. Applying the theory to numerically study the isothermal rho dependencies of the VER rate constant k(T,rho)=T-1(-1) at 150 K reveals that for both H-2/Ar and D-2/Ar, as for the solutions of Paper I, k(T,rho) can be factorized as in the isolated binary collision (IBC) model. Moreover, the molecular theory and IBC rate isotherms differ only slightly for both solutions, a result interpreted in terms of the form of the H-Ar pair correlation function. The theoretical and experimental rate isotherms at 150 K are then compared. Agreement is very good for the H-2/Ar solution, but for the D-2/Ar solution the theoretical rates are about four times too large. Finally, the isochoric T dependencies of k(T,rho) in the range 200-1000 K are found for both solutions to conform to an Arrhenius rate law.