Zinc displacement and coordination

GROMACS version: 2023.1
GROMACS modification: Yes/No

Dear users,
I have to deal with a simulation of a protein containing two zinc ions, tetrahedrally coordinated but with two different sets of ligands. In the first case, one zinc is bound by four cysteines residues - and I guess all of them should be in anioni form - whereas in the second case, one zinc is coordinated by two cysteines and two histidines. I tried different Force Fields, but, actually, any attempts to perform simulation was unsuccessful, as the two zinc are displaced during NVT equilibration.
I realized - but I am not completely sure - that the first problem could be due the fact that topology build up renders only cysteines as -CH2-SH and not S- (minus, as anion).
Do you have any experience or any hints about the solution of this problem?
thanks in advance

I was wondering if someone does use H++ server, as, I guess, It May address the aforementioned issue about cysteines state