The authors calculate the spectral density functions required for the theoretical determination of nuclear spin relaxation rates due to magnetic dipolar coupling between diffusing hydrogen spins in BCC metals. Both the like-spin (hydrogen-hydrogen) and unlike-spin (metal-hydrogen) contributions are determined using the Monte Carlo method of Faux, Ross and Sholl (1986). Results are obtained for hydrogen-to-metal ratios of 0.12, 0.3 and 0.6 for the simple hopping model (no multiple occupancy of sites) and for a multiple site-blocking model in which the hydrogen spins block all sites as far as the second or third neighbour. For the simple hopping model, the results are in very good agreement with the multiple-scattering theory of Sankey and Fedders (1980), in the high- and low-temperature limits. For the multiple site-blocking models, it is found that the BPP model (due to Bloembergen, Purcell and Pound, 1948) and Torrey model (1953), can differ significantly from the Monte Carlo result, particularly at higher concentrations. The results obtained for D/T1, where D is the tracer diffusion coefficient and T1 is the spin-lattice relaxation time, are compared to the experimental result on NbH0.6. The Monte Carlo method gives a value which is 2/3 of the experimental value if blocking to either the second or the third neighbour is assumed. Agreement with experiment may be obtained if it is assumed that the hydrogen diffuses by a combination of jumps to nearest-neighbour and second-nearest-neighbour sites.