Huge RMSD for ligand when using protein for least squares fit

GROMACS version: 2023.2
GROMACS modification: No

Dear GROMACS users,

I’ve produced an MD trajectory for a homotetrameric protein complex in which each chain is bound to individual ligands that are identical. The PBC effects have been removed from the trajectory. Upon its inspection, it is clear that all ligands remained bound to their respective binding sites throughout the simulation. However, no matter which protein groups I select for the least squares fit, I get huge RMSD values for the ligands varying between 10 and 15 nm. When I use the ligand itself for the least squares calculation, this value drops to 0.1 - 0.3 nm, which is expected.

Can anyone propose a possible explanation for this behavior and how to fix it?

Best regards,

Gustavo