Please see the documentation:
Since you call the function from every MPI rank, the components in the list are divided by the number of MPI ranks for the LAMMPS instance, since value will be later summed over all MPI ranks.
Since you mention box volume, are you certain that you are not mixing up virial and pressure?
I’m computing the global virial every step via compute virial all pressure NULL virial which should just be the virial component of the pressure. I don’t need to worry about mixing the ke term between states as it is constant between states.
However, this doesn’t explain my observations previously stated. The documentation states that I need to supply every MPI rank with the same virial tensor, which is what I am doing.
No, you are not. You are computing the global pressure but only include the virial component.
The virial part of the global pressure (in 3d) is the negative of the sum of the three diagonal terms divided by 3 and divided by the volume of the box.
To quote from the docs of the library interface function:
The components of the virial need to be stored in the order: xx, yy, zz, xy, xz, yz. In LAMMPS the virial is stored internally as stress*volume in units of pressure*volume as determined by the current units settings […]