Book excerpt: There are many tools available to measure the thermophysical properties of compounds. Experimental measurements have been evolving for many years and are very accurate at determining the properties of most compounds. However, many of the measurements are unreliable when the compound of interest is thermally unstable. Throughout the years molecular simulation techniques have been developed to understand the thermophysical properties of thermally unstable compounds. There are primarily two methods to study Vapor-Liquid Equilibrium by molecular simulation Gibbs Ensemble Monte Carlo and Molecular Dynamics. MD is a technique that allows one to simulate the vapor and the liquid in the same simulation cell. The advantage to having the vapor and liquid in the same simulation cell is that an interface forms and properties not available by GEMC can be investigated. However, the inclusion of the interface complicates the determination of the phase densities. There are two methods available in the literature to determine the phase densities from a two-phase MD simulation. The first utilizes a hyperbolic tangent function to fit the density profile across the axis normal to the interface. The second method calculates the average of a local property spatially and then determines the resulting distribution function. The distribution function is used to determine the phases from user defined phase cut-offs. These methods only work well far from the critical point and have many adjustable parameters. These adjustable parameters make it difficult to reliably obtain accurate results. This lack of reliability is one of the main driving forces behind this dissertation. In order to correct the limitations of previous methods, a new technique is presented and tested against three cases. The new technique utilizes Voronoi tessellations to calculate the volume of every molecule in the simulation cell. The molecular volumes generated can be interpreted by simple statistical parameters such as the mean and variance to determine the density on the two phase envelope. In this dissertation a new method is presented and applied to three test cases, a simple fluid, and two polyatomic cases.