GROMACS version: 2025.4-conda_forge
GROMACS modification: No
I am processing a zinc finger protein structure, in which Zn ions form four-coordinate bonds with two Cys SG atoms and two His NE2 atoms (Cys (SG) and His (NE2)). There are a total of 9 Zn ions, and I have identified the residue atomic positions involved in the four-coordinate interactions of all these 9 Zn ions.
I have also read several forum threads discussing methods for handling Zn coordination bonds and gained valuable insights ——
use only harmonic potential between the Zn and SG/NE2，bond type 6
However, after checking the syntax for bond type 6, I have some questions regarding how to determine the b0 and kb parameters for these bonds.
As for b0, I may calculate it based on the actual distances observed in the PDB structure. For example, for Zn-NE2 bonds, I could collect data from all 9 sets of four-coordinate bonds and compute the average bond length, then repeat this process to obtain the average values for the other three bond types. Are there any more recommended approaches for this step?
What about the kb parameter? I have not yet found relevant values in published literature. How should I set these parameters?
Additionally, I am also wondering whether I can use bond type 1. Why is bond type 6 used instead?