Fe2+ should already be present in both, though is generally designed to only be used in heme (at least for CHARMM), so your mileage may vary for other species.The paper above for CHARMM has new Fe2+ parameters but I have never used them so I cannot vouch for their quality.
Add the relevant parameters to ffnonbonded.itp and if you need pdb2gmx to create a topology with them, also add relevant information to the .rtp and .atp files.
Thank you for your reply. I am trying to build the protein-Fe complex with charmm36-mar2019.ff. In the charmm36-mar2019.ff directory, there are the atomtypes.atp, ffnonbonded.itp and merged.rtp files. These files all contains lines for FE. For example, in the ffnonbonded.itp file, there is a line reading:
FE 26 55.847000 0.000 A 0.115816833358 -0.00000
In the atomtypes.atp file, there is a line reading:
FE 55.84700 ; heme iron 56
In the merged.rtp file, there are lines reading:
[ HEME ]
[ atoms ]
FE FE 0.240 0
However, when I use pdb2gmx to build the topology, I get the following error:
Residue ‘FE’ not found in residue topology database
So Could you further explain how to deal with this situation?
Be advised that the CHARMM Fe2+ parameters are only developed for heme and have not been validated for His and Asp/Glu interactions like you have in your complex. You should do some validation against QM interaction energies and dimer geometries before trying to just apply the existing parameters. I suspect you will find you need NBFIX terms.
I added the following lines to the merged.rtp file:
[ FE ]
[ atoms ]
FE FE 2.000 0
Now the pdb2gmx has worked through. However, could you also give me some instruction on:
(1) how to validate for Fe2+ interactions with His and Asp/Glu?
(2) how to validate against QM interaction energies?
Do you mean the paper for the metal ion charmm ff you listed above all Just the Charmm ff paper?
And I wonder what is the most reliable method to calculate the binding free energy between the Fe2+ and the protein, the SMD method or the alchemical method?
Any review of the CHARMM force field will tell you all about how the force field is parametrized, including the proper QM methods.
Decoupling charged species is tricky. My first try would be with SMD, but given the partial covalent nature of the actual interactions of these species with proteins, obtaining something physically realistic may be extremely challenging.
I am trying to do MD for the above mentioned protein-Fe2+ complex. I first did minimization, which finished briefly with the following output:
Steepest Descents converged to machine precision in 40 steps,
but did not reach the requested Fmax < 10.
Potential Energy = -1.1980662e+06
Maximum force = 3.3884829e+08 on atom 2806
Norm of force = 1.8291110e+06
However, when I did the NVT step, I got the following error:
step 5: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Back Off! I just backed up step5b.pdb to ./#step5b.pdb.1#
Back Off! I just backed up step5c.pdb to ./#step5c.pdb.1#
Wrote pdb files with previous and current coordinates
Segmentation fault (core dumped)
Is there a way to tackle this problem?
Yeah, the normal magnitude of the system maximum force is 10^2, isn’t it?
The atom 2806 is the OE1 atom of Glu189, one of the contact atom of the Fe2+. So is the problem arose by the fact that the OE1 atom is too close to the FE2+ in the original coordinates? How to deal with this problem?