Recently I've tried to use the following parameters:

https://www.mayo.edu/research/labs/computer-aided-molecular-design/projects/zinc-protein-simulations-using-cationic-dummy-atom-cada-approach

...to represent a zinc ion coordinated by four CYM residues. I'm using GROMACS with Amber force field. I've added the above-mentioned parameters to ffbonded.itp, ffnonbonded.itp and aminoacids.rtp files. pdb2gmx recognized the new residue, as well as grompp did. But then mdrun crashed. I either got LINCS warnings all the time, or strange output files appeared (with no coordinates). I changed the interaction type between dummy atoms from 1 to 6 and the system started to work quite properly. However, I still had to change timestep to 1fs (because of short oscillation period between two particular dummy atoms - but only these two, while there are four identical dummies with identical interactions). Another problem is that geometry of the artificial zinc 'molecule' tends to be disturbed when one of the cysteines is too far, with two dummies getting almost the same coordinates and the entire molecule being planar. I thought that the system will get equilibrated during NVT on its own, but it din't and I had to manually change and restrain the coordinates of dummies and SG atoms to make it work.

To sum up, I had to make a lot of manual manipulations with my zinc-containing system to make it work, and I'm afraid that something is wrong. I have a feeling that if everything was OK, the whole thing would work more automatically. I'd be grateful for suggestions and opinions from users who used this or similar approach.