In the present work, isothermal-isobaric (NPT) and canonical (NVT) ensemble molecular dynamics (MD) simulations have been used to calculate the binary mixture properties (excess molar heats and volumes) of 1-ethyl-3-methylimidazolium methyl sulfate [EMIM][MeSO4] with six polar covalent molecules (acetic acid, acetone, chloroform, dimethyl sulfoxide, isopropyl alcohol and methanol) over the entire composition range. Excess partial molar properties as a function of mole fraction of co-solvent were calculated by using a Redlich-Kister expansion equation. The excess molar heats and volumes for [EMIM][MeSO4] with all the six co-solvents were all negative over the entire composition range. The size and structure of these molecules, as well as the nature of interaction between the mixture components seems to play a major role in the observed effects. These results have been interpreted in terms of hydrogen bonding and intermolecular interactions between ionic liquid ions and co-solvent molecules via radial distribution functions. The results showed that hydrogen bond interactions between ionic liquid and co-solvents play a major role in the volumetric and energetic properties of the mixtures. Further, radial distribution functions suggest that ionic liquid ions aggregate at higher concentrations of co-solvents. (C) 2018 Published by Elsevier B.V.