In this paper, we describe a parallel molecular dynamics simulation scheme in NPT ensemble which can be used for studying molecular systems with bond constraints. As an example of application, we choose to investigate general anesthesia whose molecular mechanism is largely unknown. In doing that, we performed two simulations. One is pure lipid bilayer simulation and the other is lipid bilayer simulation with 50 mol% halothane (which is a well-known inhaled anesthetic). From the simulations, we have found that d-spacing increases and area per molecule decreases in the presence of the halothane molecules. We have also found that the halothane molecules are distributed in the hydrophobic tail region, but near headgroup region.