Isotropic deformation of the simulation box with deform option

GROMACS version: 2018.9

Hello everyone,
I want to isotropically deform the simulation box of my system in order to calculate the bulk modulus. The box is cubic and I need to deform the box at the same rate along all three axis directions. I used the following mdp options:
pcoupl = Parrinello-Rahman
pcoupltype = anisotropic
tau_p = 1.0
ref_p = 0.0 0.0 0.0 0.0 0.0 0.0
compressibility = 0 0 0 4.5e-5 4.5e-5 4.5e-5
refcoord_scaling = com
deform = 0.00001 0.00001 0.00001 0 0 0
The temperature coupling with V-rescale is also used. But I found that at the very beginning of the simulation, the box will have off-diagonal value. It is no more a cubic box and it is not expected. Could you please tell me why the box will deform along XY, YZ and ZX? And is that possible to calculate the bulk modulus with the deform option?
Any help would be much appreciated!