The two-phase MD technique employed in this work determines the liquid and vapor phase densities from a histogram of molecular densities within phase clusters in the simulation cell using a new Monte Carlo (MC) sampling method. These equilibrium densities are then fitted in conjunction with known critical-point scaling laws to obtain the critical temperature, and the critical density. This MC post-processing method was found to be more easily implemented in code, and it is efficient and easily applied to complex, structured molecules. This method has been successfully applied and benchmarked for a simple Lennard-Jones (LJ) fluid and a structured molecule, propane. Various degrees of internal flexibility in the propane models showed little effect on the coexisting densities far from critical point, but internal flexibility (angle bending and bond vibrations) seemed to affect the saturated liquid densities in the near-critical region, changing the critical temperature by approximately 20 K. Shorter cutoffs were also found to affect the phase dome and the location of the critical point. The developed MD+MC method was then used to test the efficacy of two all-atom, site-site pair potential models (with and without point charges) developed solely from the energy landscape obtained from high-level ab initio pair interactions for the first time. Both models produced equivalent phase domes and critical loci. The model's critical temperature for methanol is 77 K too high while that for 1-propanol is 80 K too low, but the critical densities are in good agreement. These differences are likely attributable to the lack of multi-body interactions in the true pair potential models used here. Lastly, the transferability of the ab initio potential model was evaluated by applying it to 1-pentanol. An attempt has been made to separate the errors due to transferability of the potential model from errors due to the use of a true-pair potential. The results suggested a good level of transferability for the site-site model. The lack of multi-body effects appears to be dominant weakness in using the generalized ab initio potential model for determination of the phase dome and critical properties of larger alcohols.