We studied the orientational order-disorder transition of crystals made up of near-spherical Lennard-Jones dumbbells, of which reduced bond lengths are 0.225, 0.250 and 0.275. Various techniques of Monte Carlo (MC) simulations are used to calculate the chemical potentials of ordered and disordered crystals, and thereby to predict orderdisorder phase transition. First, we performed NPT MC simulations to determine crystal structure, equilibrium positions and orientations of the molecules. We then calculated the free energies of the crystals using the expanded ensemble MC simulations combined with the Einstein-molecule method and the thermodynamic integration method. The solid-solid phase equilibrium is determined from the free energy profiles of the individual phases by equating the chemical potential. The predictions of phase transition obtained from the conventional NPT MC simulation and the free energy simulation were in excellent agreement with each other, which confirms the validity of the present method of calculating the chemical potential of crystal. In addition, the Gibbs-Duhem integration was performed to obtain a complete coexistence curve between the two crystal phases. Orientational probability distributions of molecular axes were analyzed to find the characteristic behavior of rotational motion of molecule in the crystal. At sufficiently low temperature, flipping rotation of molecule in the ordered crystal is suppressed. In contrast, the flipping rotation occurs at higher temperature close to the transition while orientationally ordered structure is still maintained. In the freeenergy calculation, such a unique rotational behavior requires to use a suitable form of external rotational potential with proper symmetry number. The present study demonstrates how one can judiciously choose a correct simulation scheme for the calculation of chemical potentials of molecular crystals.