I am wondering why I am getting different simulation results when changing the skin distance from 0.3 to 2, 3, 3.5, 5…, even though I am updating the neighbor list at every time step using:
neigh_modify delay 0 every 1 check yes
Note: I am using metal units and the rebo potential.
You are NOT updating the neighbor list at every step. You are only checking if the neighbor list needs to be updated at every step. Thus the different neighbor list skin values will result in different number of steps between neighbor list updates because the check is if any particle has moved by more than half the skin distance. You can verify this by looking at the neighbor list statistics at the end of a run. To force neighbor list updating in every steps (and making your simulation significantly slower) you have to use check no. There is an optimum choice of skin since a larger skin will reduce the number of times the neighbor list needs to be updated, but will slow down the non-bonded computation due to having more neighbors and needing to compute their distances.
Also, you are not really getting “different results”. You are just getting diverging trajectories because MD is a chaotic system and thus you are subject to the so-called “butterfly effect” (scientifically referred to as Lyapunov instability).