Please, always reply to all to include the list.
Assuming your equations are correct you got to have a problem somewhere as the crash you observe while loosing atoms is not normal. I can understand that 1fs might be a bit too large to sample properly the dynamics of certain bond types including oxygen atoms. But I guess typically-small oscillation periods are about 9-10fs (for O-H for example in OPLS) which results in a timestep of 0.4-0.5fs. Now, you are using 0.08fs as timestep which is pretty small. I find hard to believe that 1000000 steps with such a timestep provides you with enough statistics about the steady state (note: I only have experience with this kind of simulations working on liquids not solids).
The other reason that I can imagine could be affecting your dynamics will be that swapping velocities in the solid could have an effect on the bond dynamics and thus your need for a much small timestep.
How do the RNMED and NEMD results compare?
For NEMD I never use fix heat but just two Langevin thermostats. They are easy to set up and you know right away what the temperature gradient should be.
Carlos