In a recently proposed computational model of open molecular systems out of equilibrium the action of different reservoirs enters as a linear sum into the Liouville-type evolution equations for the open system's statistics. The linearity of the coupling is common to different mathematical models of open systems and essentially relies on neglecting the feedback of the system onto the reservoir due to their interaction. In this paper, the range of applicability of the computational model is tested with a linear coupling to two different reservoirs, which induces a nonequilibrium situation. To this end, the density profiles of Lennard–Jones liquids in large thermal gradients are studied using nonequilibrium molecular dynamics simulations with open boundaries. The authors put in perspective the formulation of an extension of the mathematical model that can account for nonlinear effects.