When liquid mixtures containing 2 or more constituants are subject to thermal gradients, their respective concentrations are not homogeneous anymore and usually reach a steady-state where certain molecules accumulate in cold regions while others are rather more concentrated in hot regions, a phenomenon known as thermophoresis. Also called the Soret effect, it has been evidenced in a variety of systems and has many practical applications, as well as implications in the context of the origins of life where it could have served as way to accumulate otherwise dilute precursors in the absence of biological compartments. While it was discovered in the 19th century, a complete theoretical picture able to explain this effect is still missing and phenomenological approaches are often employed to account for the experimental observations. In particular, the amplitude of the resulting concentration-gradients (quantified by the Soret coefficient) depends on many factors and is not straightforwardly rationalized. All-atom molecular dynamics simulations appear as an exquisite tool to shed light on the molecular origins for this phenomenon in molecular systems. Although they have already been applied in such contexts in the past, they are not routinely employed, in particular in the case of dilute aqueous solutions, which are of particular relevance in the context of reactant accumulation in aqueous environments, but which pose significant practical challenges. Here, we propose a robust all-atom MD approach to tackle thermophoresis in dilute realistic solutions at the molecular level. We propose to rely on a recent enhanced heat-exchange algorithm to generate temperature-gradients, which, to the best of our knowledge, has not yet been applied in this context. We carefully assess the convergence of thermophoretic simulations in dilute aqueous solutions. We show that simulations typically need to be propagated on very long timescales (hundreds of nanoseconds) to generate reliable concentration-gradients. We find that the magnitude of the temperature gradient and the box sizes have little effect on the measured Soret coefficients. Practical guidelines are derived from such observations. Provided with this reliable setup, we discuss the results of thermophoretic simulations on several examples of molecular, neutral solutes, which we find in very good agreement with experimental measurements regarding the concentration-, mass-, and temperature-dependence of the Soret coefficient.