How to set up topology for metalloprotein and ligand in free energy perturbation calculations

User discussions
1

GROMACS version: 2022.5
GROMACS modification: No

I am working on a system with a protein, one zinc ion, and a ligand. In the binding pocket, the zinc ion interacts with two HID residues, one HIE residue, and one nitrogen atom on the ligand. I have used MCPB.py included in the AMBER package to generate metal-center parameters for this system and converted the parameters to Gromacs format (xxx.top and xxx.gro files) for further simulations.

I am aware that I can use separate .itp files (protein.itp, ion.itp, ligand.itp) or merge them into a single complex.itp file for MD simulations, which works well. However, for further free energy perturbation (FEP) calculations, the ligand topology needs to be separate to define the ligand using couple-moltype = LIG for ligand coupling or decoupling calculations. This separation prevents me from adding bonded parameters between the ligand atom and Zn2+.

I tried to start the numbering in the ligand .itp file from where the numbering in the protein .itp file ends, like:

protein.itp:

[ moleculetype ]
; Name            nrexcl
PROTEIN           3

[ atoms ]
; nr  type  resnr  residue  atom  cgnr  charge  mass
1     N     1      PROT     N     1     0.345   14.01
2     CA    1      PROT     CA    1     0.095   12.01
; ... other protein atoms ...
1001  Zn    50     ZN       ZN    50    2.000   65.38

[ bonds ]
; ai   aj   funct   length    force.c
1     2     1       0.152     320000  ; Example protein bond

[ angles ]
; ai   aj   ak   funct   angle   force.c
1     2     3    1       109.5   450   ; Example protein angle

[ dihedrals ]
; ai   aj   ak   al   funct   dihedral   force.c
1     2     3    4    1       180.0     4     ; Example protein dihedral

ligand.itp:

[ moleculetype ]
; Name            nrexcl
LIGAND            3

[ atoms ]
; nr   type   resnr  residue  atom   cgnr  charge  mass
1002  C1     51     LIG      C1     51    -0.120  12.01
1003  O1     51     LIG      O1     51    -0.500  16.00
; ... other ligand atoms ...

[ bonds ]
1001  1002  1       0.19618   84684.0 ; Zn - Ligand bond

[ angles ]
; ai   aj   ak   funct   angle   force.c
1001  2     1002  1       120.0   300   ; Zn - Ligand angle

However, it results in an error indicating that atom numbers in bonds are out of bounds. It seems that the atom numbering must start from 1 in every individual .itp file. Given this constraint, how can I create the relationship between the ligand and the zinc ion?