We present a numerical method for solving the Langevin dynamics model. Rather than the trajectory-wise accuracy, we focus on the consistency to the equilibrium statistics at the discrete level. A discrete fluctuation-dissipation theorem is imposed to ensure that the statistical properties are preserved.