GROMACS simulation using the AMBER99SB*-ILDN forcefield

GROMACS version:
GROMACS modification: Yes/No
Here post your question Hello GROMACS users,

I have a couple of questions:

  1. I am working on a protein that has two chains (identical) so is it ok to use one of them to generate the topology based on the force field? Although the interaction of interest can be observed by considering one chain (because the crystal structure of protein-ligand shows the binding of the ligand to one chain), is it ok to do that?

  2. What can happen If I want to use both the chains, as it seems like the force field is still generating one file, but two .itp files for each chain (eg for chain A: topol_protein_chain_A.itp and porse_protein_chain_A.itp) ?

  3. At last, I am also getting an error when generating topology: “Duplicate line found in or between hackblock and rtp entries” what does it mean ??
    (I did some search and it seems like if your input .pdb file has any gaps it shows this error, is it correct ?? Although the protein crystal structure itself does not have information for some missing loops, in this case how to proceed further ??

At last, I want to thank Dr. Lemkul for providing helpful tutorials on GROMACS.

Depends on the biology of the system. If it’s a functional dimer, you may need to simulate both protein chains. If it’s a crystallographic dimer only, then you can simulate one chain.

The topology and position restraint topologies are different. If you include both chains in your input coordinates, you will get topology and restraint files for both of them. pdb2gmx natively deals with complexes. The will basically be a series of #include statements and the chain topologies will be identical, but that’s what you want (and is, of course, correct).

This means there are issues in the force field files. If you need help troubleshooting that, please provide your full pdb2gmx screen output.

Hello Dr. Lemkul, Thank you for the reply. Below are the screenshots from the pdb2gmx screen output.

Also I using the generalized amber force field procedure to generate the ligand topology file (using acpype) but I also want to include the Hydrogen Mass Repartiotioning. I am aware of the importance of HMR but I am not sure if this is required for my system or not I am kind of confused about where to use it and where not to use it.

I am following a paper and trying to reproduce somehow the same system dynamics in which they have mentioned that “Hydrogen masses of the ligand were set to 4 u to allow for a timestep of 4 fs.”

pdb2gmx is just being noisy and such messages should probably be suppressed. Your topology was successfully written so there’s no actual problem.

This is what the -heavyh option is for. You need to repartition all hydrogen masses to be able to extend the time step, not just the ligand, AFAIK.

In case you have not come acroess, there is a merge option for pdb2gmx. This will make all chains in pdb file to be combined as a single chain so it wouldnt output 4 .itp files for each chain. There is also another option of selecting a specific chain in the pdb file, but i havent used this option myself.

Hello Dr. Lemkul,

Thank you for the response. You mentioned that “You need to repartition all hydrogen masses to be able to extend the time step”, with all you mean I have to repartition my protein molecule also along with the ligand and then separate the ligand from protein, following generating the topologies for each separately?

If yes, I am also curious to know that what affect will it have on my simulation if I follow the normal step (following your tutorial, and of course including the necessary changes) and generate the topologies and carry on with my simulation.

Hello @darrenphua,

Thank you for your reply. That’s a good thing to know. But now as far as the topology for the protein is concerned I believe that the protein topology provided by pdb2gmx is correct. I am now more concerned about how to use the -heavyh command to get the repartitioning correct, as I have never done that before. I will appreciate if I can get some helpful suggestions on that too. Thank you!

I have not made hydrogens heavy before but it seems like a great way to speed up computation time. Might wanna consider this for my MD runs as well!
Off the top of my head, i dont foresee much issues with making hydrogens heavy except in simulations where you need the usual movements of hydrogens such as during binding with proteins or ligands? Or maybe cases where hydrogen bonding is concerned?

@darrenphua as far as my knowledge goes I believe it would be a good idea to consider -heavyh on the whole protein-ligand complex and then separate the ligand from the protein (that might take care of the concern regarding the hydrogen bonding) I might be wrong too. But to me, that seems to be a logical way.

There’s no reason to think the ligand’s H atoms are special and can be made heavy while the rest of the system isn’t. Your fastest degrees of freedom then reside in the protein and you can’t use a 4-fs time step. Either repartition everything, or nothing.

Dr. @jalemkul.
Thank you so much for the quick response and clarification, that absolutely makes sense.
Next, as I have never done repartitioning, and all I can find is the theory regarding why it’s done and its application, with no concordant steps to be followed to attain. May I request, if I can get suggestion on how to proceed ? Thank you.

What about virtual hydrogen sites? I read that it gives less artifacts than repartitioning and can increase timestep to 5 fs as well.

@darrenphua I 'll read about it, thanks.