Monte Carlo computer simulation techniques are used to study the orientational phase transition in the three-dimensional Lebwohl-Lasher model which couples molecular rotors placed on a cubic lattice by the potential P₂(cosθ_ij). The orientational transition is a model of the nematic-isotropic phase transition in liquid crystals. The simulations involve the determination of energy and order parameter distribution functions which permit free energy functions to be derived. The data have been analysed by finite size scaling methods to reveal the nature of the phase transition which is found to be weakly first order with stability limits of the nematic and isotropic phases being extremely close to the equilibrium transition temperature, in good agreement with experimental studies of room temperature nematogens. Results are reported for the specific heat, the axial and biaxial susceptibilities, as well as the enthalpy and nematic order parameter discontinuity at the transition. It is shown that the inclusion of a term p₄(cosθ_ij) in the potential enhances the first-order character of the transition and displaces the stability limits further from the equilibrium transition temperature. The director fluctuations have been analysed, and it is found that, whereas the fluctuations in the isotropic phase follow ordinary Brownian motion, the fluctuations in the nematic phase correspond to fractional Brownian motion. By introducing a symmetry- breaking field, –h² cosθ²_i, a field-induced crossover between fractional and normal Brownian motion is observed in agreement with recent neutron scattering studies on deuterated p-azoxyanisole.