Issue in mass density fluctuation in NPT ensemble with GPU calculations

GROMACS version: 5.1.4 Single precision
Hardware: GTX 1080Ti, CPU E5 2673
GROMACS modification: No
I performed NPT run with GPUs. The density fluctuations is too large than normal.
BTW, the simulated system is very large, in which contains almost 40k atoms. And polar groups with large charges are included, such as COOH/CH2OH.
The problem is this: NPT run with the 4 GPUs (4 MPI rank) will result the “abnormal density fluctuations”, but CPU-only NPT run can get right results. It seems that the MPI GPU calculations has bugs.
The attached files are snapshot after long equilibration, while the above problem is still there.