2 particles communicated to PME rank 23 are more than 2/3 times the cut-off out

GROMACS version: 2020.4
GROMACS modification: Yes, with PLUMED
I know there are already a few posts about this error, but no clear solutions have been provided (or at least I didn’t find it).

In my case, I was running an NPT MD simulation for a glycosylated insulin, but I got the following fatal error:

Fatal error:
2 particles communicated to PME rank 23 are more than 2/3 times the cut-off
out of the domain decomposition cell of their charge group in dimension y.
This usually means that your system is not well equilibrated.

Before the MD simulation, the structure had been through energy minimization, 200 ps NVT equilibration and then 200 ps NPT equilibration without any error. Notably, the error happened at 209 ns and there were a few LINCS warnings as below:

Step 104531392, time 209063 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.638738, max 2.927066 (between atoms 12 and 13)
Step 104531392, time 209063 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.471350, max 2.267098 (between atoms 72 and 74)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
72 73 114.6 0.9406 0.2428 0.1010
72 74 107.8 0.7914 0.3300 0.1010
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
12 13 104.8 14.3085 0.4281 0.1090

It is surprising to me that the system appears unstable after such a long simulation. However, I’m not sure how to solve this problem. It would be nice if someone could provide some guidance.

For your information, here I attach the mdp files that I was using.

md.mdp (2.3 KB)
npt_equil.mdp (2.3 KB)
nvt_equil.mdp (2.0 KB)