A new hybrid molecular dynamics-hydrodynamics method based on the analogy with two-phase flows is implemented that takes into account the feedback of molecular dynamics on hydrodynamics consistently. The consistency is achieved by deriving a discrete system of fluctuating hydrodynamic equations whose solution converges to the locally averaged molecular dynamics field exactly in terms of the locally averaged fields. The new equations can be viewed as a generalisation of the classical continuum Landau-Lifshitz fluctuating hydrodynamics model in statistical mechanics to include a smooth transition from large-scale continuum hydrodynamics that obeys a Gaussian statistics to all-atom molecular dynamics. Similar to the classical Landau-Lifshitz fluctuating hydrodynamics model, the suggested generalised Landau-Lifshitz fluctuating hydrodynamics equations are too complex for analytical solution; hence, a computational scheme for solving these equations is suggested. The scheme is implemented in a popular open-source molecular dynamics code GROMACS (GROningen MAchine for Chemical Simulations), and numerical examples are provided for liquid argon simulations under equilibrium conditions and under macroscopic flow effects.