We implement the statistically sound G-JF thermostat for Langevin Dynamics simulations into the ESPREesSo molecular package for large-scale simulations of soft matter systems. The implemented integration method is tested against the integrator currently used by the molecular package in simulations of a fluid bilayer membrane. While the latter exhibits deviations in the sampling statistics that increase with the integration time step dt, the former reproduces near-correct configurational statistics for all dt within the stability range of the simulations. We conclude that, with very modest revisions to existing codes, one can significantly improve the performance of statistical sampling using Langevin thermostats.