The purely molecular dynamical formulation of grand canonical ensemble of Pettitt and co-workers was extended to implement the Nose-Hoover thermostat and introduce multiple fractional particles. The algorithm was applied to simulation of muVT ensembles of TIP4P water and methanol molecules at 298 K. The procedure reproduced the experimental density of water for input chemical potential of -24.0 kJ/mol, as well that of methanol for -19.0 kJ/mol. It was shown that by freezing the fractional dynamics and assigning values corresponding to desired integration points to the fractional particles, one can apply the polynomial-path method of thermodynamic integration to configurations sampled during the simulation run. Thermodynamic integration is performed concurrently with the main simulation; moreover, one obtains four integration points in a single run. The resulting free energies agreed very well with model-specific values of excess free energy calculated by scaling the interactions of all particles. On the other hand, the input chemical potential showed considerable deviation from its actual value in the case of methanol; this discrepancy was blamed on deficiencies in the treatment of ideal free energy under all existing schemes.