I have been running MD simulation with respect to free surface effect of metal.

During the series of simulation with ascending thickness of material in the direction of free surface , there is fluctuation of virial stress on surface atoms.

Can you give me the answer of question whether the fluctuation of virial stress on surface atoms is numerical error of LAMMPS?

in my experience that’s an effect of the virial formulation, it performs well in situations like steady state averages. It’s worth noting that any quantity in the phase space of the simulation has a sinusoidal error amplitude that scales with the timestep. That’s what stable error looks like.

Consider a single pair of atoms, 1 and 2, in vacuum, interacting via a harmonic spring. The virial stress tensor is proportional to F_{12} x R_{12}. At any instant, this can be positive or negative, depending on whether the spring is in tension or compression. We can only expect that the average value will be close to zero, and even then, there is lots of stat mech fine print, but for large systems, it is true that the non-periodic components of the virial will be approximately zero (when the kinetic contribution is added in). As you average over larger and larger surface areas A, and/or longer and longer times T, the magnitude of the fluctuations in the virial should decrease roughly as 1/(AxT)