A negative thermal conductivity would imply negative entropy
production. It's against the rules. However, in a finite time
simulation, there is no law that prevents the finite-time integral of
the average flux-flux autocorrelation function from taking on a
negative value. I suggest examining the way in which your
calculations depend on things like: equilibration time, sampling time,
sampling frequency, upper limit on the Green-Kubo time integral,
system size, timestep size, etc. Also consider getting a second
optinion using:

I started to study on RNEMD that is discussed in the link you provided.
Besides, I think I chose the optimized timestep size with respect to my
previous trials and literature, and the system is shown as in equilibrium
without a problem. But I will try different number of atoms as you advice.
Another problem is I dont know how to select optimized sampling time,
frequency and upper limit of GK integral(most of the papers dont include
this specific detail as well) and trying different numbers without an
reference doesnt seem like a feasible trial and error solution for me.