How to deal with Zn ion in Zinc finger protein

GROMACS version: 2025.4-conda_forge
GROMACS modification: No

I am currently conducting simulations on the binding mechanism between a zinc finger protein and a segment of DNA sequence. This zinc finger protein contains nine zinc finger domains that require nine zinc fingers for structural stabilization. My initial structure, predicted by AlphaFold3, consists of a truncated zinc finger sequence (rather than the full-length sequence), the DNA sequence, and nine zinc ions. Afterwards, I ran the most basic command for this structure:

gmx pdb2gmx -f znf263_zf1_9_hs7_c36.pdb -o znf263_zf1_9_hs7_c36.gro -water tip3p -ter -ignh

I chose NH3+ / COO- as the termini for my truncated sequence, yet the log indicates an issue with the processing of my zinc ion chain:

Processing chain 2 'B' (9 atoms, 1 residues)

Residue ZN1 has type 'Ion', assuming it is not linked into a chain.

Problem with chain definition, or missing terminal residues. This chain does not appear to contain a recognized chain molecule. If this is incorrect, you can edit residuetypes.dat to modify the behavior.
9 out of 9 lines of specbond.dat converted successfully
Opening force field file ./charmm36-feb2026_cgenff-5.0.ff/aminoacids.arn

Checking for duplicate atoms....

Now there are 1 atoms. Deleted 8 duplicates.

Generating any missing hydrogen atoms and/or adding termini.

Now there are 1 residues with 1 atoms
Chain time...

Making bonds...

No bonds

Generating angles, dihedrals and pairs...

Making cmap torsions...

There are    0 dihedrals,    0 impropers,    0 angles
             0 pairs,        0 bonds and     0 virtual sites

Total mass 65.370 a.m.u.

Total charge 2.000 e

Writing topology

I don’t know how to handle this?

I am also in similar boat my freind. Your best option is to firstly understand how your zinc ions are coordinating to your neighbouring residues/ligands. For example if your Zinc is coordinated tetrahedrally to four cysteine residues, then most likely those cysteine residues should be deprotonated, the residue identifier should also be changed from CYS to CYM(deprotonated).

Once you are clear how Zinc is coordinated in your starting structure, then you must decide how you will paramterise the zinc ions present in your structure for the MD simulation run. CharmmGUI is good alternative to generate your input files.
This is study that benchmarks bonded and non-bonded models. https:// doi.org/10.3390/ijms24065440 and https://doi.org/10.1002/jcc.27052

The above study is good starting point. Keep in mind you will have to make do with tools you have for example:

  • If you are using Gromacs as your MD engine: you cant use the special non-bonded model 12-6-4 LJ, this is only available on AMBER MD engine i believe. Stick with standard 12-6LJ model.

  • If you go down the route of paramterising your zinc ions using bonded model, MCPB.py is very good tool, but i would suggest using this tool if you have Gaussian or GAMESS-US QM software for geometery and frequency calculations. http://dx.doi.org/10.1021/acs.jcim.5b00674
    If you dont have have Gaussian license, ORCA is a free QM software but this is not natively supported in MCPB.py. So i would suggest EasyPARM: https://doi.org/10.1063/5.0301038

FYI you dont have to know everything stated above, little by little you pick up information and put the pieces together, hopefully it will make sense. This took me months to get it to work.

Dear slaa1999, thank you very much for your reply. Currently, I am stuck at an issue where nine Zn ions cannot be recognized correctly.

My structural file contains nine Zn ions, and I have extracted the relevant segment of the PDB file as follows:

Chain A at the front corresponds to the protein, chains B to J are my nine Zn ions, and chain K onwards is DNA.

ATOM   2465  OXT GLY A 306      10.724  65.388 -70.879  1.00 54.11           O  
TER    2466      GLY A 306                                                       
HETATM 2466 ZN    ZN B   1       0.577 -10.370 -14.292  1.00 81.23          ZN  
TER    2467       ZN B   1                                                       
HETATM 2467 ZN    ZN C   1       6.663 -17.401  30.574  1.00 52.72          ZN  
TER    2468       ZN C   1                                                       
HETATM 2468 ZN    ZN D   1      21.273  11.441 -31.875  1.00 70.86          ZN  
TER    2469       ZN D   1                                                       
HETATM 2469 ZN    ZN E   1     -14.967  11.850 -17.113  1.00 85.54          ZN  
TER    2470       ZN E   1                                                       
HETATM 2470 ZN    ZN F   1       7.778  26.846 -14.643  1.00 74.03          ZN  
TER    2471       ZN F   1                                                       
HETATM 2471 ZN    ZN G   1      32.403  32.856 -51.445  1.00 43.71          ZN  
TER    2472       ZN G   1                                                       
HETATM 2472 ZN    ZN H   1      17.105  27.698 -73.080  1.00 39.01          ZN  
TER    2473       ZN H   1                                                       
HETATM 2473 ZN    ZN I   1      10.810  45.327 -40.338  1.00 40.23          ZN  
TER    2474       ZN I   1                                                       
HETATM 2474 ZN    ZN J   1       5.666  52.790 -72.310  1.00 32.79          ZN  
TER    2475       ZN J   1                                                       
ATOM   2476  P    DT K   1      27.454  56.961 -76.759  1.00 63.20           P  

I did not encounter any error messages when running the pdb2gmx command, yet I found confusing content in the log file:

Residue ZN1 has type 'Ion', assuming it is not linked into a chain.

Problem with chain definition, or missing terminal residues. This chain does not appear to contain a recognized chain molecule. If this is incorrect, you can edit residuetypes.dat to modify the behavior.
9 out of 9 lines of specbond.dat converted successfully
Opening force field file ./charmm36-feb2026_cgenff-5.0.ff/aminoacids.arn
Opening force field file ./charmm36-feb2026_cgenff-5.0.ff/na.arn

Checking for duplicate atoms....

Now there are 1 atoms. Deleted 8 duplicates.

Generating any missing hydrogen atoms and/or adding termini.

Now there are 1 residues with 1 atoms
Chain time...

Back Off! I just backed up topol_Ion_chain_B.itp to ./#topol_Ion_chain_B.itp.1#

Making bonds...

No bonds

Generating angles, dihedrals and pairs...

Making cmap torsions...

There are    0 dihedrals,    0 impropers,    0 angles
             0 pairs,        0 bonds and     0 virtual sites

Total mass 65.370 a.m.u.

Total charge 2.000 e

Writing topology

Back Off! I just backed up posre_Ion_chain_B.itp to ./#posre_Ion_chain_B.itp.1#

As you can see, eight of my nine Zn ions were deleted, leaving only chain B, which means merely one Zn ion remains.

Including chain 1 in system: 4817 atoms 306 residues

Including chain 2 in system: 1 atoms 1 residues

Including chain 3 in system: 2414 atoms 75 residues

Including chain 4 in system: 2348 atoms 75 residues

I have no idea why this occurs, but I clearly need to retain all nine Zn ions in my simulation.

I would highly recommend you to use CharmmGUI for now. There are extensive videos on how to use CharmmGUI webserver. This will generate your input file.