Hybrid Monte Carlo (HMC) is a rigorous sampling method that uses molecular dynamics (MD) as a global Monte Carlo move to search the phase space of molecules. However, the acceptance rate of HMC, and thus the efficiency of the method, decays exponentially with increased system size or timestep. The Shadow Hybrid Monte Carlo (SHMC) method was previously introduced to overcome this performance degradation by sampling from the shadow Hamiltonian in lieu of the actual Hamiltonian. A numerical integrator for a Hamiltonian system of ordinary differential equations exactly solves a shadow Hamiltonian if and only if the integrator is symplectic. High order approximations to the shadow Hamiltonian are efficiently computed, thus, allowing significant acceptance of the MD step in SHMC. In order to remove the bias of sampling from the shadow Hamiltonian, a reweighting step is performed on the collected data. The increase in efficiency of SHMC comes at the expense of the variance of reweighted observables. In this dissertation, I discuss the development of an analytical framework for the analysis of shadow hybrid Monte Carlo methods. From this analysis, I show that the performance of SHMC is limited by the need to generate momenta for the MD step from a non-separable shadow Hamiltonian. In addition, I present a method for determining the optimal parameterization of SHMC. This includes the ability to ascertain the timestep which will allow for maximum efficiency of the method. The reweighting step usually causes an increase in the variance of observed data, but it is now possible to ascertain the magnitude of the variance beforehand. Finally, I introduce a new algorithm, Separable Shadow Hybrid Monte Carlo (S2HMC), which is based on a separable formulation of the shadow Hamiltonian. S2HMC is able to quickly generate momenta from the correct distribution, thereby overcoming a limitation of SHMC. S2HMC has the added advantage that no additional parameters are required. Through analysis and numerical experiments I show that S2HMC is at least twice more efficient than HMC for all systems larger than a few thousand atoms.