  In a recently proposed computational model of open molecular systems out of
equilibrium [Ebrahimi Viand et al. J.Chem.Phys. 153, 101102 (2020)], 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, we test the range of applicability of
the computational model with a linear coupling to two different reservoirs,
which induces a nonequilibrium situation. To this end, we studied the density
profiles of Lennard-Jones liquids in large thermal gradients using
nonequilibrium molecular dynamics simulations with open boundaries. We put in
perspective the formulation of an extension of the mathematical model that can
account for nonlinear effects.
