Orientational ordering of C-60 molecules within monolayer and multilayer islands is a regularly observed phenomenon in scanning tunneling microscopy (STM) studies. Here we simulate the orientational ordering seen in STM images via a novel combination of Monte Carlo and Huckel theory methods and compare to experimental data. A measure of the repulsive interaction energy between two adjacent C-60 molecules is precalculated by estimating and processing the electron density distribution between them. Many combinations of molecular orientations are considered to encompass all the details of the molecular orbitals. Precalculated intermolecular interaction energies are inputted into a simulated C-60 island. Here, the center position of each molecule is fixed, but the molecules are allowed to rotate freely around their centers. A minimum in the total island free energy is sought by sequentially picking molecules at random and rotating them according to their neighbors. Results show significant correlation with experimentally observed features in both mono- and multilayered islands on a variety of different substrates.