I’m simulating a Cytochrome P450 enzyme using GROMACS with the CHARMM36 force field. I want to ensure that the cysteine residue (thiolate) properly forms a coordinate bond with the Fe atom in the HEME group, as is required for P450 function. The distance between SG and Fe is 0.2039nm in the hem_swiss.pdb.
Here’s what I did:
I used gmx pdb2gmx -f hem_swiss.pdb -o swiss.gro -p topol.top -ter -merge all
I confirmed that both the HEM group and the protein are present in the swiss.gro and topol.top files.
The Fe atom from HEM appears as in swiss.gro file:
- 501HEM FE 1 1.999 -0.131 -1.334
However, the cysteine residue that is supposed to ligate Fe is still labeled as CYS (not CYS2):
345CYS SG 5496 2.152 -0.243 -1.258
And no bond between Fe and SG was created (confirmed by checking [ bonds ] in the topology).
⚠️ Problem:
Because this Fe–SG bond is missing, the distance between Fe and SG drifts from ~0.2–0.3 nm up to 0.6–0.8 nm after 120 ns of simulation, meaning the coordination is lost.
❓ My question:
What exact steps or conditions are required to ensure that pdb2gmx will:
Rename CYS to CYS2
Form a bond between the Fe of HEME and the SG atom of the cysteine
Write that bond into the [ bonds ] section of the topology
I already checked that:
The specbond.dat file contains:
CYS SG 1 HEM FE 2 0.25 CYS2 HEME
The residuetypes.dat file includes:
CYS2 protein