The high-quality Fermi LAT observations of gamma-ray pulsars have opened a new window to understanding the generation mechanisms of high-energy emission from these systems. The high statistics allow for careful modeling of the light curve features as well as for phase-resolved spectral modeling. We model the LAT light curves of four bright LAT pulsars using simulated high-energy light curves. The model light curves and phase-dependent radii of curvature are generated using geometrical representations of the outer gap and slot gap/two-pole caustic emission models, within the context of both the vacuum retarded dipole and force-free magnetosphere models. These simulated light curves are compared with observed LAT light curves via maximum likelihood using a Markov Chain Monte Carlo method to explore the phase space of fitted parameters such as magnetic inclination, viewing angle, maximum emission radius and gap width. We find that the observed light curves can be fit within both the vacuum dipole and force-free fields, but the force-free magnetosphere produces phase lags between the gamma-ray and radio peaks larger than those observed. We have also used the measured phase-dependent spectral cutoff energies to estimate the accelerating parallel electric field dependence on emission radius for each pulsar, under the assumptions that the high-energy emission is dominated by curvature radiation and the geometry (radius of emission and minimum radius of curvature of the magnetic field lines) is determined by the best fitting light curves for each model.