shear modulus of single crystal with different crystal orientations


I'm trying to calculate the shear modulus for a single crystal using the lammps version 15May15. Therefore I deform the simulation box using the command "change_box all xy delta ${delta} remap units box" and let the stress per atom be calculated.
For a crystal orientation in the (100) directions it works fine and gives the value stated for the potential I'm using. But when I try to do the same with an orientation other than that, the resulting shear modulus differs from the one theoretically expected for the given orientation. Depending on the chosen orientation the difference is up the 20%.

So I wonder if any of you has an idea what the problem might be, and how to fix it.

Thank you in advance,
Lisa Neier