Semicrystalline polymers are an important class of materials for their prevalence in today's markets and their desirable properties. These properties depend on the early stages of the polymer crystallization process where a crystal nucleates from the polymer melt. This nucleation process is conventionally understood via an extension of Classical Nucleation Theory to polymers (CNTP). However, recent experimental and simulation evidence points to nucleation mechanisms that do not agree with the predictions of CNTP. Specifically, these experiments suggest a previously unrecognized role of nematic phases in mediating the melt"“crystal transtion. To explain these observations, several new theories of nucleation alternate to CNTP have emerged in the literature, all of which suggest specific modifications to the free energy landscape (FEL) near-equilibrium. To address these theoretical controversies, this dissertation aimed to study the equilibrium phase behavior of polymers via Monte Carlo (MC) simulations. Simulating equilibrium phase behavior of polymer melts is not a trivial task due to the large free energy barriers involved. Throughout this research, we employed a combination of strategies to speed up these molecular simulations. First, we employed a domain decomposition to divide the simulation box into multiple independent simulations that execute independent MC trajectories in parallel. The novel GPU-accelerated MC algorithm successfully and accurately simulated the phase behavior of bead spring chains. Additionally, it sped up MC simulations of Lennard Jones chains by up to 10 times. In its current form, the GPU-accelerated algorithm did not achieve significant speedups to improve outcomes of simulating large polymer melts with detailed potentials. We recommended various strategies to improving the current algorithm. This reality motivated the use of biased MC simulations to study the phase behavior of polymers more expediently without the need for GPU acceleration. Specifically, the latter part of the Dissertation employed Wang Landau MC (WLMC) simulations to build phase diagrams and expanded ensemble density of states (EXEDOS) simulations to construct FELs. Phase diagrams from WLMC simulations divided volume-temperature space into melt, nematic and crystal phases. Then, FELs from EXEDOS simulations at equilibrium provided direct access to the relative stability and minimum free energy paths between coexistant states. By employing a two-dimensional EXEDOS sampling in both crystal and nematic order for hard bead semiflexible oligomers with a stepwise bending stiffness, we built FELs that show that the crystalline transition cooperatively and simultaneously formed crystal and nematic order. This nucleation mechanism was not in agreement with predictions from CNTP or newer theoretical formulations. To investigate the sensitivity of the phase behavior to the employed polymer model, we then employed WLMC simulations to build phase diagrams for a number of different polymer models to ascertain their impact on the resulting nucleation mechanism. We found that the phase behavior was sensitive to the form of the bending stiffness potential used. Chains with a stepwise bending stiffness yielded the previously mentioned cooperative and simultaneous crystal and nematic ordering. In contrast, chains with a harmonic bending stiffness potential crystallized via a two-step nucleation process, first forming a nematic phase that nucleates the crystal. The latter nucleation mechanism was in line with predictions from new theories of nucleation that incorporate the nematic phase as a precursor. Furthermore, we found that it is important to correct for excluded volume differences when comparing chains with soft and hard beads or chains with differing bending stiffnesses.